0

I am a beginner in R and have two rasters with different extents. I need to resample the rasters and I use terra package for this but the output is a raster with all NAs. Does anyone have any idea how to fix it? I found a similar post but the solution that was suggested for that problem does not work in my case.

Thanks

> rast2
class       : SpatRaster 
dimensions  : 276, 573, 3  (nrow, ncol, nlyr)
resolution  : 0.0833333, 0.08333333  (x, y)
extent      : -123.25, -75.50002, 26, 49  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
source      : memory 
names       : GID,         yi,         yr 
min values  :  NA,  0.8363878,  0.1122796 
max values  :  NA, 15.7337789, 13.2661868 
> proyield2
class       : SpatRaster 
dimensions  : 335, 527, 1  (nrow, ncol, nlyr)
resolution  : 0.09027645, 0.09027645  (x, y)
extent      : -24.27065, 23.30504, 0.291846, 30.53446  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
source      : memory 
name        : corn_relDiff_local 
min value   :          -26.83018 
max value   :           27.11203 
reyield <- resample(proyield2, rast2, method = 'bilinear')

> reyield
class       : SpatRaster 
dimensions  : 276, 573, 1  (nrow, ncol, nlyr)
resolution  : 0.0833333, 0.08333333  (x, y)
extent      : -123.25, -75.50002, 26, 49  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
source      : memory 
name        : corn_relDiff_local 
min value   :                NaN 
max value   :                NaN 
mbh
  • 3
  • 3
  • can you please provide a link to the original file? – Robert Hijmans Oct 19 '22 at 17:31
  • The original source of the data is: https://zenodo.org/record/3905243, I downloaded yieldDifferentials.zip that has several layers. yields.csv is another data that I have changed into raster and called it rast2. – mbh Oct 19 '22 at 18:30

1 Answers1

4

From the original file we get this

library(terra)
r <- rast("corn_relDiff_local.tif")
r
#class       : SpatRaster 
#dimensions  : 296, 470, 1  (nrow, ncol, nlyr)
#resolution  : 9900, 9900  (x, y)
#extent      : -2376000, 2277000, 257400, 3187800  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=aea +lat_0=0 +lon_0=0 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs 
#source      : corn_relDiff_local.tif 
#name        : corn_relDiff_local 
#min value   :          -28.11839 
#max value   :           28.08758 

p <- project(r, "+proj=longlat")
p
#class       : SpatRaster 
#dimensions  : 335, 527, 1  (nrow, ncol, nlyr)
#resolution  : 0.09027645, 0.09027645  (x, y)
#extent      : -24.27065, 23.30504, 0.291846, 30.53446  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#source(s)   : memory
#name        : corn_relDiff_local 
#min value   :          -26.83018 
#max value   :           27.11203 

Which puts the USA in West Africa.

w <- geodata::world(path=".")
plot(p); lines(w)

enter image description here

So the definition of the coordinate reference system of the original file is wrong.

That is quite clear from this

#coord. ref. : +proj=aea +lat_0=0 +lon_0=0 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs 

A more appropriate value for +lon_0=0 would be +lon_0=-96. Changing that would put the USA in Canada. By also changing the latitude of origin to 23, you get everything in the right place:

crs(r) <- "+proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"
p <- project(r, "+proj=longlat")

w <- geodata::world(path=".")
plot(p); lines(w)

enter image description here

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63