如何将经纬度坐标中的栅格投影到 UTM，以便在 tmap 中绘图？

Posted 2023-03-23

【中文标题】如何将经纬度坐标中的栅格投影到 UTM，以便在 tmap 中绘图？

【英文标题】：How can I project a raster in lat-long coordinates to UTM, for plotting in tmap?

【发布时间】：2021-12-25 15:47:40

【问题描述】：

我有一张带有高程栅格的地图，它采用纬度/经度坐标，但我现在需要以东/北显示地图。我怎样才能做到这一点？

我尝试将栅格从 lat/long 重新投影到 UTM，但这会扭曲地图（我假设为 reasons discussed in this SO post）。

一个最小的工作示例遵循我用来在纬度/经度坐标中生成地图的代码的描述。我使用手动定义的边界框extent 从FedData 下载感兴趣的高程栅格。然后我使用tmap 绘制了栅格。

# - Libraries ----
library(proj4)
library(FedData)
library(rgdal)
library(tmap)

# extent defined manually
extent <- rgeos::readWKT("POLYGON((-118.25 36.45, -118.25 36.6, -118.25 36.9, -118.8 36.9, -118.8 36.45, -118.25 36.45))")

# define shape polygon on lat/long coordinates
proj4string(extent) <- "+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs"

# - Get national elevation raster ----
# download National Elevation Database elevation data in extent
# default resolution is 1 arcsecond (res="1);
# to get 1/3 arcsecond (res="13)

ned_raster<-FedData::get_ned(template=extent, label="ned_raster", res="1", force.redo = T)

# - Plot with tmap ----
tm_shape(ned_raster,unit.size=1)+      
  tm_graticules() +
  tm_raster(legend.show=FALSE) 

---

更新 在 Robert Hijmans 的回答之后，我正在更新问题，并提供一些额外的说明。我在这里包含了我用于重新投影的代码：

ned_raster2 <- raster::projectRaster(ned_raster, crs=CRS("+proj=utm +init=epsg:26711 +zone=11 +north +no_defs"))

输出与您在答案中包含的地图中的输出相当。通过使用tmaptools::bb 修改边界框，我可以使用tmap 剪辑地图，以便投影正确但不会“出现”扭曲：

# - Project from angular to planar ----
ned_raster2 <- raster::projectRaster(ned_raster, crs=CRS("+proj=utm +init=epsg:26711 +zone=11 +north +no_defs"))

# redefine extent/bounding box
e2 = tmaptools::bb(ned_raster2)
e2 = e2*c(1.01,1.001,.99,.999)

# - Plot with tmap ----
tm_shape(ned_raster2,unit.size=1,bbox=e2)+      
  tm_graticules(n.x=5) +
  tm_raster(legend.show=FALSE) 

据此，我怎样才能在这张地图上包含适当的轴（东/北）？

【参考方案1】：

我现在需要以东/北显示地图。

我了解您希望使用平面（笛卡尔）坐标而不是角度（经度/纬度）坐标。这意味着您需要选择合适的坐标参考系 (CRS)。你说你使用了 UTM，这可能是合理的，但你发现结果是“扭曲的”。您应该显示您使用的代码，因为您可能犯了一个错误。总会有一些失真，但如果您正确指定 CRS，那么对于这样的小区域不会有问题。 （否则，请解释“扭曲”及其重要性。）

例如

library(raster)
library(FedData)

# extent defined manually
e <- as(extent(-118.8, -118.25, 36.45, 36.9), "SpatialPolygons")
crs(e) <- "+proj=longlat"

ned <- FedData::get_ned(template=e, label="ned_raster", res="1", force.redo = T)
ned2 <- projectRaster(ned, crs="+proj=utm +zone=11")
plot(ned2)

对于较大的区域，您不能使用 UTM，您需要使用另一个 CRS。

或者，您可以从投影栅格创建网格线（经纬网），存储它们的坐标，然后将它们转换回 longlat。那会有点奇怪。通常，人们可能希望在其他投影（平面）地图上显示经纬度坐标。

【讨论】：

感谢您对我的问题的清晰解释和重述。我已经用我用来从角度投影到平面的代码更新了我的示例，并尝试在我的问题中添加一些内容。