0

I used polygon to extract the nc file, which is the temperature data of ERA5, but when using stack and using the exact_extract function, there were warnings and errors.

> rasdata <- stack("./data/era5.temperature.19590106.nc")
Warning message:
In .rasterObjectFromCDF(x, type = objecttype, band = band, ...) :
  "level" set to 1 (there are 19 levels)
> extrad <- exact_extract(rasdata,map,progress = FALSE) 
Error: [names<-] incorrect number of names

The era5 file can be download in https://drive.google.com/file/d/1IPwNuIa7Lk1FTAxTgLX7A3m4UsM2KRZc/view?usp=sharing

Improve this question
2
  • What's your map object here?
    – Spacedman
    Commented Nov 6, 2022 at 9:59
  • Are you adding the ncdf4 library before reading the NetCDF file? Why are you ignoring the warning coming from raster::stack? Your problem is not with exact_extract but, coming from an incorrect read of the source file. Admittedly, it is not clear that ncdf4 is required for correct import of cdf files but this is, in fact, the case. Alternately, you could try this using the terra library with terra::rast to see if this mitigates the cdf import.
    – Jeffrey Evans
    Commented Nov 6, 2022 at 16:40

1 Answer 1

Reset to default
-1

exactextractr::exact_extract() requires x to be of type RasterLayer, RasterStack, RasterBrick or SpatRaster and y of type sf, sfc, SpatialPolygonsDataFrame or SpatialPolygons. Basically, you can use {raster} and {terra} for raster data input and {sf} and {sp} for vector data.

Just like this:

library(terra)
#> terra 1.6.33
library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE
library(geodata)
library(exactextractr)

# read raster
r <- rast("era5.temperature.19590106.nc")
r
#> class       : SpatRaster 
#> dimensions  : 241, 281, 456  (nrow, ncol, nlyr)
#> resolution  : 0.25, 0.25  (x, y)
#> extent      : 69.875, 140.125, -0.125, 60.125  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> source      : era5.temperature.19590106.nc 
#> varname     : t (Temperature) 
#> names       : t_lev~100_1, t_lev~150_1, t_lev~200_1, t_lev~250_1, t_lev~300_1, t_lev~350_1, ... 
#> unit        :           K,           K,           K,           K,           K,           K, ... 
#> time        : 1959-01-06 00:00:00 to 1959-01-06 23:00:00 UTC

# subset to one layer, just for this demo
r_sub <- r["t_level=100_10"]

# get some random polygon vector for demonstration purposes
adm <- gadm("Cambodia", level = 0, path = tempdir()) |> sf::st_as_sf()
adm
#> Simple feature collection with 1 feature and 2 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 102.3338 ymin: 9.91361 xmax: 107.6277 ymax: 14.69027
#> Geodetic CRS:  WGS 84
#>   GID_0  COUNTRY                       geometry
#> 1   KHM Cambodia MULTIPOLYGON (((104.5359 10...

# extract and inspect
vals <- exact_extract(r_sub, adm)

head(vals[[1]])
#>      value coverage_fraction
#> 1 193.1998      0.0369466245
#> 2 193.1374      0.1065797433
#> 3 193.2218      0.0737281293
#> 4 193.3062      0.0003664213
#> 5 193.3852      0.0259982608
#> 6 193.4476      0.0130981132

Created on 2022-11-06 with reprex v2.0.2

Improve this answer
3
  • Sorry, this just does not address the question. The op's code clearly indicates that the object is being read in as a raster stack object, which is returning a warning. If map was not the correct object class the error would clearly indicate this. Also, when providing a worked example you do not need to show object details, it makes the answer much less readable. You can accomplish much the same thing by simply returning the object so this information is show when working through the example eg., ( r <- rast("era5.temperature.19590106.nc") )
    – Jeffrey Evans
    Commented Nov 6, 2022 at 16:47
  • Thanks for your feedback, Jeffrey! I'd say this depends. Indeed, I blatantly ignored the fact that there were some issues when importing the file with {raster}, BUT I was also able to show that there are no general problems with exact_extract() when accessing ERA5 netcdf data (and you might assume this is the main concern here).
    – dimfalk
    Commented Nov 6, 2022 at 17:59
  • 1
    Also, you're right, I do not need to show the content of r and adm to make this work, but I might show from my point of view. This makes the example maybe a little less readable, but also a lot more transparent since you do not need to execute this yourself anymore to see the details, e.g. when you are reading this from your phone. Matter of taste, I'd say.
    – dimfalk
    Commented Nov 6, 2022 at 18:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.