1

A previous script that I had is no longer working after the sf package major change. I am using sf 0.9.6 and raster 3.0.2. I have tried to modify my code but cannot run raster::extract().

My previous code now breaks on the second command and does not generate extent_program. I receive: Error in sp::CRS(SRS_string = from$wkt): unused argument (SRS_string = from$wkt).

# Read in polygons from 'my_shapefile.shp'

sf_program_fallow <- 'my_shapefile.shp' %>% sf::st_read() %>% sf::st_transform(crs = 32611) 

# Generate extent of polygons (now generates error)
extent_program <- extent(sf::as_Spatial(sf_program_fallow))

# Calculate zonal statistics (note: 'my_ndvi_RasterLayer' created earlier)
ndvi_fallowed_fields <- raster::extract(x = 'my_ndvi_RasterLayer', 
                       y = sf_program_fallow,
                       method = 'simple',
                       fun = mean,
                       na.rm = T,
                       sp = TRUE)

Updating the first line of code to the following produces extent_program but not ndvi_fallowed_fields. I get the same error: Error in sp::CRS(SRS_string = from$wkt): unused argument (SRS_string = from$wkt)

sf_program_fallow <- 'my_shapefile.shp' %>% sf::st_read() %>% sf::st_transform(., sf::st_crs(32611))
5
  • "my_shapefile.shp" should be string, with quotes, unless you define it as my_shapefile.shp = "my_shapefile.shp"
    – Elio Diaz
    Jan 5, 2021 at 23:40
  • your argument sp = TRUE, according to help is for spatial* objects, sf_program_fallow is an sf object; I'm not sure .img files have crs info and, also according to extract manual, x should be a raster object, not sure if it can be read from disk straight into the function
    – Elio Diaz
    Jan 6, 2021 at 1:12
  • @ElioDiaz thanks so much for catching that. The 'my_ndvi_Raster.img' is actually a S4 RasterLayer.
    – ENIAC-6
    Jan 6, 2021 at 1:41
  • is it working now? just answer your own question or add the fix to OP
    – Elio Diaz
    Jan 6, 2021 at 2:00
  • @ElioDiaz, thanks, the code above is a truncated version of a much longer script so I haven't actually changed anything. The issue I believe lies in the recent sf change regarding coordinate reference systems.
    – ENIAC-6
    Jan 6, 2021 at 2:36

1 Answer 1

2

Two things solved this issue:

  1. Updated all packages following the instructions from: https://inbo.github.io/tutorials/tutorials/spatial_crs_coding/#defining-a-crs-with-sp
  • Windows users: to make sure that the rgdal and sf packages use GDAL 3 / PROJ ≥ 6 – which is now highly advised: please update all your packages in a pure R console (not RStudio) that uses your latest installed R version, by running update.packages(ask = FALSE, checkBuilt = TRUE).
  1. Following instructions from https://cran.r-project.org/web/packages/rgdal/vignettes/PROJ6_GDAL3.html#Migration_to_PROJ6GDAL3 I replaced crs = 32611 with SRS_string = "EPSG:32611".

The following is final code that ran:

# Updated section:
crs_wgs84utm11 <- sp::CRS(SRS_string = "EPSG:32611")
sf_program_fallow <- 'my_shapefile.shp' %>% sf::st_read() %>% sf::st_transform(crs = crs_wgs84utm11)

# unchanged:
extent_program <- extent(sf::as_Spatial(sf_program_fallow))

# Calculate zonal statistics (note: 'my_ndvi_RasterLayer' created earlier)
ndvi_fallowed_fields <- raster::extract(x = 'my_ndvi_RasterLayer', 
                   y = sf_program_fallow,
                   method = 'simple',
                   fun = mean,
                   na.rm = T,
                   sp = TRUE)

Your Answer

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

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