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 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 at 1:12
  • @ElioDiaz thanks so much for catching that. The 'my_ndvi_Raster.img' is actually a S4 RasterLayer. – ENIAC-6 Jan 6 at 1:41
  • is it working now? just answer your own question or add the fix to OP – Elio Diaz Jan 6 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 at 2:36
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, privacy policy and cookie policy

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