Wrong mean for a raster with units of meters

[furu-taka02]

HI, I tested if "exactextract" worked correctly on Daymet data, which is a raster data with its unit of meters. I prepared a Shapefile of a polygon which covers exactly 2 Daymet cells as test data (the polygon in this Shapefile has WGS84 as its CRS), and compared the result of exactextract and that calculated manually (since it covers exactly 2 cells, I can compute the mean by extracting the precipitation values of these 2 cells). The true (manually computed) result is (10.0 + 10.0)/2 = 10.0 mm, but the result from exactextract is 29.16989 mm (whether or not I re-projected the Shapefile into the same units of meters through the CRS used in the Daymet raster data).

Is this a bug of "exactextractr"? Or is "exactextractr" not applicable to rasters in meters as its unit?

[dbaston]

Can you share the data used and command that you ran?

[furu-taka02]

I used Daymet dataset "daymet_v3_prcp_2017_na.nc4" obtained from https://daac.ornl.gov/cgi-bin/dsviewer.pl?ds_id=1328. (the size of the Daymet data is big, ~1.2GB)

The Shapefile I used can be downloaded from my Google Drive, https://drive.google.com/open?id=1hLXHUnSm9OKm0-Mc4akR4NLOsXw5FvVG. The Shapefile above is created as follows. The coordinates in meters of the 2 Daymet cells I used is (x, y) = (1419250, -795500), (1419250, -796500), (1420750, -796000), (1421250, -795500). This (x, y) coordinates are transformed into (longitude, latitude) coordinates by the CRS "Lambert Conformal Conic". You can see the detail of the CRS from the Daymet dataset I specified above, but for reference, I write the proj code of the CRS as follows: +proj=lcc +lat_1=25.0 +lat_2=60.0 +lon_0=-100.0 +lat_0=42.5 +x_0=0.0 +y_0=0.0 +ellps=WGS84 By the CRS written above, the original (x, y) coordinates are transformed to (longitude, latitude) = (-84.054301, 33.741580), (-84.056422, 33.732421), (-84.034499, 33.728874), (-84.032375, 33.738033) These (lon, lat) coordinates are used for making a polygon contained in the Shapefile above.

The R script I ran is as follows.

library(raster)
library(sf)
library(exactextractr)

test_data_path <- './Daymet_2cells_testData_for_exactextract.shp'
test_shapefile <- st_read(test_data_path)

test_shapefile <- st_transform(test_shapefile, 
                               crs = "+proj=lcc +lat_1=25.0 +lat_2=60.0 +lon_0=-100.0 +lat_0=42.5 +x_0=0.0 +y_0=0.0 +ellps=WGS84")

netcdf_filepath <- './daymet_v3_prcp_2017_na.nc4'
netcdf_variable_name <- 'prcp' #precipitation is 'prcp' in Daymet's netcdf file
DaymetPrecip_Data <- stack(netcdf_filepath, varname=netcdf_variable_name) #import multi-band netcdf file

Mean_Precip_jan1_2017 <- exact_extract(DaymetPrecip_Data[[1]], test_shapefile, 'mean')

print(Mean_Precip_jan1_2017)

The true precipitation value of the target 2 Daymet cells can be seen by the following Python code.

from netCDF4 import Dataset     # for reading netcdf file of Daymet data

daymet_ds = Dataset('./daymet_v3_prcp_2017_na.nc4', 'r') # reads the netCDF file

precp = daymet_ds.variables['prcp'][0,:,:] # get the precipitation raster data of the first day of 2017 (Jan 1)
print(precp[5780][5980])
print(precp[5780][5981])

As you can see in the Python code above, as the test Daymet cells, I used 2 cells with indices [5780][5980] and [5780][5981].

[dbaston]

It looks like the CRS of the Daymet file is not read correctly by the raster package. If you look at DaymetPrecip_Data you'll see the PROJ string:

+proj=lcc +lon_0=-100 +lat_0=42.5 +x_0=0 +y_0=0 +lat_1=25 +a=6378137 +rf=298.257232666016 

which has only one standard parallel (25 degrees), whereas the transformed shapefile has two (25 degrees and 60 degrees).

A workaround is to force the Daymet CRS to what it should be:

crs(DaymetPrecip_Data) <- crs(test_shapefile)

After which exact_extract(DaymetPrecip_Data[[1]], test_shapefile, 'mean') returns 9.999985.

[furu-taka02]

Thank you, I got the same result!