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

В этом примере USGS DEM используется для квадратного квадратика длиной 7,5 дуги, расположенного в Белых горах Нью-Гэмпшира, США. Импорт данных и ссылки на ячейки карты с помощью readgeoraster функция. Получение дополнительных метаданных с помощью georasterinfo функция.

[Z,R] = readgeoraster('MtWashington-ft.grd','OutputType','double');
info = georasterinfo('MtWashington-ft.grd');

Заменить отсутствующие данные на NaN значения.

m = info.MissingDataIndicator;
Z = standardizeMissing(Z,m);

Шаг 2: Получение информации о проекции

Получение информации о спроецированной системе координат путем запроса ProjectedCRS свойства ссылочного объекта. Результатом является projcrs объект. Затем получите эллипсоид для системы координат.

p = R.ProjectedCRS;
ellipsoid = p.GeographicCRS.Spheroid

ellipsoid = 

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

От Name имущества projcrs можно указать, что DEM имеет сетку в системе координат универсального поперечного меркатора (UTM).

p.Name

ans = 

    "UTM Zone 19, Northern Hemisphere"

Чтобы найти зону UTM, сначала найдите центр DEM в координатах UTM. Затем преобразуйте координаты в широту-долготу.

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

latCenter =

          44.3124367104673


lonCenter =

         -71.3126432693478

Найдите зону UTM для DEM с помощью utmzone функция.

utmZone = utmzone(latCenter,lonCenter)

utmZone =

    '19T'

Используйте зону и эллипсоид для создания осей карты.

figure
axesm('utm','zone',utmZone,'geoid',ellipsoid)
axis off
gridm
mlabel on
plabel on
framem on

Примечание: если вы можете визуально разместить приблизительное местоположение Нью-Гэмпшира на карте мира, то вы можете подтвердить этот результат с помощью 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 = R.CellExtentInWorldY;  % In meters, consistent with p.LengthUnit
latcrad = deg2rad(latCenter);   % latCenter in radians

% Change in latitude, in degrees
dLat = rad2deg(meridianfwd(latcrad,rng,ellipsoid) - latcrad)

dLat =

      0.000269984939366415

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

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
      310380
      320730
      320730


yCorners =

     4901880
     4916040
     4916040
     4901880

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

[latCorners, lonCorners] = projinv(p,xCorners, yCorners)
hCorners = geoshow(latCorners,lonCorners,'DisplayType','polygon',...
    'FaceColor','none','EdgeColor','red');

latCorners =

          44.2474175605687
          44.3747915486804
          44.3774240601986
          44.2500384686392


lonCorners =

         -71.3749065609587
         -71.3800513603087
         -71.2502438233865
         -71.2453790282992

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


lonWest =

                 -71.38025


latNorth =

                   44.3775


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.247375

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

nRows =

   521


nCols =

   540

Rlatlon = georasterref(W,[nRows nCols],'cells');
Rlatlon.GeographicCRS = p.GeographicCRS

Rlatlon = 

  GeographicCellsReference with properties:

             LatitudeLimits: [44.24725 44.3775]
            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'
              GeographicCRS: [1x1 geocrs]
                  AngleUnit: 'degree'

Rlatlon полностью определяет количество и местоположение каждой ячейки в выходной сетке.

Шаг 6: Сопоставление местоположения каждой выходной точки сетки с UTM X-Y

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

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

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

[XI,YI] = projfwd(p,lat,lon);

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

Шаг 7: Повторная выборка исходного DEM

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

На этом этапе становится очевидным значение проецирования из сетки широты-долготы в систему координат карты 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, выравнивается с UTM x и y.) Также обратите внимание на области, заполненные NaN, вдоль кромок сетки. Границы этих областей выглядят слегка зазубренными, на уровне одного интервала сетки, из-за эффектов скругления во время интерполяции. Переместите красный четырехугольник и синий четырехугольник наверх, чтобы убедиться, что они не скрыты растровым дисплеем.

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

Восстановите исходный формат отображения выходных данных.

format(currentFormat)

Кредиты

MtWashington-ft.grd (и вспомогательные файлы):

  United States Geological Survey (USGS) 7.5-minute Digital Elevation
  Model (DEM) for the Mt. Washington quadrangle, with elevation in
  meters. http://edc.usgs.gov/products/elevation/dem.html

  For more information, run:

  >> type MtWashington-ft.txt