CODEWITHSUNDEEP
rggplot2geospatialshapefiler-sf
Sophie Williams
Sophie Williams

Reputation: 348

Error when trying to crop spatial extent using package sf

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?

Upvotes: 1

Views: 2261

Answers (1)

dieghernan
dieghernan

Reputation: 3402

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:

  1. sf::st_make_valid(geol) before cropping.
  2. 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)

Upvotes: 1

Related Questions