解除数字高程模型(DEM)的投影
本示例展示了如何将USGS DEM转换为具有可比空间分辨率的常规经纬度网格。美国地质调查局(USGS) 30米数字高程模型(dem)是使用UTM坐标系统的常规网格(栅格数据)。在应用程序中使用这样的dem可能需要重新投影和重新采样。您可以很容易地将这里展示的方法应用于投影地图坐标系，而不是UTM、其他dem和大多数常规数据网格类型。
首先，设置输出显示格式为longG因此，输出显示更多小数位。获取当前输出格式，以便在示例的末尾恢复它。
CurrentFormat = Get（0，'格式'）;格式longG
步骤1:导入DEM及其元数据
本例使用USGS DEM对位于美国新罕布什尔州白山的一个7.5弧分四边形正方形进行分析。控件导入数据和映射单元格引用对象readgeoraster功能。使用额外的元数据使用georasterinfo功能。
(Z, R) = readgeoraster (“MtWashington-ft.grd”那'OutputType'那“双”）;信息= georasterinfo (“MtWashington-ft.grd”）;
用替换丢失的数据南价值观。
m = info.missingDataindicator;z =标准化（z，m）;
第二步:获取投影信息
通过查询来获取有关投影坐标参考系统的信息ProjectedCRS属性的。结果是projcrs.对象。然后，得到坐标系的椭球面。
p = r.projectedcrs;ellipsoid = p.geographiccrs.phoid.
椭球= referenceEllipsoid with defining properties: Code: 7008 Name: 'Clarke 1866' LengthUnit: 'meter' SemimajorAxis: 6378206.4 SemiminorAxis: 6356583.8 inverse扁平化:294.978698213898偏心:0.0822718542230038和附加属性:扁平化第三扁平化平均半径表面积体积
步骤3:确定哪个UTM区域使用和构建一个映射轴
来自名称财产的财产projcrs.对象，您可以告诉DEM在通用横向符号（UTM）坐标系中网格网。
p.name.
ans = "北半球，UTM 19区"
要找到UTM区域，首先在UTM坐标中找到DEM的中心。然后，将坐标转换为纬度-经度。
[m，n] =尺寸（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找到utm区域utmzone功能。
Utmzone = Utmzone（LatCenter，Loncenter）
Utmzone ='19t'
使用区域和椭圆体创建地图轴。
图轴（'UTM'那“区域”，Utmzone，'Geoid'椭球)轴从gridm mlabel.上plabel.上弗拉姆姆上
注意：如果您可以在线地将新罕布什尔州的大致位置放在世界地图上，那么您可以使用此结果确认此结果Utmzoneui.吉..
Utmzoneui（实际Zone）
第4步：在地图轴上显示原始DEM
使用mapshow（而不是地球机或者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 (,'mlabellocation'，0.05，'mlabelround'甘氨胆酸,2)setm (,“PLabelLocation”，0.05，“PLabelRound”甘氨胆酸,2)setm (,“PLineLocation”，0.025，'mlinelocation'，0.025）
当以这种较大的尺度观察时，沿着栅格的边缘出现窄的楔形区域。这些是哪里Z.包含值NaN，表示缺少实际数据。默认情况下，它们接收颜色表中的第一个颜色，在本例中是深绿色。出现这些空数据区域的原因是，尽管DEM是在UTM坐标中网格化的，但它的数据限制是由经纬度四边形定义的。每个楔形的窄角对应于该区域内UTM坐标系的非零“网格赤纬”。(常数x线只沿着该区域的中央子午线精确地朝南运行。在其他地方，它们与当地的经线略有夹角。)
第5步：定义输出纬度 - 经度网格
下一步是在纬度-经度上定义一组规则间隔的网格点，以与原始数据集相同的空间分辨率覆盖DEM内的区域。
首先，您需要确定输入DEM中行与行之间的纬度变化(即向北移动30米)。
rng = r.cellextentinworldy;%以米为单位，与p一致latcrad =函数(latCenter);弧度中的％latcenter纬度，度数的变化DLAT = RAD2DEG（Meridianfwd（Latcrad，RNG，Ellipsoid） - 拉拉德）
DLAT = 0.000269984939366415
实际间距可以略微舍入以定义用于输出（纬度 - 经度）网格的网格间隔。
gridSpacing = 1/4000;换句话说，每度4000个样品
要设置输出的范围（纬度 - 经度）网格，请首先查找UTM地图坐标中的DEM的角。
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”那'多边形'那......“FaceColor”那'没有任何'那'Edgecolor'那'红色的'）;
latCorners = 44.2474175605687 44.3747915486804 44.3774240601986 44.2500384686392 lonCorners = -71.3749065609587 -71.3800513603087 -71.2502438233865 -71.2453790282992
接下来，向外舍入以定义一个完全包围DEM的输出纬度 - 经度四边形，并用网格间距的倍数对准。
Latsouth = Gridspacing *地板（Min（Latcorners）/ Gridspacing）Lonwest = Gridspacing *地板（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”那'多边形'那“FaceColor”那'没有任何'那'Edgecolor'那“蓝”）;snapnow;
Latsouth = 44.24725 Lonwest = -71.38025 Latnorth = 44.3775 Loneast = -71.24525
最后，为输出网格构造一个地理栅格参考对象。它支持金宝app纬度-经度以及行和列下标之间的转换。在这种情况下，使用世界文件矩阵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（lonease - lonwest）/ gridspacing）
nrows = 521 ncols = 540
rlatlon = georasterref（w，[nrows ncols]，'细胞'）;Rlatlon。GeographicCRS = p.GeographicCRS
Rlatlon = GeopraphicCellsReference有属性：LATITudelimits：[44.24725 44.3775] LogitudeLimits：[-71.38025 -71.24525] rastersize：[521 540] rastersize：rastersize：'cells'columnsstartfrom：'南'Rowsstartfrom：'West'CelleXtentinlitude：1/4000 CelleXtentinlongitude：1 /4000 RASTEREXTENTINLATICE：0.13025 RASTEREXTENTINLONGE：0.135 XintRinsiclimits：[0.5 540.5]金林林利：[0.5 521.5] COOLINATESYSTEMTYPE：'地理'地理（地理）：[1x1 GeOcrs] AnallUnit：'学位'
Rlatlon完全定义输出网格中每个单元格的数字和位置。
步骤6：将每个输出网格点位置映射到UTM x-y
最后，您已准备好使用地图投影，将其应用于输出网格中的每个点的位置。首先计算存储在2-D阵列中的那些点的纬度和长度。
[行，cols] = ndgrid（1：nrows，1：ncols）;[LAT，LON] =内在地理（Rlatlon，Cols，Rows）;
然后将投影应用于每个纬度对对，具有与纬度度量阵列相同的形状和大小的UTM X-Y位置阵列。
[xi，yi] = projfwd（p，lat，lon）;
在此刻，XI（I，J）和易(i, j)指定与输出网格的第i行和第j列对应的网格点的UTM坐标。
第7步：重新取决于原始DEM
最后一步是使用MATLABinterp2功能执行双线性重采样。
在此阶段，从纬度经度网格投射到UTM地图坐标系的值变得明显：这意味着重采样可以在常规X-Y网格中进行，制作interp2适用的。反向方法，将每个（x，y）指向纬度经度，可能似乎更直观，但它会导致要插值的不规则点数 - 更难的任务，要求使用更昂贵的昂贵griddata功能或一些粗糙的等价物。
[行，cols] = ndgrid（1：m，1：n）;[x，y] = intrinsictoworld（r，cols，行）;方法='双线性';extrapval =南;Zlatlon = interp2 (X, Y, Z, XI,咦,方法,extrapval);
使用投影轴查看结果地球机，它将在飞行中重新投入。请注意，它填充了蓝色矩形，这与纬度和经度线对齐。（相比之下，概述原始DEM的红色矩形，与UTM X和Y对齐。）还注意到沿网格边缘的纳米填充的区域。由于在插值期间的圆截止效应，这些区域的边界在单个网格间距的水平上显得略微锯齿。将红色四边形和蓝色四边形移动到顶部，以确保它们不会被光栅显示隐藏。
Geoshow（Zlatlon，Rlatlon，“DisplayType”那“texturemap”）Uistack（[Hcorners HQuad]，'最佳'）
恢复原始输出显示格式。
格式(currentFormat)
学分
MtWashington-ft。GRD(及支持文件金宝app):
美国地质调查（USGS）7.5分钟数字海拔模型（DEM）为Mt.华盛顿四边形，高涨以米。http://edc.usgs.gov/下载188bet金宝搏products/elevation/dem.html.
有关更多信息，请运行：
>>型mtwashington-ft.txt
也可以看看
demcmap.|Georasterref.|内在地理|IntrinsictOWORLD.|refmatToMapRasterReference
你也可以从以下列表中选择一个网站:
