CODEWITHSUNDEEP
rrasterterra
dauby gilles
dauby gilles

Reputation: 67

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

Upvotes: 2

Views: 876

Answers (2)

Eric Lino
Eric Lino

Reputation: 479

Whenever you are reprojecting a raster, there is a resampling process linked to it. When you have the argument 'GDAL=TRUE', it uses the Nearest Neighbor method by default, so it can make sense if your raster is mostly composed by NaNs and few 1s. Alternatively, when you have 'GDAL=FALSE', it uses a bilinear interpolation. If you want to keep the 1s, you could consider using 'maximum' method.

Upvotes: -1

Robert Hijmans
Robert Hijmans

Reputation: 47026

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)

Upvotes: 1

Related Questions