Unexpected `project` results #1356

elgabbas · 2023-11-27T18:12:24Z

Hello,

I would like to use terra::project() to project CHELSA data into a custom grid at EPSG:3035 projection. However, the values of some of the projected maps are not as expected. Here is an example:

download the tif file

require(dplyr)
require(raster)
require(terra)

# "https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V2/GLOBAL/climatologies/1981-2010/bio/CHELSA_gdd5_1981-2010_V.2.1.tif" %>%
#   download.file(destfile = "CHELSA_gdd5_1981-2010_V.2.1.tif", mode = "wb", quiet = TRUE)

reference grid [attached]

GridR <- terra::rast("Grid.tif")
plot(GridR)

reading as raster, then convert to spatraster

# No unscaling of the values (but okay for now)
Rstr1 <- "CHELSA_gdd5_1981-2010_V.2.1.tif" %>%
  raster::raster() %>%
  terra::rast()
plot(Rstr1)

plot(terra::crop(Rstr1, terra::ext(-30, 50, 35, 75)))

projecting to the reference grid

tictoc::tic()
Rstr1_proj <- terra::project(Rstr1, GridR, method = "average", threads = TRUE)
Rstr1_proj
tictoc::toc()
# 73.64 sec elapsed
# 
# class       : SpatRaster 
# dimensions  : 404, 469, 1  (nrow, ncol, nlyr)
# resolution  : 10000, 10000  (x, y)
# extent      : 2630000, 7320000, 1380000, 5420000  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs 
# source(s)   : memory
# varname     : Grid_10_Land_trimed 
# name        : CHELSA_gdd5_1981.2010_V.2.1 
# min value   :                        1.00 
# max value   :                    66754.57 
plot(Rstr1_proj)

slow, but values are within the expected range

reading directly as spatraster

Rstr2 <- terra::rast("CHELSA_gdd5_1981-2010_V.2.1.tif")
plot(Rstr2)

# comparable values 
plot(terra::crop(Rstr2, terra::ext(-30, 50, 35, 75)))

projecting to the reference grid

tictoc::tic()
Rstr2_proj <- terra::project(Rstr2, GridR, method = "average", threads = TRUE)
Rstr2_proj
tictoc::toc()
# 6.47 sec elapsed [much faster]
# 
# class       : SpatRaster 
# dimensions  : 404, 469, 1  (nrow, ncol, nlyr)
# resolution  : 10000, 10000  (x, y)
# extent      : 2630000, 7320000, 1380000, 5420000  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs 
# source(s)   : memory
# varname     : Grid_10_Land_trimed 
# name        : CHELSA_gdd5_1981-2010_V.2.1 
# min value   :                         0.1 
# max value   :                 214748364.8 
plot(Rstr2_proj)

The output is not expected

It seems that the NA values in the original tif file [NoData Value=2147483647 for CHELSA_gdd5_1981-2010_V.2.1.tif] are considered in the calculation of new values.

plot(terra::crop(Rstr2, terra::ext(-30, 50, 35, 75)), colNA = "red")

description of input files

describe("CHELSA_gdd5_1981-2010_V.2.1.tif")
 [1] "Driver: GTiff/GeoTIFF"                                                    
 [2] "Files: CHELSA_gdd5_1981-2010_V.2.1.tif"                                   
 [3] "Size is 43200, 20880"                                                     
 [4] "Coordinate System is:"                                                    
 [5] "GEOGCRS[\"WGS 84\","                                                      
 [6] "    DATUM[\"World Geodetic System 1984\","                                
 [7] "        ELLIPSOID[\"WGS 84\",6378137,298.257223563,"                      
 [8] "            LENGTHUNIT[\"metre\",1]]],"                                   
 [9] "    PRIMEM[\"Greenwich\",0,"                                              
[10] "        ANGLEUNIT[\"degree\",0.0174532925199433]],"                       
[11] "    CS[ellipsoidal,2],"                                                   
[12] "        AXIS[\"geodetic latitude (Lat)\",north,"                          
[13] "            ORDER[1],"                                                    
[14] "            ANGLEUNIT[\"degree\",0.0174532925199433]],"                   
[15] "        AXIS[\"geodetic longitude (Lon)\",east,"                          
[16] "            ORDER[2],"                                                    
[17] "            ANGLEUNIT[\"degree\",0.0174532925199433]],"                   
[18] "    USAGE["                                                               
[19] "        SCOPE[\"Horizontal component of 3D system.\"],"                   
[20] "        AREA[\"World.\"],"                                                
[21] "        BBOX[-90,-180,90,180]],"                                          
[22] "    ID[\"EPSG\",4326]]"                                                   
[23] "Data axis to CRS axis mapping: 2,1"                                       
[24] "Origin = (-180.000138888850017,83.999860415149996)"                       
[25] "Pixel Size = (0.008333333300000,-0.008333333300000)"                      
[26] "Metadata:"                                                                
[27] "  AREA_OR_POINT=Area"                                                     
[28] "Image Structure Metadata:"                                                
[29] "  COMPRESSION=DEFLATE"                                                    
[30] "  INTERLEAVE=BAND"                                                        
[31] "  PREDICTOR=2"                                                            
[32] "Corner Coordinates:"                                                      
[33] "Upper Left  (-180.0001389,  83.9998604) (180d 0' 0.50\"W, 83d59'59.50\"N)"
[34] "Lower Left  (-180.0001389, -90.0001389) (180d 0' 0.50\"W, 90d 0' 0.50\"S)"
[35] "Upper Right ( 179.9998597,  83.9998604) (179d59'59.49\"E, 83d59'59.50\"N)"
[36] "Lower Right ( 179.9998597, -90.0001389) (179d59'59.49\"E, 90d 0' 0.50\"S)"
[37] "Center      (  -0.0001396,  -3.0001392) (  0d 0' 0.50\"W,  3d 0' 0.50\"S)"
[38] "Band 1 Block=43200x1 Type=Int32, ColorInterp=Gray"                        
[39] "  NoData Value=2147483647"                                                
[40] "  Offset: 0,   Scale:0.1"

describe("Grid.tif")
 [1] "Driver: GTiff/GeoTIFF"                                                    
 [2] "Files: D:/BioDT_IAS/Data/Maps/Grid/Grid_10_Land_trimed.tif"               
 [3] "Size is 469, 404"                                                         
 [4] "Coordinate System is:"                                                    
 [5] "PROJCRS[\"unknown\","                                                     
 [6] "    BASEGEOGCRS[\"unknown\","                                             
 [7] "        DATUM[\"Unknown based on GRS80 ellipsoid\","                      
 [8] "            ELLIPSOID[\"GRS 1980\",6378137,298.257222101004,"             
 [9] "                LENGTHUNIT[\"metre\",1],"                                 
[10] "                ID[\"EPSG\",7019]]],"                                     
[11] "        PRIMEM[\"Greenwich\",0,"                                          
[12] "            ANGLEUNIT[\"degree\",0.0174532925199433,"                     
[13] "                ID[\"EPSG\",9122]]]],"                                    
[14] "    CONVERSION[\"Lambert Azimuthal Equal Area\","                         
[15] "        METHOD[\"Lambert Azimuthal Equal Area\","                         
[16] "            ID[\"EPSG\",9820]],"                                          
[17] "        PARAMETER[\"Latitude of natural origin\",52,"                     
[18] "            ANGLEUNIT[\"degree\",0.0174532925199433],"                    
[19] "            ID[\"EPSG\",8801]],"                                          
[20] "        PARAMETER[\"Longitude of natural origin\",10,"                    
[21] "            ANGLEUNIT[\"degree\",0.0174532925199433],"                    
[22] "            ID[\"EPSG\",8802]],"                                          
[23] "        PARAMETER[\"False easting\",4321000,"                             
[24] "            LENGTHUNIT[\"metre\",1],"                                     
[25] "            ID[\"EPSG\",8806]],"                                          
[26] "        PARAMETER[\"False northing\",3210000,"                            
[27] "            LENGTHUNIT[\"metre\",1],"                                     
[28] "            ID[\"EPSG\",8807]]],"                                         
[29] "    CS[Cartesian,2],"                                                     
[30] "        AXIS[\"(E)\",east,"                                               
[31] "            ORDER[1],"                                                    
[32] "            LENGTHUNIT[\"metre\",1]],"                                    
[33] "        AXIS[\"(N)\",north,"                                              
[34] "            ORDER[2],"                                                    
[35] "            LENGTHUNIT[\"metre\",1]]]"                                    
[36] "Data axis to CRS axis mapping: 1,2"                                       
[37] "Origin = (2630000.000000000000000,5420000.000000000000000)"               
[38] "Pixel Size = (10000.000000000000000,-10000.000000000000000)"              
[39] "Metadata:"                                                                
[40] "  AREA_OR_POINT=Area"                                                     
[41] "Image Structure Metadata:"                                                
[42] "  COMPRESSION=LZW"                                                        
[43] "  INTERLEAVE=BAND"                                                        
[44] "Corner Coordinates:"                                                      
[45] "Upper Left  ( 2630000.000, 5420000.000) ( 31d33'21.33\"W, 67d 7' 8.99\"N)"
[46] "Lower Left  ( 2630000.000, 1380000.000) (  8d11'22.24\"W, 33d39'30.22\"N)"
[47] "Upper Right ( 7320000.000, 5420000.000) ( 70d52' 4.44\"E, 59d 6'23.44\"N)"
[48] "Lower Right ( 7320000.000, 1380000.000) ( 41d24'16.96\"E, 29d53' 3.45\"N)"
[49] "Center      ( 4975000.000, 3400000.000) ( 19d50'43.02\"E, 53d18'28.11\"N)"
[50] "Band 1 Block=469x4 Type=Float32, ColorInterp=Gray"                        
[51] "  Min=1.000 Max=1.000 "                                                   
[52] "  Minimum=1.000, Maximum=1.000, Mean=-9999.000, StdDev=-9999.000"         
[53] "  NoData Value=-3.4e+38"                                                  
[54] "  Metadata:"                                                              
[55] "    STATISTICS_MAXIMUM=1"                                                 
[56] "    STATISTICS_MEAN=-9999"                                                
[57] "    STATISTICS_MINIMUM=1"                                                 
[58] "    STATISTICS_STDDEV=-9999"

Session info

sessioninfo::session_info()
─ Session info ───────────────────────────────────────────────────────────
 setting  value
 version  R version 4.3.2 (2023-10-31 ucrt)
 os       Windows 10 x64 (build 19045)
 system   x86_64, mingw32
 ui       RStudio
 language (EN)
 collate  English_United States.utf8
 ctype    English_United States.utf8
 tz       Europe/Berlin
 date     2023-11-27
 rstudio  2023.09.1+494 Desert Sunflower (desktop)
 pandoc   NA

─ Packages ───────────────────────────────────────────────────────────────
 package     * version date (UTC) lib source
 class         7.3-22  2023-05-03 [1] CRAN (R 4.3.0)
 classInt      0.4-10  2023-09-05 [1] CRAN (R 4.3.2)
 cli           3.6.1   2023-03-23 [1] CRAN (R 4.3.0)
 codetools     0.2-19  2023-02-01 [2] CRAN (R 4.3.2)
 DBI           1.1.3   2022-06-18 [1] CRAN (R 4.3.0)
 dplyr       * 1.1.3   2023-09-03 [1] CRAN (R 4.3.2)
 e1071         1.7-13  2023-02-01 [1] CRAN (R 4.3.0)
 fansi         1.0.5   2023-10-08 [1] CRAN (R 4.3.2)
 generics      0.1.3   2022-07-05 [1] CRAN (R 4.3.0)
 glue          1.6.2   2022-02-24 [1] CRAN (R 4.3.0)
 KernSmooth    2.23-22 2023-07-10 [1] CRAN (R 4.3.1)
 lattice       0.21-9  2023-10-01 [2] CRAN (R 4.3.2)
 lifecycle     1.0.4   2023-11-07 [1] CRAN (R 4.3.2)
 magrittr      2.0.3   2022-03-30 [1] CRAN (R 4.3.0)
 pillar        1.9.0   2023-03-22 [1] CRAN (R 4.3.0)
 pkgconfig     2.0.3   2019-09-22 [1] CRAN (R 4.3.0)
 proxy         0.4-27  2022-06-09 [1] CRAN (R 4.3.0)
 R6            2.5.1   2021-08-19 [1] CRAN (R 4.3.0)
 raster      * 3.6-26  2023-10-14 [1] CRAN (R 4.3.2)
 Rcpp          1.0.11  2023-07-06 [1] CRAN (R 4.3.1)
 rlang         1.1.2   2023-11-04 [1] CRAN (R 4.3.2)
 rstudioapi    0.15.0  2023-07-07 [1] CRAN (R 4.3.1)
 sessioninfo   1.2.2   2021-12-06 [1] CRAN (R 4.3.0)
 sf            1.0-14  2023-07-11 [1] CRAN (R 4.3.1)
 sp          * 2.1-1   2023-10-16 [1] CRAN (R 4.3.2)
 terra       * 1.7-55  2023-10-13 [1] CRAN (R 4.3.2)
 tibble        3.2.1   2023-03-20 [1] CRAN (R 4.3.0)
 tictoc        1.2     2023-04-23 [1] CRAN (R 4.3.0)
 tidyselect    1.2.0   2022-10-10 [1] CRAN (R 4.3.0)
 units         0.8-4   2023-09-13 [1] CRAN (R 4.3.2)
 utf8          1.2.4   2023-10-22 [1] CRAN (R 4.3.2)
 vctrs         0.6.4   2023-10-12 [1] CRAN (R 4.3.2)

 [1] C:/Users/*****/AppData/Local/R/win-library/4.3
 [2] C:/Program Files/R/R-4.3.2/library

reference grid

Thanks

The text was updated successfully, but these errors were encountered:

rhijmans · 2023-11-27T18:48:02Z

Here is a simplified version of your example that does not require an additional download.

library(terra)
url <- "https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V2/GLOBAL/climatologies/1981-2010/bio/CHELSA_gdd5_1981-2010_V.2.1.tif"
f <- basename(url)
if (!file.exists(f)) download.file(url, f, mode = "wb")

# Reference grid created with `as.character(GridR)`
GridR <- rast(ncols=469, nrows=404, nlyrs=1, xmin=2630000, xmax=7320000, ymin=1380000, ymax=5420000, names=c('Grid'), crs='PROJCRS[\"unknown\",BASEGEOGCRS[\"unknown\",DATUM[\"Unknown based on GRS80 ellipsoid\",ELLIPSOID[\"GRS 1980\",6378137,298.257222101004,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",7019]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]]],CONVERSION[\"Lambert Azimuthal Equal Area\",METHOD[\"Lambert Azimuthal Equal Area\",ID[\"EPSG\",9820]],PARAMETER[\"Latitude of natural origin\",52,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",10,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"False easting\",4321000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",3210000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1]]]')

And the projections

# bad
r1 <- rast(f)
r1_p <- project(r1, GridR, method = "average", threads = TRUE)

# good if first reading all values
r1a <- rast(f) * 1
r1a_p <- project(r1a, GridR, method = "average", threads = TRUE)

# good when coercing from a RasterLayer
r2 <- rast(raster::raster(f))
r2_p <- terra::project(r2, GridR, method = "average", threads = TRUE)

This is mysterious. Especially that it works when coercing from a RasterLayer! I will have a look.

mdsumner · 2023-11-30T19:41:54Z

fwiw, looks ok through 'GDALWarp()` api (with 'by_util'), also no need to download with the /vsicurl protocol

library(terra)
url <- "https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V2/GLOBAL/climatologies/1981-2010/bio/CHELSA_gdd5_1981-2010_V.2.1.tif"

f <- file.path("/vsicurl", url)

GridR <- rast(ncols=469, nrows=404, nlyrs=1, xmin=2630000, xmax=7320000, ymin=1380000, ymax=5420000, names=c('Grid'), crs='EPSG:3035')

r1 <- rast(f)
r1_p <- project(r1, GridR, method = "average", threads = TRUE, by_util = TRUE)

elgabbas · 2023-11-30T23:17:40Z

Thanks @mdsumner and @rhijmans for your feedback. This seems to be a temporary solution to avoid extremely high values shown in the above example.

However, there are value differences when coercing raster to spatraster vs using by_util = TRUE.

# raster coerced to spatraster
R1 <- terra::rast(raster::raster(f)/10)
(R1_p <- terra::project(R1, GridR, method = "average", threads = TRUE))
# class       : SpatRaster 
# dimensions  : 404, 469, 1  (nrow, ncol, nlyr)
# resolution  : 10000, 10000  (x, y)
# extent      : 2630000, 7320000, 1380000, 5420000  (xmin, xmax, ymin, ymax)
# coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035) 
# source(s)   : memory
# name        : CHELSA_gdd5_1981.2010_V.2.1 
# min value   :                       0.100 
# max value   :                    6675.457

# read as spatraster, `by_util = TRUE`
R2 <- rast(f)
(R2_p <- terra::project(R2, GridR, method = "average", threads = TRUE, by_util = TRUE))
# class       : SpatRaster 
# dimensions  : 404, 469, 1  (nrow, ncol, nlyr)
# resolution  : 10000, 10000  (x, y)
# extent      : 2630000, 7320000, 1380000, 5420000  (xmin, xmax, ymin, ymax)
# coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035) 
# source(s)   : memory
# name        : CHELSA_gdd5_1981-2010_V.2.1 
# min value   :                         0.0 
# max value   :                      6707.1

The difference between both approaches:

R1_p - R2_p
# class       : SpatRaster 
# dimensions  : 404, 469, 1  (nrow, ncol, nlyr)
# resolution  : 10000, 10000  (x, y)
# extent      : 2630000, 7320000, 1380000, 5420000  (xmin, xmax, ymin, ymax)
# coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035) 
# source(s)   : memory
# name        : CHELSA_gdd5_1981.2010_V.2.1 
# min value   :                   -1230.521 
# max value   :                    1370.416

Please note differences in the pattern of NA values (e.g. north Italy and Iceland) when using by_util = TRUE and also the value range of both maps

plot(c(R1_p, R2_p, R1_p - R2_p), axes = F, colNA = "blue")

Any explanation/suggestion? Thanks..

rhijmans · 2023-12-04T04:05:37Z

This bug occurs because the input has a scale that is 0.1 (not 1). This is why if you first make a copy it works (that removes the scale). With a SpatRaster created from a RasterLayer, a copy is first made by project. That is why it works.

elgabbas · 2023-12-04T21:58:12Z

Thanks @rhijmans for your explanation.

For another CHELSA map (which also has a scale factor of 0.1), this problem did not occur.

library(terra)

url <- "https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V2/GLOBAL/climatologies/1981-2010/bio/CHELSA_bio12_1981-2010_V.2.1.tif"
f <- basename(url)
if (!file.exists(f)) download.file(url, f, mode = "wb")

# Reference grid created with `as.character(GridR)`
GridR <- rast(ncols=469, nrows=404, nlyrs=1, xmin=2630000, xmax=7320000, ymin=1380000, ymax=5420000, names=c('Grid'), crs='PROJCRS[\"unknown\",BASEGEOGCRS[\"unknown\",DATUM[\"Unknown based on GRS80 ellipsoid\",ELLIPSOID[\"GRS 1980\",6378137,298.257222101004,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",7019]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]]],CONVERSION[\"Lambert Azimuthal Equal Area\",METHOD[\"Lambert Azimuthal Equal Area\",ID[\"EPSG\",9820]],PARAMETER[\"Latitude of natural origin\",52,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",10,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"False easting\",4321000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",3210000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1]]]')

r1 <- rast(f)
r1_p <- project(r1, GridR, method = "average", threads = TRUE)

r1a <- rast(f) * 1
r1a_p <- project(r1a, GridR, method = "average", threads = TRUE)

r2 <- rast(raster::raster(f))
r2_p <- terra::project(r2, GridR, method = "average", threads = TRUE) * 0.1

plot(c(r1_p, r1a_p, r2_p))

The difference between both files is mainly in the NA value (using terra::describe()). For the CHELSA_gdd5_1981-2010_V.2.1.tif the NoData value is 2147483647, while for the CHELSA_bio12_1981-2010_V.2.1.tif the NoData value is -99999".

In the first example, the problematic Rstr2_proj object has a maximum value of 214748364.8 which is very similar to the NoData value (after unscaling). Do you think that grid cells with NoData are considered in calculating the new values (i.e., NA cells were considered as 2147483647 in the calculation of the average)?

Among all current bio variables in CHLESA, extremely high values in projected maps existed in the following three maps (all having a NoData value of 2147483647)

https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V2/GLOBAL/climatologies/1981-2010/bio/CHELSA_gdd5_1981-2010_V.2.1.tif
https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V2/GLOBAL/climatologies/1981-2010/bio/CHELSA_gdd10_1981-2010_V.2.1.tif
https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V2/GLOBAL/climatologies/1981-2010/bio/CHELSA_gsp_1981-2010_V.2.1.tif

Thanks

elgabbas · 2023-12-04T22:55:57Z

The following example suggests that the problem is how the NoData is treated.

library(terra)
library(stringr)
library(dplyr)

# download file
url <- "https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V2/GLOBAL/climatologies/1981-2010/bio/CHELSA_gdd5_1981-2010_V.2.1.tif"
f <- basename(url)
if (!file.exists(f)) download.file(url, f, mode = "wb")

# Reference grid
GridR <- rast(ncols=469, nrows=404, nlyrs=1, xmin=2630000, xmax=7320000, ymin=1380000, ymax=5420000, names=c('Grid'), crs='PROJCRS[\"unknown\",BASEGEOGCRS[\"unknown\",DATUM[\"Unknown based on GRS80 ellipsoid\",ELLIPSOID[\"GRS 1980\",6378137,298.257222101004,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",7019]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]]],CONVERSION[\"Lambert Azimuthal Equal Area\",METHOD[\"Lambert Azimuthal Equal Area\",ID[\"EPSG\",9820]],PARAMETER[\"Latitude of natural origin\",52,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",10,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"False easting\",4321000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",3210000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1]]]')

# read the tif file, then save it as a different name
f_Edited <- stringr::str_replace(f, ".tif", "_edited.tif")

f %>% 
  terra::rast() %>% 
  # write the raster as a tif file, keep scale = 0.1 and offset = 0
  terra::writeRaster(
    filename = f_Edited, NAflag = -99999, overwrite = TRUE, progress = TRUE, scale = 0.1, offset = 0)
# output file size is larger, but this is okay

Now, projecting the original and edited tiff file

# This is problematic as before
r1 <- rast(f)
r1_p <- project(r1, GridR, method = "average", threads = TRUE)

# This works as expected
r3 <- rast(f_Edited)
r3_p <- project(r3, GridR, method = "average", threads = TRUE)

plot(c(r1_p, r3_p))

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Unexpected `project` results #1356

Unexpected `project` results #1356

elgabbas commented Nov 27, 2023 •

edited

Loading

rhijmans commented Nov 27, 2023 •

edited

Loading

mdsumner commented Nov 30, 2023

elgabbas commented Nov 30, 2023 •

edited

Loading

rhijmans commented Dec 4, 2023

elgabbas commented Dec 4, 2023 •

edited

Loading

elgabbas commented Dec 4, 2023 •

edited

Loading

Unexpected project results #1356

Unexpected project results #1356

Comments

elgabbas commented Nov 27, 2023 • edited Loading

download the tif file

reference grid [attached]

reading as raster, then convert to spatraster

projecting to the reference grid

reading directly as spatraster

projecting to the reference grid

description of input files

Session info

rhijmans commented Nov 27, 2023 • edited Loading

mdsumner commented Nov 30, 2023

elgabbas commented Nov 30, 2023 • edited Loading

rhijmans commented Dec 4, 2023

elgabbas commented Dec 4, 2023 • edited Loading

elgabbas commented Dec 4, 2023 • edited Loading

Unexpected `project` results #1356

Unexpected `project` results #1356

elgabbas commented Nov 27, 2023 •

edited

Loading

rhijmans commented Nov 27, 2023 •

edited

Loading

elgabbas commented Nov 30, 2023 •

edited

Loading

elgabbas commented Dec 4, 2023 •

edited

Loading

elgabbas commented Dec 4, 2023 •

edited

Loading