read_srtm30plus

 READ_SRTM30PLUS  Read SRTM30+ ~1km) World topo/bathy into Matlab via WMS
 Usage:   [lon,lat,depth]=read_srtm30plus(lon_range,lat_range,res,iplot);
  Inputs: lon_range=[lon_min lon_max] (decimal degrees, east positive)
          lat_range=[lat_min lat_max] (decimal degrees, north positive)
          res = optional resolution of requested grid (arc seconds) [default=30]
          iplot = option to control plot. (1 for plot) [default=0]
          
  Outputs:lon = vector of longitudes
          lat = vector of latitudes
          z  = topography/bathymetry in meters (up positive)
               NOTE: SRTM30+ data are integers, therefore flat shallow
               areas will have a stair-step behavior, and gentle sloping
               coastal areas will have poorly defined coastlines (z=0).
          url = the url sent to the WMS

  Examples: [lon,lat,z]=read_srtm30plus([10 20],[40 46],30,1); 
               Reads and plots 30" topo/bathy for the Adriatic Sea

            [lon,lat,z,url]=read_srtm30plus([-71 -65],[40 46],60);
               Reads topo/bathy averaged to 60" for the Gulf of Maine (no
             plot)

0001 function [lon,lat,z,url]=read_srtm30plus(lon_range,lat_range,res,iplot);
0002 % READ_SRTM30PLUS  Read SRTM30+ ~1km) World topo/bathy into Matlab via WMS
0003 % Usage:   [lon,lat,depth]=read_srtm30plus(lon_range,lat_range,res,iplot);
0004 %  Inputs: lon_range=[lon_min lon_max] (decimal degrees, east positive)
0005 %          lat_range=[lat_min lat_max] (decimal degrees, north positive)
0006 %          res = optional resolution of requested grid (arc seconds) [default=30]
0007 %          iplot = option to control plot. (1 for plot) [default=0]
0008 %
0009 %  Outputs:lon = vector of longitudes
0010 %          lat = vector of latitudes
0011 %          z  = topography/bathymetry in meters (up positive)
0012 %               NOTE: SRTM30+ data are integers, therefore flat shallow
0013 %               areas will have a stair-step behavior, and gentle sloping
0014 %               coastal areas will have poorly defined coastlines (z=0).
0015 %          url = the url sent to the WMS
0016 %
0017 %  Examples: [lon,lat,z]=read_srtm30plus([10 20],[40 46],30,1);
0018 %               Reads and plots 30" topo/bathy for the Adriatic Sea
0019 %
0020 %            [lon,lat,z,url]=read_srtm30plus([-71 -65],[40 46],60);
0021 %               Reads topo/bathy averaged to 60" for the Gulf of Maine (no
0022 %             plot)
0023 
0024 % Full SRTM30PLUS info:  http://topex.ucsd.edu/WWW_html/srtm30_plus.html
0025 
0026 % Rich Signell (rsignell@usgs.gov)
0027 % Revised  Thu Jan  5 09:38:28 EST 2006
0028 %   - remove unneccesary "ind=find(z>32768)" as returned values are already short integer.
0029 %   - Leave z values as short integer in matlab instead of converting to double to save memory.
0030 if(nargin<3);iplot=0;res=30;end
0031 if(nargin<4);iplot=0;end
0032 if(nargin>4);help read_srtm30plus;end
0033 ds=res/3600;
0034 lon_range=round(lon_range/ds)*ds; 
0035 lat_range=round(lat_range/ds)*ds; 
0036 lon=(lon_range(1)+ds/2):ds:(lon_range(2)-ds/2);  % 30 sec resolution
0037 lat=(lat_range(1)+ds/2):ds:(lat_range(2)-ds/2);  
0038 nx=length(lon);
0039 ny=length(lat);
0040 %wms_server='http://maps.gdal.org/cgi-bin/mapserv_dem?';
0041 %wms_layer='srtmplus_raw';
0042 wms_server='http://onearth.jpl.nasa.gov/wms.cgi?';
0043 wms_layer='srtmplus';
0044 wms_proj='EPSG:4326';
0045 wms_format='image/geotiff&version=1.1.1';
0046 wms_max_width=4096;  %maximum image width allowed (WMS server dependent)
0047 wms_max_height=4096; %maximum image height allowed (WMS server dependent)
0048 if(nx>wms_max_width);
0049     sprintf('Maximum image width %d exceeded.  Please try again.',wms_max_width)
0050     lon=nan;lat=nan;z=nan;
0051     return
0052 end
0053 if(ny>wms_max_height);
0054     sprintf('Maximum image height %d exceeded.  Please try again.',wms_max_height)
0055     lon=nan;lat=nan;z=nan;
0056     return
0057 end
0058 % The WMS query specifies a lon,lat range and also the desired size of the image in
0059 % pixels.  We request an image size that gives us the original 30 second data values
0060 url=sprintf('%srequest=GetMap&styles=short_int&layers=%s&bbox=%d,%d,%d,%d&SRS=%s&TRANSPARENT=FALSE&width=%d&height=%d&format=%s',...
0061     wms_server,wms_layer,lon_range(1),lat_range(1),lon_range(2),lat_range(2),wms_proj,nx,ny,wms_format)
0062 pix=imread(url,'tif');
0063 z=flipud(pix);
0064 if iplot
0065   clf
0066   imagesc(lon,flipud(lat),z);set(gca,'ydir','normal');colorbar
0067   set(gca,'tickdir','out');
0068   xfac=cos(mean(lat(:))*pi/180);set(gca,'DataAspectRatio', [1 xfac 1] );
0069   if(max(z(:))*min(z(:))<0),
0070     hold on;
0071     contour(lon,lat,z,[0 0],'k-');
0072     hold off;
0073     title('SRTM30+ World Topo/Bathy');
0074   end
0075 end