Documentation Center

  • Trials
  • Product Updates

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.

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® path if the Mapping Toolbox™ 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)
info = 

            Filename: '9129CATD.DDF'
               Title: 'MOUNT WASHINGTON, NH - 24000'
           ProfileID: 'SDTS RASTER PROFILE'
      ProfileVersion: 'DRAFT VERSION JULY 1997'
             MapDate: ''
    DataCreationDate: '19980811'
     HorizontalDatum: 'North American 1927'
        MapRefSystem: 'UTM'
          ZoneNumber: 19
         XResolution: 30
         YResolution: 30
        NumberOfRows: 472
        NumberOfCols: 345
     HorizontalUnits: 'METERS'
       VerticalUnits: 'METERS'
        MinElevation: 367
        MaxElevation: 1909

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))
R = 

  MapCellsReference with properties:

            XWorldLimits: [310380.84375 320730.84375]
            YWorldLimits: [4901891.5 4916051.5]
              RasterSize: [472 345]
    RasterInterpretation: 'cells'
        ColumnsStartFrom: 'north'
           RowsStartFrom: 'west'
      CellExtentInWorldX: 30
      CellExtentInWorldY: 30
    RasterExtentInWorldX: 10350
    RasterExtentInWorldY: 14160
        XIntrinsicLimits: [0.5 345.5]
        YIntrinsicLimits: [0.5 472.5]
      TransformationType: 'rectilinear'
    CoordinateSystemType: 'planar'


Step 2: Assign a Reference Ellipsoid

The value of

info.HorizontalDatum
ans =

North American 1927

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)
clarke66 = 

referenceEllipsoid with defining properties:

                 Code: 7008
                 Name: 'Clarke 1866'
           LengthUnit: 'meter'
        SemimajorAxis: 6378206.4
        SemiminorAxis: 6356583.8
    InverseFlattening: 294.978698213898
         Eccentricity: 0.0822718542230038

  and additional properties:

    Flattening
    ThirdFlattening
    MeanRadius
    SurfaceArea
    Volume

Step 3: Determine which UTM Zone to Use and Construct a Map Axes

From the MapRefSystem field in the SDTS info struct,

info.MapRefSystem
ans =

UTM

you can tell that the DEM is gridded in a Universal Transverse Mercator (UTM) coordinate system.

The 'ZoneNumber' field

info.ZoneNumber
ans =

    19

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]
provisionalZone =

19N

Then, construct a UTM axes based on this provisional zone

figure('Color','white')
ax = 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)
xCenter =

              315555.84375


yCenter =

                 4908971.5

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)
latCenter =

          44.3125403747455


lonCenter =

         -71.3126367639007

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)
actualZone =

19T

Finally, use the result to reset the zone of the axes constructed earlier.

setm(ax, '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(ax, 'MapLatLimit', [44.2 44.4], 'MapLonLimit', [-71.4 -71.2])
setm(ax, 'MLabelLocation', 0.05, 'MLabelRound', -2)
setm(ax, 'PLabelLocation', 0.05, 'PLabelRound', -2)
setm(ax, '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)
dLat =

      0.000269984934404749

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])'
xCorners =

              310380.84375
              310380.84375
              320730.84375
              320730.84375


yCorners =

                 4901891.5
                 4916051.5
                 4916051.5
                 4901891.5

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');
latCorners =

          44.2475212269387
          44.3748952132925
          44.3775277222457
          44.2501421324565


lonCorners =

          -71.374900167689
         -71.3800449718219
         -71.2502372050341
         -71.2453724066789

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;
latSouth =

                   44.2475


lonWest =

                 -71.38025


latNorth =

                  44.37775


lonEast =

                 -71.24525

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]
W =

                   0.00025                         0                -71.380125
                         0                   0.00025                 44.247625

nRows = round(   (latNorth - latSouth)     / gridSpacing)
nCols = round(wrapTo360(lonEast - lonWest) / gridSpacing)
nRows =

   521


nCols =

   540

Rlatlon = georasterref(W, [nRows nCols], 'cells')
Rlatlon = 

  GeographicCellsReference with properties:

             LatitudeLimits: [44.2475 44.37775]
            LongitudeLimits: [-71.38025 -71.24525]
                 RasterSize: [521 540]
       RasterInterpretation: 'cells'
           ColumnsStartFrom: 'south'
              RowsStartFrom: 'west'
       CellExtentInLatitude: 1/4000
      CellExtentInLongitude: 1/4000
     RasterExtentInLatitude: 0.13025
    RasterExtentInLongitude: 0.135
           XIntrinsicLimits: [0.5 540.5]
           YIntrinsicLimits: [0.5 521.5]
       CoordinateSystemType: 'geographic'
                  AngleUnit: 'degree'


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.

figure(get(ax,'Parent'))
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
Was this topic helpful?