Этот пример показывает, как преобразовать USGS DEM в обычную сетку долготы широты, имеющую сопоставимое пространственное разрешение. Американская Геологическая служба (USGS), 30-метровые Цифровые Модели Повышения (демократы) являются обычными сетками (растровые данные), которые используют систему координат UTM. Используя таких демократов в приложениях может потребовать перепроектирования и передискретизации их. Можно с готовностью применяться, подход показывают здесь спроектированным системам координат карты кроме UTM и другим демократам и большинству типов обычных сеток данных.
Этот пример использует 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'
Значение
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
От поля MapRefSystem в информационном struct 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)
Используйте 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, запущенного точно между севером и югом только вдоль центрального меридиана зоны. В другом месте они следуют за небольшим углом относительно локальных меридианов.)
Следующий шаг должен задать расположенный с равными интервалами набор узлов решетки в долготе широты, которая покрывает область в 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 полностью задает номер и местоположение каждой ячейки в выходной сетке.
Наконец, вы готовы использовать проекцию карты, применяя его к местоположению каждой точки в выходной сетке. Сначала вычислите широты и долготы тех точек, сохраненных в 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 столбцу выходной сетки.
Последний шаг должен использовать функцию 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
demcmap | georasterref | intrinsicToGeographic | intrinsicToWorld | refmatToMapRasterReference