Шаг 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 с координатной сеткой в системе координат Universal, поперечной меркаторской (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.3124


lonCenter =

  -71.3126

Найдите зону 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 (i.e., путем перемещения к северу на 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 =

   2.6998e-04

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

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.2474
   44.3748
   44.3774
   44.2500


lonCorners =

  -71.3749
  -71.3801
  -71.2502
  -71.2454

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


lonWest =

  -71.3803


latNorth =

   44.3775


lonEast =

  -71.2452

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

currentFormat = format; % since R2021a
format longG
W = [gridSpacing    0              lonWest + gridSpacing/2; ...
     0              gridSpacing    latSouth + gridSpacing/2]
format(currentFormat)

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

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

[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

Последний шаг должен использовать 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')

Кредиты

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