3

From a set of 420 files (monthly) containing three columns (longitude, latitude, cluster) I have extracted shapefiles defined by cluster value. To do this I have read some posts in StackOverflow and Gis.Stackechange and used this R code

Then I have tried to use this code to plot shapefiles in monthly panels with this code but had received an error message from one of the shapefiles in the fortify step

      datoshp <- readOGR(dsn=".", layer=shpname[cluster])

      #fortify for ggplot mapping
      datoshp@data$id <- rownames(datoshp@data)
      datoshp.f <- fortify(datoshp, region = "id")

Full R code available here

The code has run fine for a number of files but it has crashed when reading the shapefile corresponding to cluster 4 in January 1993 with this error code in the fortify step

datoshp.f <- fortify(datoshp, region = "id") Error: TopologyException: Input geom 0 is invalid: Hole lies outside shell at or near point 11 35 at 11 35

I have found some info in other posts about problems with intersecting shapefiles and projections. But as the first set of shapefiles have run fine, I assume there is no problem with projection (maybe I'm wrong) as they were created from the same code.

Original lon-lat data and the shapefiles can be found in a Dropbox folder at https://www.dropbox.com/sh/e95gjjxu2563zvy/AADDypr810i0p1W2S8xpre1Za?dl=0

These are my first attempts with R and spatial data so I am not an expert about shapefiles or projections.

4
  • Try bringing the corrupted shapefile into QGIS to repair the geometry. This may be helpful: smsgis.co.za/using-qgis-fix-geometry-errors
    – Aaron
    Commented Feb 8, 2017 at 13:44
  • 1
    Can you edit down your question? You only need to show us what is going wrong with that one file. You should be able to create that error in about two or three lines of R code (readOGR, fortify). All the other detail is nice, but ultimately irrelevant and distracting.
    – Spacedman
    Commented Feb 8, 2017 at 13:45
  • @Spacedman Question edited. Sometimes it's difficult to know what would be relevant or not for the question. Hope it's better now.
    – pacomet
    Commented Feb 8, 2017 at 14:58
  • @Aaron I will try QGIS but I would prefer an R based anwer. I've never used QGIS. Thanks
    – pacomet
    Commented Feb 8, 2017 at 15:04

1 Answer 1

3

Your question boils down to getting this error with your fourth shapefile:

> library(raster); library(sp); library(rgdal)
> s4 = shapefile("SST-shape-199301-4.shp")
> s4$id = rownames(s4@data)
> f4 = fortify(s4, region="id")
Error: TopologyException: Input geom 0 is invalid: Hole lies outside shell at or near point 11 35 at 11 35

Interestingly you don't get it when you fortify without regions, even though its only one region:

> f4 = fortify(s4)
Regions defined for each Polygons
>

By plotting s4 and looking at coordinate 11 35 we can see something:

> plot(s4, col="red")
> points(11,35, cex=4)

enter image description here

Holes in polygons are okay, but they can cause problems if they go the wrong way round or if the hole is caused by the polygon looping itself. There's also something odd here:

enter image description here

which could be a number of things.

One way of fixing these kinds of errors in R is to use a zero-width buffer using gBuffer from the rgeos package.

> s4b = rgeos::gBuffer(s4, width=0)
Warning message:
In rgeos::gBuffer(s4, width = 0) :
  Spatial object is not projected; GEOS expects planar coordinates

> s4b$id = s4$id
> s4b$id
[1] 0
Levels: 0
> f4b = fortify(s4b, region="id")
> 

Ignore the warning message (which you might not get if you set the projection correctly) and you've now got something that fortifies without error.

If you plot(s4b, col="red") you'll notice that the second of the above anomalies has gone, but the first one (the hole) is still there.

The real solution of course if to find out why you got invalid geometries in the first place!

3
  • Tried your suggestion for another shapefile giving the same problem and works fine. Now it is launched for the whole bunch of files.
    – pacomet
    Commented Feb 9, 2017 at 8:36
  • To look for the real solution I have to rerun the clustering process, used CLARA algorithm for the clustering. Maybe there is an option to avoid isolated points inside a cluster, will give it a look.
    – pacomet
    Commented Feb 9, 2017 at 8:37
  • Nice trick with the 0-width buffer! (Though kind of crazy that there's an internal representation that's more powerful than the standard, and we only have indirect access to it via weird tricks).
    – mdsumner
    Commented Feb 9, 2017 at 10:47

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.