Шаг 1: Импортируйте DEM и его Метаданные

Этот пример использует USGS DEM в квадрате 7.5 минут дуги четырехугольника, расположенном в Уайт-Маунтинсе Нью-Гэмпшира, США. Набор данных хранится в Пространственном Стандарте Передачи данных (STDS) формат и расположен в папке

fullfile(matlabroot,'toolbox','map','mapdata','sdts');

Эта папка находится на пути MATLAB®, если Mapping Toolbox™ установлен, таким образом, это достаточно, чтобы относиться к набору данных одним только именем файла.

stdsfilename = '9129CATD.ddf';

Можно использовать sdtsinfo команда, чтобы получить основные метаданные о DEM.

info = sdtsinfo(stdsfilename)

info = 

  struct with fields:

            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

и можно использовать sdtsdemread импортировать DEM в 2D массив MATLAB (Z), наряду с его матрицей привязки (refmat), 3-на-2 матрица, которая сопоставляет индексы строки и столбца Z сопоставлять X и Y (UTM "движения на восток" и "northings" в этом случае).

[Z,refmat] = sdtsdemread(stdsfilename);

Можно преобразовать refmat к объекту растровой привязки карты, который обеспечивает более полный и самодокументирующий способ закодировать пространственную информацию о ссылке.

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'

Шаг 2: присвойте ссылочный эллипсоид

Значение

info.HorizontalDatum

ans =

    'North American 1927'

указывает на использование североамериканской Данной величины 1 927. Кларк 1 866 эллипсоидов является стандартным ссылочным эллипсоидом для этой данной величины.

ellipsoidName = 'clarke66';

Можно также проверять значение HorizontalUnits поле ,

mapUnits = info.HorizontalUnits;

который указывает, что горизонтальные координаты DEM находятся в модулях метров и используют оба данные, чтобы создать 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

Шаг 3: Определите который Зона UTM Использовать и Создать Карту Оси

От MapRefSystem поле в информационной структуре SDTS,

info.MapRefSystem

ans =

    'UTM'

можно сказать, что DEM с координатной сеткой в системе координат Universal, поперечной меркаторской (UTM).

ZoneNumber поле

info.ZoneNumber

ans =

    19

указывает, какая продольная зона UTM использовалась. Mapping Toolbox utm функция, однако, также требует широтной зоны; это не обеспечивается в метаданных, но можно вывести их из размерностей сетки и матрицы привязки.

UTM включает 60 продольных зон каждый охват 6 градусов долготы и 20 широтных зон в пределах от 80 градусов на юг к 84 градусам на север. Продольные зоны определяются числами в пределах от 1 - 60. Широтные зоны определяются буквами в пределах от C к X (исключение I и O). В данном полушарии (южный или Северный), все широтные зоны занимают разделяемую систему координат. Кроме определения полушария, тулбокс просто использует широтную зону, чтобы помочь установить пределы карты по умолчанию.

Так, можно запустить при помощи первой широтной зоны в северном полушарии, зоны N (для широт между нулем и 8 градусами на север) как временная зона.

longitudinalZone = sprintf('%d',info.ZoneNumber);
provisionalLatitudinalZone = 'N';
provisionalZone = [longitudinalZone provisionalLatitudinalZone]

provisionalZone =

    '19N'

Затем создайте оси UTM на основе этой временной зоны

figure
axesm('mapprojection','utm','zone',provisionalZone,'geoid',clarke66)
axis off
gridm
mlabel on
plabel on
framem on

Чтобы найти фактическую зону, можно определить местоположение центра DEM в координатах UTM,

[M,N] = size(Z);
xCenterIntrinsic = (1 + N)/2;
yCenterIntrinsic = (1 + M)/2;
[xCenter, yCenter] = intrinsicToWorld(R,xCenterIntrinsic,yCenterIntrinsic)

xCenter =

              315555.84375


yCenter =

                 4908971.5

затем преобразуйте долготу широты, использовав в своих интересах факт что xCenter и yCenter будет то же самое в зоне 19 Н, как они в фактическую зону.

[latCenter, lonCenter] = minvtran(xCenter,yCenter)

latCenter =

          44.3125403747455


lonCenter =

         -71.3126367639007

Затем с utmzone функция, можно использовать координаты долготы широты, чтобы определить фактическую зону UTM для DEM.

actualZone = utmzone(latCenter,lonCenter)

actualZone =

    '19T'

Наконец, используйте результат сбросить зону осей, созданных ранее.

setm(gca,'zone',actualZone)

Примечание: если можно визуально поместить приблизительно район Нью-Гэмпшира на мировой карте, затем можно подтвердить этот результат с utmzoneui ГРАФИЧЕСКИЙ ИНТЕРФЕЙС ПОЛЬЗОВАТЕЛЯ.

  utmzoneui(actualZone)

Шаг 4: отобразите исходный DEM на осях карты

Используйте mapshow (а не geoshow или meshm) отобразить DEM на осях карты, потому что данные с координатной сеткой в карте (x-y) координаты.

mapshow(Z,R,'DisplayType','texturemap')
demcmap(Z)

DEM покрывает такую небольшую часть этой карты, которую может быть трудно видеть (взгляд между 44 и 44 градусами на север и 72 и 71 градусом на запад), потому что пределы карты устанавливаются, чтобы покрыть целую зону UTM. Можно сбросить их (а также сетка карты и параметры метки), чтобы получить более внимательное рассмотрение.

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)

Когда это просматривается в этом более широком масштабе, узкие области, имеющие форму клина единого цвета появляются вдоль ребра сетки. Это места где Z содержит значение NaN, который указывает на отсутствие фактических данных. По умолчанию они получают первое, раскрашивают таблицу цветов, которая в этом случае темно-зеленая. Эти области пустых данных возникают, потому что несмотря на то, что DEM с координатной сеткой в координатах UTM, его пределы данных заданы четырехугольником долготы широты. Узкий угол каждого клина соответствует ненулевому "наклону сетки" системы координат UTM в этой части зоны. (Линии постоянного x, запущенного точно между севером и югом только вдоль центрального меридиана зоны. В другом месте они следуют за небольшим углом относительно локальных меридианов.)

Шаг 5: задайте Выходную сетку долготы широты

Следующий шаг должен задать расположенный с равными интервалами набор узлов решетки в долготе широты, которая покрывает область в DEM приблизительно при том же пространственном разрешении как исходный набор данных.

Во-первых, необходимо определить, как широта изменяется между строками во входе DEM (т.е. путем перемещения к северу на 30 метров).

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

Фактический интервал может быть округлен немного, чтобы задать интервал сетки, который будет использоваться в выходе (долгота широты) сетка.

gridSpacing = 1/4000;   % In other words, 4000 samples per degree

Чтобы установить степень выхода (долгота широты) сетка, запустите путем нахождения углов DEM в координатах карты UTM.

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

Затем преобразуйте углы в долготу широты. Отобразите углы долготы широты на карте (через проекцию UTM), чтобы проверять, что результаты разумны.

[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

Затем вокруг исходящий, чтобы задать выходной четырехугольник долготы широты, который полностью заключает DEM и выравнивается со множителями интервала сетки.

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

Наконец, создайте географический растровый объект привязки для выходной сетки. Это поддерживает преобразования между долготой широты и индексами строки и столбца. В этом случае использование матрицы файла привязки, W, включает точную спецификацию интервала сетки.

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 полностью задает номер и местоположение каждой ячейки в выходной сетке.

Шаг 6: сопоставьте каждое Выходное местоположение узла решетки с UTM X-Y

Наконец, вы готовы использовать проекцию карты, применяя его к местоположению каждой точки в выходной сетке. Сначала вычислите широты и долготы тех точек, сохраненных в 2D массивах.

[rows,cols] = ndgrid(1:nRows, 1:nCols);
[lat,lon] = intrinsicToGeographic(Rlatlon,cols,rows);

Затем примените проекцию к каждой паре долготы широты, массивам UTM x-y местоположения, имеющие ту же форму и размер как массивы долготы широты.

[XI,YI] = mfwdtran(lat,lon);

На данном этапе XI(i,j) и YI(i,j) задайте координату UTM узла решетки, соответствующего i-th строке и j-th столбцу выходной сетки.

Шаг 7: передискретизируйте исходный DEM

Последний шаг должен использовать interp2 MATLAB функция, чтобы выполнить билинейную передискретизацию.

На данном этапе значение проектирования от сетки долготы широты в систему координат карты UTM становится очевидным: это означает, что передискретизация может произойти в обычной сетке X-Y, делая interp2 применимый. Противоположный подход, не проектируя каждого (X, Y) точка в долготу широты, может казаться более интуитивным, но это привело бы к неправильному массиву точек, которые будут интерполированы - намного более трудная задача, требуя использования намного более дорогостоящего griddata функционируйте или некоторый грубый эквивалент.

[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);

Просмотрите результат в спроектированных осях с помощью geoshow, который повторно спроектирует его на лету. Заметьте, что это заполняет синий прямоугольник, который выравнивается с линиями широты и долготы. (В отличие от этого красный прямоугольник, который обрисовывает в общих чертах исходный DEM, выравнивается с X и Y UTM.) Также замечают NaN-заполненные области вдоль ребер сетки. Контуры этих областей кажутся немного зубчатыми, на уровне одного интервала сетки, из-за эффектов округления во время интерполяции. Переместите красный четырехсторонний и синий четырехугольник в верхнюю часть, чтобы гарантировать, что они не скрыты растровым дисплеем.

geoshow(Zlatlon,Rlatlon,'DisplayType','texturemap')
uistack([hCorners hquad],'top')

format(currentFormat)

Кредиты

9129CATD.ddf (и вспомогательные файлы):

  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