Error when trying to crop spatial extent using package sf

Question

I am trying to plot a geological map of Australia, however I only want the southeastern corner so want to crop the shape file. I am trying to do this using the sf package however I am getting an error when trying to crop it using this tutorial : https://www.r-bloggers.com/2019/04/zooming-in-on-maps-with-sf-and-ggplot2/

Anyone know where I am going wrong?


Error in s2_geography_from_wkb(x, oriented = oriented, check = check) : 
  Evaluation error: Found 321 features with invalid spherical geometry.
[106] Loop 0 is not valid: Edge 632 has duplicate vertex with edge 705
[2405] Loop 0 is not valid: Edge 44 has duplicate vertex with edge 10856
[2530] Loop 0 is not valid: Edge 1344 has duplicate vertex with edge 1369
[3321] Loop 0 is not valid: Edge 0 has duplicate vertex with edge 1035
[3425] Loop 0 is not valid: Edge 5841 has duplicate vertex with edge 6409
[3817] Loop 0 is not valid: Edge 699 has duplicate vertex with edge 711
[4207] Loop 0 is not valid: Edge 0 has duplicate vertex with edge 41
[7980] Loop 0 is not valid: Edge 697 has duplicate vertex with edge 801
[8806] Loop 0 is not valid: Edge 498 has duplicate vertex with edge 558
[9000] Loop 0 is not valid: Edge 7388 has duplicate vertex with edge 8363
...and 311 more.

### Code

library(sf)
library(ggplot2)

## Load sf

geol <- st_read("GeologicUnitPolygons1M.shp")

## Crop to only southeastern Australia

geol_cropped <- st_crop(geol, xmin = 138, xmax =155,ymin = -32, ymax = -43)

## plot 

ggplot() + 
  geom_sf(data = geol_cropped, fill = geol$LITHOLOGY) + 
  ggtitle("Geological Map") + 
  coord_sf()

I tried to dput the head of geol but it was too many characters - any idea how I can post a sample of the data?

The error messages are not due to your error; they come from spherical geometry faults in your shapefile. For possible approaches consider this answer https://stackoverflow.com/questions/68478179/how-to-resolve-spherical-geometry-failures-when-joining-spatial-data/ — Jindra Lacko, Aug 06 '21 at 10:44
Hi there Jindra, many thanks for your comment. I have tried the following : ```geol$geometry <- geol$geometry %>% s2::s2_rebuild() %>% sf::st_as_sfc()``` it runs for over half an hour and still hasn't finished, is this usual? Even the full shapefile only takes maybe a few minutes to load. — Sophie Williams, Aug 06 '21 at 11:27
If the shapefile is big and complex - and the fact that it takes minutes to load would indicate so - then the repairing of spherically faulty geometries can take long (like overnight long). On the other hand it is an one off time investment, you can easily save it afterwards and have it ready for future work. — Jindra Lacko, Aug 06 '21 at 12:25
you may want to try the other approach - turning the s2 processing off - first. It is a quick and dirty approach, may be deprecated in the future and does not resolve the underlying issue (it just reverts to a less sophisticated fallback method) but if you require just a quick plot it may do the trick — Jindra Lacko, Aug 06 '21 at 12:28

score 1 · Answer 1 · answered Aug 06 '21 at 11:42

I tried to reproduce the error with giscoR package but I couldn't. Maybe your shapefile is not valid. A couple of things I would try:

sf::st_make_valid(geol) before cropping.
If it just a zooming matter, you don't really need to crop. You can limit the map to the desired extent with coord_sf().

See a reprex, hope you find a solution:

library(ggplot2)
library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(giscoR)

# Cropping

# My shape comes to giscoR since I don't have access to the original shp

geol <- gisco_get_countries(
  country = "Australia",
  resolution = 1,
  epsg = 4326
)


# Try this

geol <- st_make_valid(geol)

geol_cropped <- st_crop(geol, xmin = 138, xmax = 155, ymin = -32, ymax = -43)
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries

## plot cropped

ggplot(geol_cropped) +
  geom_sf() +
  ggtitle("Geological Map cropped") +
  coord_sf()


# Another approach, no crop an limit the plot

ggplot(geol) +
  geom_sf() +
  ggtitle("Geological Map with limits") +
  coord_sf(
    xlim = c(138, 155),
    ylim = c(-32, -43)
  )

^{Created on 2021-08-06 by the reprex package (v2.0.1)}

Hi there, many thanks for your answer - i'll have a go. The shapefile is accessible from the Australian government so you would have thought it would have been valid, but yes it potentially could be the issue, I am not sure! Link is here to shp file https://data.gov.au/dataset/ds-dga-48fe9c9d-2f10-49d2-bd24-ac546662c4ec/details - it is quite large though — Sophie Williams, Aug 06 '21 at 11:48

Error when trying to crop spatial extent using package sf

1 Answers1