grdcut

 GRDCUT:  reads subregion of GMT grid file into Matlab arrays x,y,z

   Usage: [x,y,z]=grdcut(grdfile, region);
    where 
          x = east coordinate vector (eg. longitude)
          y = north coordinate vector (eg. latitude)
          z = matrix of gridded values (eg. bathy grid)

          grdfile = gmt filename.
          region = [west east south north] 

   Example:
           [x,y,z]=grdcut('foo.grd');
           contour(x,y,z)

 joevans@usgs.gov

0001 function [x,y,z]=grdcut(grdfile, region)
0002 % GRDCUT:  reads subregion of GMT grid file into Matlab arrays x,y,z
0003 %
0004 %   Usage: [x,y,z]=grdcut(grdfile, region);
0005 %    where
0006 %          x = east coordinate vector (eg. longitude)
0007 %          y = north coordinate vector (eg. latitude)
0008 %          z = matrix of gridded values (eg. bathy grid)
0009 %
0010 %          grdfile = gmt filename.
0011 %          region = [west east south north]
0012 %
0013 %   Example:
0014 %           [x,y,z]=grdcut('foo.grd');
0015 %           contour(x,y,z)
0016 %
0017 % joevans@usgs.gov
0018   
0019 if ( nargin ~= 2 )
0020     help grdcut;
0021     return;
0022 end
0023 
0024 %
0025 % Error checking on the region vector.
0026 if ( region(1) >= region(2) )
0027     fprintf ( 2, 'West longitude must be numerically smaller than east longitude.\n' );
0028     help grdcut;
0029     return;
0030 end
0031 if ( region(3) >= region(4) )
0032     fprintf ( 2, 'South longitude must be numerically smaller than north longitude.\n' );
0033     help grdcut;
0034     return;
0035 end
0036 
0037 west_longitude = region(1);
0038 east_longitude = region(2);
0039 south_latitude = region(3);
0040 north_latitude = region(4);
0041 
0042 cdfid=mexcdf('open',grdfile,'nowrite');
0043 
0044 oldopts=mexcdf('setopts',0);
0045 
0046 x_range=mexcdf('varget',cdfid,'x_range',0,2);
0047 y_range=mexcdf('varget',cdfid,'y_range',0,2);
0048 spacing=mexcdf('varget',cdfid,'spacing',0,2);
0049 
0050 dims=mexcdf('varget',cdfid,'dimension',0,2);
0051 
0052 nx=dims(1);
0053 ny=dims(2);
0054 
0055 x=x_range(1)+[0:(nx-1)]*spacing(1);
0056 y=y_range(1)+[0:(ny-1)]*spacing(2);
0057 
0058 
0059 %
0060 % Now restrict the x vector to the longitudes given in the
0061 % region vector.
0062 xinds = find ( (x>=west_longitude) & (x<=east_longitude) );
0063 x = x(xinds);
0064 yinds = find ( (y>=south_latitude) & (y<=north_latitude) );
0065 y = y(yinds);
0066 
0067 z = zeros(length(yinds),length(xinds));
0068 n = length(yinds);
0069 for i = 1:n
0070     z_row = mexcdf ( 'varget', cdfid, 'z', (ny - (yinds(i)-1))*nx, nx );
0071     z(i,:) = z_row(xinds)';
0072 end
0073 %z = flipud(z);
0074 
0075 
0076 mexcdf('close',cdfid);
0077 
0078