search

rspatial / terra

R package for spatial data handling https://rspatial.github.io/terra/reference/terra-package.html
GNU General Public License v3.0
543 stars 90 forks source link

Unexpected `project` results #1356

Open elgabbas opened 11 months ago

elgabbas commented 11 months ago

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)

image

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)

image

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

image

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 image

reading directly as spatraster

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

image

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

image

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 image

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")

image

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

rhijmans commented 11 months ago

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 commented 11 months ago

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 commented 11 months ago

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")

image

Any explanation/suggestion? Thanks..

rhijmans commented 11 months ago

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 commented 11 months ago

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))

image

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 commented 11 months ago

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))

image