ncdisp('geo_em.d01.nc')

clear all
close all

hgt = ncread('geo_em.d01.nc','HGT_M');
lat = ncread('geo_em.d01.nc','XLAT_M');
lon = ncread('geo_em.d01.nc','XLONG_M');

lon  = double(lon(:,:,1));
lat  = double(lat(:,:,1));

x_min = 100;    
x_max = 800;
y_min = 275;
y_max = 1010;

states = shaperead('usastatelo', 'UseGeoCoords', true)

figure(1)
ax=usamap([lat(x_min,y_min),lat(x_max,y_max)],[lon(x_min,y_max),lon(x_max,y_min)])
surfm(lat(x_min:x_max,y_min:y_max),lon(x_min:x_max,y_min:y_max),hgt(x_min:x_max,y_min:y_max))
geoshow(states,'FaceColor','none')

bounds = [lat(x_min,y_min),lat(x_max,y_max),lon(x_min,y_max),lon(x_max,y_min)]