www.gusucode.com > map 案例源码 matlab代码程序 > map/mapexunprojectdem.m
%% Un-Projecting a Digital Elevation Model (DEM) % % This example shows how to how to convert a USGS DEM into a regular % latitude-longitude grid having comparable spatial resolution. U.S. % Geological Survey (USGS) 30-meter Digital Elevation Models (DEMs) are % regular grids (raster data) that use the UTM coordinate system. Using % such DEMs in applications may require reprojecting and resampling them. % You can readily apply the approach show here to projected map coordinate % systems other than UTM and to other DEMs and most types of regular data % grids. % Copyright 2005-2014 The MathWorks, Inc. %% Step 1: Import the DEM and its Metadata % % This example uses a USGS DEM for a quadrangle 7.5-arc-minutes square % located in the White Mountains of New Hampshire, USA. The data set is % stored in the Spatial Data Transfer Standard (STDS) format and is % located in the folder fullfile(matlabroot,'toolbox','map','mapdata','sdts'); %% % This folder is on the MATLAB(R) path if the Mapping Toolbox(TM) is installed, % so it suffices to refer to the data set by filename alone. stdsfilename = '9129CATD.ddf'; %% % You can use the |sdtsinfo| command to obtain basic metadata about the % DEM. info = sdtsinfo(stdsfilename) %% % and you can use |sdtsdemread| to import the DEM into a 2-D MATLAB array % (|Z|), along with its referencing matrix (|refmat|), a 3-by-2 matrix that % maps the row and column subscripts of |Z| to map x and y (UTM "eastings" % and "northings" in this case). [Z,refmat] = sdtsdemread(stdsfilename); %% % You can convert |refmat| to a map raster reference object, which provides % a more complete and self-documenting way of encoding spatial referencing % information. currentFormat = get(0,'format'); format long g R = refmatToMapRasterReference(refmat, size(Z)) %% Step 2: Assign a Reference Ellipsoid % % The value of info.HorizontalDatum %% % indicates use of the North American Datum of 1927. The Clarke 1866 % ellipsoid is the standard reference ellipsoid for this datum. ellipsoidName = 'clarke66'; %% % You can also check the value of the |HorizontalUnits| field, mapUnits = info.HorizontalUnits; %% % which indicates that the horizontal coordinates of the DEM are in units % of meters, and use both pieces of information to construct a % |referenceEllipsoid|. clarke66 = referenceEllipsoid(ellipsoidName, mapUnits) %% Step 3: Determine which UTM Zone to Use and Construct a Map Axes % % From the |MapRefSystem| field in the SDTS info struct, info.MapRefSystem %% % you can tell that the DEM is gridded in a Universal Transverse Mercator % (UTM) coordinate system. % % The 'ZoneNumber' field info.ZoneNumber %% % indicates which longitudinal UTM zone was used. The Mapping Toolbox % |utm| function, however, also requires a latitudinal zone; this is not % provided in the metadata, but you can derive it from the referencing % matrix and grid dimensions. % % UTM comprises 60 longitudinal zones each spanning 6 degrees of longitude % and 20 latitudinal zones ranging from 80 degrees South to 84 degrees % North. Longitudinal zones are designated by numbers ranging from 1 to % 60. Latitudinal zones are designated by letters ranging from C to X % (omitting I and O). In a given hemisphere (Southern or Northern), all % the latitudinal zones occupy a shared coordinate system. Aside from % determining the hemisphere, the toolbox merely uses latitudinal zone to % help set the default map limits. % % So, you can start by using the first latitudinal zone in the Northern % Hemisphere, zone N (for latitudes between zero and 8 degrees North) as a % provisional zone. longitudinalZone = sprintf('%d',info.ZoneNumber); provisionalLatitudinalZone = 'N'; provisionalZone = [longitudinalZone provisionalLatitudinalZone] %% % Then, construct a UTM axes based on this provisional zone figure axesm('mapprojection','utm','zone',provisionalZone,'geoid',clarke66) axis off gridm mlabel on plabel on framem on %% % To find the actual zone, you can locate the center of the DEM in UTM % coordinates, [M,N] = size(Z); xCenterIntrinsic = (1 + N)/2; yCenterIntrinsic = (1 + M)/2; [xCenter, yCenter] = intrinsicToWorld(R,xCenterIntrinsic,yCenterIntrinsic) %% % then convert latitude-longitude, taking advantage of the fact that % xCenter and yCenter will be the same in zone 19N as they are into the % actual zone. [latCenter, lonCenter] = minvtran(xCenter,yCenter) %% % Then, with the |utmzone| function, you can use the latitude-longitude % coordinates to determine the actual UTM zone for the DEM. actualZone = utmzone(latCenter,lonCenter) %% % Finally, use the result to reset the zone of the axes constructed % earlier. setm(gca,'zone',actualZone) %% % Note: if you can visually place the approximately location of New % Hampshire on a world map, then you can confirm this result with the % |utmzoneui| GUI. % % utmzoneui(actualZone) %% Step 4: Display the Original DEM on the Map Axes % % Use |mapshow| (rather than |geoshow| or |meshm|) to display the DEM on % the map axes because the data are gridded in map (x-y) coordinates. mapshow(Z,R,'DisplayType','texturemap') demcmap(Z) %% % The DEM covers such a small part of this map that it may be hard to see % (look between 44 and 44 degrees North and 72 and 71 degrees West), % because the map limits are set to cover the entire UTM zone. You can % reset them (as well as the map grid and label parameters) to get a closer % look. setm(gca,'MapLatLimit',[44.2 44.4],'MapLonLimit',[-71.4 -71.2]) setm(gca,'MLabelLocation',0.05,'MLabelRound',-2) setm(gca,'PLabelLocation',0.05,'PLabelRound',-2) setm(gca,'PLineLocation',0.025,'MLineLocation',0.025) %% % When it is viewed at this larger scale, narrow wedge-shaped areas of % uniform color appear along the edge of the grid. These are places where % |Z| contains the value NaN, which indicates the absence of actual data. % By default they receive the first color in the color table, which in this % case is dark green. These null-data areas arise because although the DEM % is gridded in UTM coordinates, its data limits are defined by a % latitude-longitude quadrangle. The narrow angle of each wedge % corresponds to the non-zero "grid declination" of the UTM coordinate % system in this part of the zone. (Lines of constant x run precisely % north-south only along the central meridian of the zone. Elsewhere, they % follow a slight angle relative to the local meridians.) %% Step 5: Define the Output Latitude-Longitude Grid % % The next step is to define a regularly-spaced set of grid points in % latitude-longitude that covers the area within the DEM at about % the same spatial resolution as the original data set. % % First, you need to determine how the latitude changes between rows in the % input DEM (i.e., by moving northward by 30 meters). rng = info.YResolution; % In meters, consistent with clarke66 latcrad = deg2rad(latCenter); % latCenter in radians % Change in latitude, in degrees dLat = rad2deg(meridianfwd(latcrad,rng,clarke66) - latcrad) %% % The actual spacing can be rounded slightly to define the grid spacing to % be used for the output (latitude-longitude) grid. gridSpacing = 1/4000; % In other words, 4000 samples per degree %% % To set the extent of the output (latitude-longitude) grid, start by % finding the corners of the DEM in UTM map coordinates. xCorners = R.XWorldLimits([1 1 2 2])' yCorners = R.YWorldLimits([1 2 2 1])' %% % Then convert the corners to latitude-longitude. Display the % latitude-longitude corners on the map (via the UTM projection) to check % that the results are reasonable. [latCorners, lonCorners] = minvtran(xCorners, yCorners) hCorners = geoshow(latCorners,lonCorners,'DisplayType','polygon',... 'FaceColor','none','EdgeColor','red'); %% % Next, round outward to define an output latitude-longitude quadrangle % that fully encloses the DEM and aligns with multiples of the grid % spacing. latSouth = gridSpacing * floor(min(latCorners)/gridSpacing) lonWest = gridSpacing * floor(min(lonCorners)/gridSpacing) latNorth = gridSpacing * ceil( max(latCorners)/gridSpacing) lonEast = gridSpacing * ceil( max(lonCorners)/gridSpacing) qlatlim = [latSouth latNorth]; qlonlim = [lonWest lonEast]; dlat = 100*gridSpacing; dlon = 100*gridSpacing; [latquad, lonquad] = outlinegeoquad(qlatlim, qlonlim, dlat, dlon); hquad = geoshow(latquad, lonquad, ... 'DisplayType','polygon','FaceColor','none','EdgeColor','blue'); snapnow; %% % Finally, construct a geographic raster referencing object for the output % grid. It supports transformations between latitude-longitude and the row % and column subscripts. In this case, use of a world file matrix, W, % enables exact specification of the grid spacing. W = [gridSpacing 0 lonWest + gridSpacing/2; ... 0 gridSpacing latSouth + gridSpacing/2] %% nRows = round((latNorth - latSouth) / gridSpacing) nCols = round(wrapTo360(lonEast - lonWest) / gridSpacing) %% Rlatlon = georasterref(W,[nRows nCols],'cells') %% % |Rlatlon| fully defines the number and location of each cell in the % output grid. %% Step 6: Map Each Output Grid Point Location to UTM X-Y % % Finally, you're ready to make use of the map projection, applying it to % the location of each point in the output grid. First compute the % latitudes and longitudes of those points, stored in 2-D arrays. [rows,cols] = ndgrid(1:nRows, 1:nCols); [lat,lon] = intrinsicToGeographic(Rlatlon,cols,rows); %% % Then apply the projection to each latitude-longitude pair, arrays of UTM % x-y locations having the same shape and size as the latitude-longitude % arrays. [XI,YI] = mfwdtran(lat,lon); %% % At this point, |XI(i,j)| and |YI(i,j)| specify the UTM coordinate of the % grid point corresponding to the i-th row and j-th column of the output % grid. %% Step 7: Resample the Original DEM % % The final step is to use the MATLAB |interp2| function to perform % bilinear resampling. % % At this stage, the value of projecting from the latitude-longitude grid % into the UTM map coordinate system becomes evident: it means that the % resampling can take place in the regular X-Y grid, making |interp2| % applicable. The reverse approach, un-projecting each (X,Y) point into % latitude-longitude, might seem more intuitive but it would result in an % irregular array of points to be interpolated -- a much harder task, % requiring use of the far more costly |griddata| function or some rough % equivalent. [rows,cols] = ndgrid(1:M,1:N); [X,Y] = intrinsicToWorld(R,cols,rows); method = 'bilinear'; extrapval = NaN; Zlatlon = interp2(X,Y,Z,XI,YI,method,extrapval); %% % View the result in the projected axes using |geoshow|, which will % re-project it on the fly. Notice that it fills the blue rectangle, which % is aligned with lines of latitude and longitude. (In contrast, the red % rectangle, which outlines the original DEM, aligns with UTM x and y.) % Also notice NaN-filled regions along the edges of the grid. The % boundaries of these regions appear slightly jagged, at the level of a % single grid spacing, due to round-off effects during interpolation. % Move the red quadrilateral and blue quadrangle to the top, to ensure that % they are not hidden by the raster display. geoshow(Zlatlon,Rlatlon,'DisplayType','texturemap') uistack([hCorners hquad],'top') %% format(currentFormat) %% Credits % % 9129CATD.ddf (and supporting files): % % United States Geological Survey (USGS) 7.5-minute Digital Elevation % Model (DEM) in Spatial Data Transfer Standard (SDTS) format for the % Mt. Washington quadrangle, with elevation in meters. % http://edc.usgs.gov/products/elevation/dem.html % % For more information, run: % % >> type 9129.txt %