В этом примере показано, как преобразовать 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
поле в информационной структуре 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