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

Question

我有一张带有 lat/long 坐标中的高程栅格的地图，但我现在需要在 easting/northing 中显示地图。我怎样才能做到这一点？

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

一个最小的工作示例遵循我用来在 lat/long 坐标中生成地图的代码的描述。我使用手动定义的边界框 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)

据此，如何在此地图上包含适当的轴（在 easting/northing 中）？

Answer 1

I now need to display the map in easting/northing.

我了解到您想使用平面（笛卡尔）坐标而不是 angular (lon/lat) 坐标。这意味着您需要选择合适的坐标参考系统 (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。那会有点奇怪。通常，人们可能希望在其他投影（平面）地图上显示 longlat 坐标。

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

How can I project a raster in lat-long coordinates to UTM, for plotting in tmap?

maps

r

map-projections

tmap