CODEWITHSUNDEEP
rgeospatialterra
robert
robert

Reputation: 21

Projecting a raster in terra fails

I have a single simple raster in EPSG:7532 that I am trying to project to EPSG:4326 but is failing

The source data is a Lidar point clould that I am able to process using the lidR package. The data source is in the link below

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²

convert to a 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

The point of failure is the attempt to project to geographic format:

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

Upvotes: 1

Views: 752

Answers (1)

dimfalk
dimfalk

Reputation: 1413

This issue is not resulting from reprojection from EPSG:7532 to EPSG:4326 per se, but seems rather connected to the fact, that your SpatRaster object created via rasterize_canopy() comes with a vertical datum, apparently causing problems downstream:

VERTCRS["NAVD88 height - Geoid18 (m)",
        VDATUM["North American Vertical Datum 1988"],
        CS[vertical,1],
            AXIS["up",up,
                LENGTHUNIT["meter",1]],
        GEOIDMODEL["GEOID18"],
        ID["EPSG",5703]]]

The quick & dirty solution would be to simply override crs definition by EPSG:7532 and dropping references in Z dimension, although this does not feel 100 % right. On the other hand, I'm not sure how terra handles vertical crs information and if it is possible to keep this information at all.

library(lidR)
library(terra)
#> terra 1.6.49

l1 = readLAS("USGS_LPC_WI_BrownRusk_2020_B20_02531702.laz")
#> Warning: There are 53206 points flagged 'withheld'.

dsm <- rasterize_canopy(l1, res = 1.0, pitfree(c(0,2,5,10,15), c(0, 1.5))) 

crs(dsm) <- "epsg:7532"

dsm_4326 <- project(dsm, "epsg:4326", method="bilinear")
dsm_4326
#> class       : SpatRaster 
#> dimensions  : 1235, 1727, 1  (nrow, ncol, nlyr)
#> resolution  : 1.093727e-05, 1.093727e-05  (x, y)
#> extent      : -88.0786, -88.05972, 44.49092, 44.50443  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source(s)   : memory
#> name        :        Z 
#> min value   : 185.8360 
#> max value   : 295.1357

From my personal point of view, it would be better to have both crs (xy, z) listed in the SpatRaster summary with reproject() being able to address z-transformations separately (or have e.g. terra::transform()), e.g. when you wanted to just transform your heights but still keep EPSG:7532 as your xy reference system.

coord. ref. xy: NAD83(2011) / WISCRS Brown (m) (EPSG:7532) 
coord. ref. z:  NAVD88 height - Geoid18 (m) (EPSG:5703)

Upvotes: 0

Related Questions