Шаг 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);

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

[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