未投影数字高程模型(DEM)gydF4y2Ba

打开脚本gydF4y2Ba

这个示例演示了如何将USGS DEM转换为具有可比空间分辨率的常规经纬度网格。美国地质调查局(USGS) 30米数字高程模型(DEMs)是使用UTM坐标系统的规则网格(光栅数据)。在应用程序中使用这些DEMs可能需要重新投影和重新采样。您可以很容易地将本文介绍的方法应用于UTM之外的投影映射坐标系统以及其他DEMs和大多数类型的常规数据网格。gydF4y2Ba

步骤1:导入DEM及其元数据gydF4y2Ba

这个例子使用了美国地质勘测局的DEM，它是位于美国新罕布什尔州白山中的一个7.5弧分的四边形广场。数据集以空间数据传输标准(STDS)格式存储，并位于文件夹中gydF4y2Ba

fullfile (matlabrootgydF4y2Ba“工具箱”gydF4y2Ba,gydF4y2Ba“地图”gydF4y2Ba,gydF4y2Ba“mapdata”gydF4y2Ba,gydF4y2Ba“项”gydF4y2Ba);gydF4y2Ba

如果安装了Mapping Toolbox™，则此文件夹位于MATLAB®路径上，因此仅通过文件名引用数据集就足够了。gydF4y2Ba

stdsfilename =gydF4y2Ba9129 catd.ddfgydF4y2Ba;gydF4y2Ba

你可以使用gydF4y2BasdtsinfogydF4y2Ba命令获取有关DEM的基本元数据。gydF4y2Ba

信息= sdtsinfo (stdsfilename)gydF4y2Ba
info = struct with fields: Filename: '9129CATD。地区指定基金的标题:“华盛顿山,NH - 24000“ProfileID:“sdt光栅概要文件“ProfileVersion:“1997年7月草案”MapDate:“DataCreationDate:“19980811”HorizontalDatum:“北美1927”MapRefSystem:“UTM”ZoneNumber: 19 XResolution: 30 YResolution: 30 NumberOfRows: 472 NumberOfCols: 345 HorizontalUnits:“米”VerticalUnits:“米”MinElevation: 367 MaxElevation: 1909gydF4y2Ba

你可以使用gydF4y2BasdtsdemreadgydF4y2Ba将DEM导入二维MATLAB阵列(gydF4y2BaZgydF4y2Ba)及其引用矩阵(gydF4y2BarefmatgydF4y2Ba的行和列下标的3×2矩阵gydF4y2BaZgydF4y2Ba映射x和y(本例中为UTM“eastings”和“northing”)。gydF4y2Ba

[Z, refmat] = sdtsdemread (stdsfilename);gydF4y2Ba

您可以转换gydF4y2BarefmatgydF4y2Ba映射光栅参考对象，它提供了一种更完整和自文档化的编码空间参考信息的方法。gydF4y2Ba

currentFormat = get (0,gydF4y2Ba“格式”gydF4y2Ba);格式gydF4y2Ba长gydF4y2BaggydF4y2BaR = refmatToMapRasterReference(refmat, size(Z))gydF4y2Ba
R = MapCellsReference属性:XWorldLimits: [310380.84375 - 320730.84375] YWorldLimits: [4901891.5 - 4916051.5] RasterSize: 345年[472]RasterInterpretation:“细胞”ColumnsStartFrom:“北”RowsStartFrom:“西方”CellExtentInWorldX: 30 CellExtentInWorldY: 30 RasterExtentInWorldX: 10350 RasterExtentInWorldY: 14160 XIntrinsicLimits: [0.5 - 345.5] YIntrinsicLimits: [0.5 - 472.5] TransformationType:“直线”CoordinateSystemType:“平面”gydF4y2Ba

步骤2:指定一个参考椭球gydF4y2Ba

的价值gydF4y2Ba

info.HorizontalDatumgydF4y2Ba
北美1927年gydF4y2Ba

表示使用1927年的北美基准。克拉克1866椭球面是这个基准的标准参考椭球面。gydF4y2Ba

ellipsoidName =gydF4y2Ba“clarke66”gydF4y2Ba;gydF4y2Ba

也可以检查的值gydF4y2BaHorizontalUnitsgydF4y2Ba领域,gydF4y2Ba

mapUnits = info.HorizontalUnits;gydF4y2Ba

这说明DEM的水平坐标以米为单位，并利用这两条信息来构造agydF4y2BareferenceEllipsoidgydF4y2Ba。gydF4y2Ba

clarke66 = referenceEllipsoid(椭球名称，mapUnits)gydF4y2Ba
clarke66 =参考椭球定义性质:代码:7008名称:'Clarke 1866'长度单位:'米'半轴:6378206.4半轴:6356583.8反向展平:294.978698213898离心率:0.0822718542230038附加性质:展平第三展平平均半径表面积体积gydF4y2Ba

步骤3:确定要使用哪个UTM区域并构造映射轴gydF4y2Ba

从gydF4y2BaMapRefSystemgydF4y2BaSDTS信息结构中的字段，gydF4y2Ba

info.MapRefSystemgydF4y2Ba
ANS = 'UTM'gydF4y2Ba

可以看出，DEM是在通用横向墨卡托(UTM)坐标系下进行网格划分的。gydF4y2Ba

的gydF4y2BaZoneNumbergydF4y2Ba场gydF4y2Ba

info.ZoneNumbergydF4y2Ba
ans = 19gydF4y2Ba

指示使用哪个纵向UTM区域。映射的工具箱gydF4y2BautmgydF4y2Ba然而，功能也需要一个纬度带;这在元数据中没有提供，但是您可以从引用的矩阵和网格维度中获得它。gydF4y2Ba

UTM包括60个经度为6度的纵向区域和20个纬向区域，从南纬80度到北纬84度。纵向区域由1到60之间的数字指定。纬向区是由C到X(略去I和O)的字母来表示的。在一个给定的半球(南方或北方)，所有的纬向区都有一个共同的坐标系统。除了确定半球之外，工具箱仅使用纬度带来帮助设置默认地图限制。gydF4y2Ba

所以，你可以从北半球的第一个纬向区开始，N区(北纬0度到8度之间)作为临时区。gydF4y2Ba

longitudinalZone = sprintf (gydF4y2Ba' % d 'gydF4y2Ba,info.ZoneNumber);provisionalLatitudinalZone =gydF4y2Ba“N”gydF4y2Ba;临时带[纵向带临时带]gydF4y2Ba
provisionalZone = ' 19 n 'gydF4y2Ba

然后，根据这个临时区域构造一个UTM轴gydF4y2Ba

图axesm (gydF4y2Ba'mapprojection'gydF4y2Ba,gydF4y2Ba“设备”gydF4y2Ba,gydF4y2Ba“区域”gydF4y2BaprovisionalZone,gydF4y2Ba“大地水准面”gydF4y2Baclarke66)轴gydF4y2Ba从gydF4y2Bagridm mlabelgydF4y2Ba在gydF4y2BaplabelgydF4y2Ba在gydF4y2BaframemgydF4y2Ba在gydF4y2Ba

要找到实际的区域，你可以在UTM坐标下找到DEM的中心，gydF4y2Ba

[M, N] =大小(Z);xCenterIntrinsic = (1 + N)/2;yCenterIntrinsic = (1 + M)/2;[xCenter, yCenter] = intrinsicToWorld(R,xCenter intrinsic,yCenter intrinsic)gydF4y2Ba
xCenter = 315555.84375 yCenter = 4908971.5gydF4y2Ba

然后转换经纬度，利用这个事实gydF4y2BaxCentergydF4y2Ba和gydF4y2BayCentergydF4y2Ba在区域19N和实际区域中是一样的。gydF4y2Ba

[latCenter, lonCenter] = minvtran(xCenter,yCenter)gydF4y2Ba
latCenter = 44.3125403747455 lonCenter = -71.3126367639007gydF4y2Ba

然后,gydF4y2BautmzonegydF4y2Ba函数，可以使用经纬度坐标确定DEM的实际UTM区域。gydF4y2Ba

actualZone = utmzone (latCenter lonCenter)gydF4y2Ba
actualZone = ' 19 t 'gydF4y2Ba

最后，使用结果重新设置先前构造的轴的区域。gydF4y2Ba

setm (gca),gydF4y2Ba“区域”gydF4y2BaactualZone)gydF4y2Ba

注意:如果你能在世界地图上直观地显示出新罕布什尔州的大致位置，那么你就可以用gydF4y2BautmzoneuigydF4y2BaGUI。gydF4y2Ba

utmzoneui（actualZone）gydF4y2Ba

步骤4:在地图轴线上显示原始DEMgydF4y2Ba

使用gydF4y2BamapshowgydF4y2Ba(而不是gydF4y2BageoshowgydF4y2Ba或gydF4y2BameshmgydF4y2Ba)显示在地图轴上的DEM，因为数据是在地图(x-y)坐标中网格化的。gydF4y2Ba

mapshow (Z, R,gydF4y2Ba“DisplayType”gydF4y2Ba,gydF4y2Ba“texturemap”gydF4y2Ba)demcmap (Z)gydF4y2Ba

DEM只覆盖了地图的一小部分，所以可能很难看到(北纬44度到44度之间，西经72度到71度之间)，因为地图的界限被设置为覆盖整个UTM区域。您可以重置它们(以及地图网格和标签参数)以获得更近距离的查看。gydF4y2Ba

setm (gca),gydF4y2Ba“MapLatLimit”gydF4y2Ba(44.2 - 44.4),gydF4y2Ba“MapLonLimit”gydF4y2Ba甘氨胆酸,[-71.4 - -71.2])setm (,gydF4y2Ba“MLabelLocation”gydF4y2Ba,0.05,gydF4y2Ba“MLabelRound”gydF4y2Ba甘氨胆酸,2)setm (,gydF4y2Ba“PLabelLocation”gydF4y2Ba,0.05,gydF4y2Ba“PLabelRound”gydF4y2Ba甘氨胆酸,2)setm (,gydF4y2Ba“PLineLocation”gydF4y2Ba,0.025,gydF4y2Ba“MLineLocation”gydF4y2Ba,0.025)gydF4y2Ba

当它在这个更大的尺度上被观察的时候，狭窄的楔形区域的均匀颜色沿着网格的边缘出现。这些地方gydF4y2BaZgydF4y2Ba包含值NaN，表示实际数据的缺失。默认情况下，它们接收颜色表中的第一种颜色，在本例中是深绿色。这些空数据区产生的原因是，尽管DEM是在UTM坐标中进行网格划分的，但它的数据限制是由经纬度四边形定义的。每个楔块的窄角对应于该区域UTM坐标系的非零“格点赤纬”。常数x线只沿着该区域的中央子午线沿南北向走。在其他地方，它们遵循一个相对于当地经脉的小角度。)gydF4y2Ba

步骤5:定义输出经纬度网格gydF4y2Ba

下一步是在经纬度上定义一组规则间隔的网格点，这些点覆盖DEM内的区域，其空间分辨率与原始数据集大致相同。gydF4y2Ba

首先，您需要确定输入DEM(即，向北移动30米)。gydF4y2Ba

rng = info.YResolution;gydF4y2Ba以米为单位，与克拉克66一致gydF4y2Balatcrad =函数(latCenter);gydF4y2Ba以弧度表示的% latCentergydF4y2Ba纬度变化百分比，角度变化百分比gydF4y2BadLat = rad2deg(meridianfwd(latcrad,rng,clarke66) - latcrad)gydF4y2Ba
dLat = 0.000269984934404749gydF4y2Ba

可以稍微四舍五入实际的间距，以定义用于输出(经纬度)网格的网格间距。gydF4y2Ba

gridSpacing = 1/4000;gydF4y2Ba换句话说，每度4000个样品gydF4y2Ba

要设置输出(经纬度)网格的范围，首先在UTM地图坐标中查找DEM的角。gydF4y2Ba

xCorners = R。XWorldLimits([1 1 2 2])' yCorners = R。YWorldLimits([1,2,2,1])'gydF4y2Ba
xCorners = 310380.84375 310380.84375 320730.84375 320730.84375 yCorners = 4901891.5 4916051.5 4916051.5 4901891.5gydF4y2Ba

然后将拐角转换为经纬度。在地图上显示经纬度角(通过UTM投影)，以检查结果是否合理。gydF4y2Ba

[latCorners, lonCorners] = minvtran(xCorners, yCorners) hCorners = geoshow(latCorners,lonCorners，gydF4y2Ba“DisplayType”gydF4y2Ba,gydF4y2Ba“多边形”gydF4y2Ba,gydF4y2Ba…gydF4y2Ba“FaceColor”gydF4y2Ba,gydF4y2Ba“没有”gydF4y2Ba,gydF4y2Ba“EdgeColor”gydF4y2Ba,gydF4y2Ba“红色”gydF4y2Ba);gydF4y2Ba
latCorners = 44.24752129387 44.3775277222457 44.2501421324565 lonCorners = -71.374900167689 -71.3800449718219 -71.2502372050341 - 71.24537066789gydF4y2Ba

接下来，向外四舍五入定义一个输出经纬度四边形，它完全包围DEM并与网格间距的倍数对齐。gydF4y2Ba

latSouth =网格间距*地板(最小(网格角)/网格间距)lonWest =网格间距*地板(最小(网格角)/网格间距)latNorth =网格间距* ceil(最大(网格角)/网格间距)lonEast =网格间距* ceil(最大(网格角)/网格间距)qlatlim = [latSouth latNorth];qlonlim = [lonWest lonEast];dlat = 100 * gridSpacing;dlon = 100 * gridSpacing;[latquad, lonquad] = outlinegeoquad(qlatlim, qlonlim, dlat, dlon);geoshow(latquad, lonquad，gydF4y2Ba…gydF4y2Ba“DisplayType”gydF4y2Ba,gydF4y2Ba“多边形”gydF4y2Ba,gydF4y2Ba“FaceColor”gydF4y2Ba,gydF4y2Ba“没有”gydF4y2Ba,gydF4y2Ba“EdgeColor”gydF4y2Ba,gydF4y2Ba“蓝”gydF4y2Ba);snapnow;gydF4y2Ba
latSouth = 44.2475 lonWest = -71.38025 latNorth = 44.37775 lonEast = -71.24525gydF4y2Ba

最后，为输出网格构造一个引用对象的地理栅格。它支持金宝app经纬度和行、列下标之间的转换。在这种情况下，使用世界文件矩阵W可以精确指定网格间距。gydF4y2Ba

W =[格网间距0 lonWest +格网间距/2;gydF4y2Ba…gydF4y2Ba0网格间距latSouth +网格间距/2]gydF4y2Ba
W = 0.00025 0 -71.380125 0 0.00025 44.247625gydF4y2Ba
NROWS = ROUND（（latNorth  -  latSouth）/ gridSpacing）NCOLS = ROUND（wrapTo360（lonEast  -  lonWest）/ gridSpacing）gydF4y2Ba
nRows = 521 nCols = 540gydF4y2Ba
Rlatlon = georasterref(W，[nRows nCols])，gydF4y2Ba“细胞”gydF4y2Ba)gydF4y2Ba
Rlatlon = GeographicCellsReference与属性：LatitudeLimits：[44.2475 44.37775] LongitudeLimits：-71.38025 -71.24525] RasterSize：[521 540] RasterInterpretation： '细胞' ColumnsStartFrom： '南' RowsStartFrom： '西方' CellExtentInLatitude：1/4000 CellExtentInLongitude：1 /4000 RasterExtentInLatitude：0.13025 RasterExtentInLongitude：0.135 XIntrinsicLimits：[0.5 540.5] YIntrinsicLimits：[0.5 521.5] CoordinateSystemType： '地理' AngleUnit： '度'gydF4y2Ba

RlatlongydF4y2Ba完全定义输出网格中每个单元格的数量和位置。gydF4y2Ba

步骤6:将每个输出网格点位置映射到UTM X-YgydF4y2Ba

最后，你就可以利用地图投影，它应用于输出网格的每个点的位置。第一计算这些点，存储在2-d阵列的纬度和经度。gydF4y2Ba

[rows,cols] = ndgrid(1:nRows, 1:nCols);(纬度、经度)= intrinsicToGeographic (Rlatlon关口,行);gydF4y2Ba

然后将投影应用到每个经纬度对，即具有与经纬度数组相同形状和大小的UTM x-y位置数组。gydF4y2Ba

[XI, YI] = mfwdtran(纬度、经度);gydF4y2Ba

在这一点上,gydF4y2Ba习(i, j)gydF4y2Ba和gydF4y2Ba易(i, j)gydF4y2Ba指定输出网格的第i行和第j列对应的网格点的UTM坐标。gydF4y2Ba

第七步:重新采样原始DEMgydF4y2Ba

最后一步是使用MATLABgydF4y2Bainterp2gydF4y2Ba函数执行双线性重采样。gydF4y2Ba

在这个阶段，将经纬度网格投影到UTM地图坐标系统的价值变得很明显:这意味着重新采样可以在常规的X-Y网格中进行gydF4y2Bainterp2gydF4y2Ba适用。相反的方法，将每个(X,Y)点不投影到经纬度，可能看起来更直观，但它会导致一个不规则的点数组被插值——这是一个更困难的任务，需要使用成本高得多的方法gydF4y2BagriddatagydF4y2Ba函数或一些粗略的等价。gydF4y2Ba

(行,关口)= ndgrid (1: M, 1: N);(X, Y) = intrinsicToWorld (R,关口,行);方法=gydF4y2Ba双线性的gydF4y2Ba;extrapval =南;Zlatlon = interp2 (X, Y, Z, XI,咦,方法,extrapval);gydF4y2Ba

使用投影轴查看结果gydF4y2BageoshowgydF4y2Ba，这将再次投射它的飞行。注意，它填补了蓝色矩形，其与纬度和经度的线对准。（与此相反，红色矩形，其概述了原始DEM，用UTM对准X和Y）。还要注意沿网格的边缘的NaN填充区域。这些区域的边界出现轻微锯齿状，在一个单一的栅格间距的水平，由于插补时的舍入效果。移动红色四边形，蓝色四边形顶端，以确保它们不被光栅显示隐藏。gydF4y2Ba

geoshow (Zlatlon Rlatlon,gydF4y2Ba“DisplayType”gydF4y2Ba,gydF4y2Ba“texturemap”gydF4y2Ba)uistack ([hCorners hquad),gydF4y2Ba“高级”gydF4y2Ba)gydF4y2Ba

格式(currentFormat)gydF4y2Ba

学分gydF4y2Ba

9129 catd。ddf(及支持文件金宝app):gydF4y2Ba

美国地质调查局(USGS) 7.5分钟数字高程模型(DEM)的空间数据传输标准(SDTS)格式的华盛顿山四边形，以米为单位。http://edc.usgs.gov/下载188bet金宝搏products/elevation/dem.htmlgydF4y2Ba
有关更多信息，请运行:gydF4y2Ba
> > 9129.三种类型gydF4y2Ba

另请参阅gydF4y2Ba

|gydF4y2Ba|gydF4y2Ba|gydF4y2Ba|gydF4y2Ba

映射工具的文档gydF4y2Ba
金宝app

尝试MATLAB、Sim金宝appulink和其他产品下载188bet金宝搏gydF4y2Ba

得到审判现在gydF4y2Ba