Change in raster spatial extent following projection

[MarianaDevault]

Hi all,

After downloading and merging 290 granules of MOD17A3Hv006 data from EARTHDATA for year 2010, I am now trying to handle these layers in terra.

Merged_Map
class       : SpatRaster 
dimensions  : 36000, 86400, 2  (nrow, ncol, nlyr)
resolution  : 463.3127, 463.3127  (x, y)
extent      : -20015109, 20015109, -6671703, 10007555  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs 
source      : All_granules_original_resolution.tif 
names       : Npp_500m, Npp_QC_500m 
min values  :   0.0000,           0 
max values  :   3.2766,         254

First thing I did was to redefine the extent of the merged product using

extent_model_equal_earth <- ext(-20001600, 19998400, -10000800, 9999200)
new_extent_NPP_500_m <- terra::extend(Merged_Map, extent_model_equal_earth)
new_extent_NPP_500_m
#class       : SpatRaster 
#dimensions  : 5697, 11394, 2  (nrow, ncol, nlyr)
#resolution  : 3510.546, 3510.546  (x, y)
#extent      : -20002474, 19996692, -10001362, 9998221  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#source(s)   : memory
#names       : Npp_500m, Npp_QC_500m 
#min values  :   0.0000,           0 
#max values  :   3.2766,         254

Because I don't neet such a fine resolution (500 m), I aggregated the new_extent_NPP_500_m raster (to reach 10 km x 10 km) doing

new_extent_NPP_500_m_10_km_Aggregated <- aggregate(new_extent_NPP_500_m, fact=20, fun=mean)
new_extent_NPP_500_m_10_km_Aggregated
#class       : SpatRaster 
#dimensions  : 2160, 4320, 2  (nrow, ncol, nlyr)
#resolution  : 9266.254, 9266.254  (x, y)
#extent      : -20015109, 20015109, -10007555, 10007555  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs 
#source(s)   : memory
#names       : Npp_500m, Npp_QC_500m 
#min values  :   0.0000,       0.565 
#max values  :   3.2766,     254.000

Now, I need to project it to Eckert IV, which is not giving me the correct product, given that I am getting the image below:

geo_proj="+proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m"
raster_model_eckert <- rast(xmin=-20001600, xmax=19998400, ymin=-10000800, ymax=9999200, nrows=2000, ncols=4000, crs=geo_proj)
raster_model_eckert
#class       : SpatRaster 
#dimensions  : 2000, 4000, 1  (nrow, ncol, nlyr)
#resolution  : 10000, 10000  (x, y)
#extent      : -20001600, 19998400, -10000800, 9999200  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

new_extent_NPP_10_km_Aggregated_Projected <- terra::project(new_extent_NPP_500_m_10_km_Aggregated, geo_proj, method="bilinear", mask=TRUE, overwrite=TRUE)
new_extent_NPP_10_km_Aggregated_Projected
#class       : SpatRaster 
#dimensions  : 4820, 9638, 2  (nrow, ncol, nlyr)
#resolution  : 3510.546, 3510.546  (x, y)
#extent      : -16916704, 16917943, -8460232, 8460601  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#source(s)   : memory
#names       : Npp_500m, Npp_QC_500m 
#min values  :   0.0000,   0.7023965 
#max values  :   3.2766, 254.0000000

Can anyone help me figure out why I am getting such a deformed map and why the raster extent has changed following projection?

[rhijmans]

This site is for bug report and feature request. For general questions, please use stackoverflow.com or gis.stackexchange.com.

But as you are here now: It is not clear why you are saying that the results are not good. The image is too small to be sure but it looks good to me. Is it because the blockiness in the ocean? That is probably because there were no tiles there (no land).

It seems to me that you take some unnecessary steps. Given SpatRaster x and template eck you can directly use project. Using extend is not going to change anything. If you first use aggregate that could improve things if you use "bilinear interpolation" in the projection; but if you use "average" you would be better off not aggregating.

x <- rast(nrow=36000, ncol=86400, nlyr=2, ext=ext(-20015109, 20015109, -6671703, 10007555), 
crs = "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs")

eck_proj="+proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m"
eck <- rast(xmin=-20001600, xmax=19998400, ymin=-10000800, ymax=9999200, nrows=2000, ncols=4000, crs=eck_proj)

y <- project(x, eck, "average")