Search code examples
rrasterterra

terra raster values lost after projection

values of a SpatRaster are lost after changing the coordinate reference system. I do not see any reasons why.

library(terra)

ext <-
  terra::ext(
    9757195,
    9853641,
    734695,
    799794 
  )

r <-
  terra::rast(ext,
              resolution = 2000,
              crs = "EPSG:6933")

I create a SpatVector points geometry to then overlay with my raster and identify in which cells of the raster the points fall. This operation is done in a projected CRS.

coord_vec <- data.frame( x = c(9849641, 9761195), y = c(795794.8, 738695.7))

coord_vec <- terra::vect(coord_vec, 
                         crs =  "EPSG:6933", geom=c("x", "y"))
r2_ <-
  terra::rasterize(x = coord_vec, y = r)

I want to get back to Geodetic coordinate system but then values are lost.

r2_proj <- terra::project(x = r2_,
               y = "epsg:4326")

r2_ spatraster before projection is

> r2_
class       : SpatRaster 
dimensions  : 33, 48, 1  (nrow, ncol, nlyr)
resolution  : 2000, 2000  (x, y)
extent      : 9757195, 9853195, 734695, 800695  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / NSIDC EASE-Grid 2.0 Global (EPSG:6933) 
source      : memory 
name        : lyr.1 
min value   :     1 
max value   :     1

After projection, values are lost.

> r2_proj 
class       : SpatRaster 
dimensions  : 27, 52, 1  (nrow, ncol, nlyr)
resolution  : 0.01927436, 0.01927436  (x, y)
extent      : 101.1252, 102.1275, 5.768228, 6.288636  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : memory 
name        : lyr.1 
min value   :   NaN 
max value   :   NaN

This workflow has been tested for many datasets of points and extent, so this unexpected output seems to be generated by these values of points and extent.

When I set gdal to FALSE, I then obtain non-null values, hence it seems to result from the GDAL-warp algorithm.

terra::project(x = r2_,
               y = "epsg:4326", gdal = F)

> terra::project(x = r2_,
+                y = "epsg:4326", gdal = F)
class       : SpatRaster 
dimensions  : 27, 52, 1  (nrow, ncol, nlyr)
resolution  : 0.01927436, 0.01927436  (x, y)
extent      : 101.1252, 102.1275, 5.768228, 6.288636  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : memory 
name        : lyr.1 
min value   :   0.5 
max value   :   0.5
Solution

  • That is mysterious. Perhaps you can raise an issue.

    But for what you want to achieve, it would be better to do it like this anyway. That is, avoid transforming raster data, and rasterize to the raster you want.

    library(terra)
ext <-  ext( 9757195, 9853641, 734695, 799794  )
r <- rast(ext, resolution = 2000, crs = "EPSG:6933")
 
coord_vec <- data.frame( x = c(9849641, 9761195), y = c(795794.8, 738695.7))
coord_vec <- vect(coord_vec, crs="EPSG:6933", geom=c("x", "y"))

v <- project(coord_vec, "epsg:4326")
# project the raster template by using rast(r)
r2 <- project(rast(r), "epsg:4326")
r2 <- rasterize(v, r2)