我在EPSG：7532中有一个简单栅格，尝试将其投影到EPSG：4326，但失败了
源数据是一个我可以使用lidR包处理的激光雷达点云。数据源在下面的链接中
https://rockyweb.usgs.gov/vdelivery/Datasets/Staged/Elevation/LPC/Projects/WI_BrownRusk_2020_B20/WI_Brown_2_2020/LAZ/USGS_LPC_WI_BrownRusk_2020_B20_02531702.laz
l1 = readLAS("USGS_LPC_WI_BrownRusk_2020_B20_02531702.laz")
> l1
class : LAS (v1.4 format 6)
memory : 2.5 Gb
extent : 25349, 26849, 170258, 171758 (xmin, xmax, ymin, ymax)
coord. ref. : NAD83(2011) / WISCRS Brown (m) + NAVD88 height - Geoid18 (m)
area : 2.25 km²
points : 35.57 million points
density : 15.79 points/m²
density : 12.89 pulses/m²转换为spatRaster：
dsm <- rasterize_canopy(l1, res = 1.0, pitfree(c(0,2,5,10,15), c(0, 1.5)))
> dsm
class : SpatRaster
dimensions : 1500, 1501, 1 (nrow, ncol, nlyr)
resolution : 1, 1 (x, y)
extent : 25349, 26850, 170258, 171758 (xmin, xmax, ymin, ymax)
coord. ref. : NAD83(2011) / WISCRS Brown (m) (EPSG:7532)
source : memory
name : Z
min value : 185.836
max value : 333.709失败之处在于试图投影到地理格式：
dsm_test <- terra::project(dsm, "EPSG:4326", method="bilinear")
> dsm_test <- terra::project(dsm, "EPSG:4326", method="bilinear")
Error: [project] cannot get output boundaries
In addition: Warning messages:
1: In x@ptr$warp(SpatRaster$new(), y, method, mask, FALSE, opt) :
GDAL Error 1: PROJ: vgridshift: could not find required grid(s).
2: In x@ptr$warp(SpatRaster$new(), y, method, mask, FALSE, opt) :
GDAL Error 1: PROJ: pipeline: Pipeline: Bad step definition: proj=vgridshift (File not found or invalid)
3: In x@ptr$warp(SpatRaster$new(), y, method, mask, FALSE, opt) :
GDAL Error 1: Too many points (961 out of 961) failed to transform, unable to compute output bounds.
A similar topic here, but seems different.
https://stackoverflow.com/questions/72404897/what-is-causing-this-raster-reprojection-error
1条答案
按热度按时间qvtsj1bj1#
此问题并非由从EPSG：7532到EPSG：4326的重投影本身引起，而是与以下事实有关：通过
rasterize_canopy()创建的SpatRaster对象带有垂直基准面，这显然会导致下游问题：最快的解决方案是简单地覆盖EPSG：7532中的crs定义，并删除Z维中的引用，尽管这感觉不是100%正确。另一方面，我不确定terra如何处理垂直crs信息，以及是否有可能保留这些信息。
从我个人的Angular 来看，最好在SpatRaster摘要中列出两个crs（xy，z），
reproject()能够单独处理z变换（或具有terra::transform()），例如，当您只想变换高度但仍将EPSG：7532作为xy参考系时。