0

I am using raster package erase function as per my previous post solution for clipping and dissolving overlapping polygons - Dissolve Overlapping Polygons using difference and union in R

For some of the polygons I am getting below error with erase funtion:

Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, unaryUnion_if_byid_false, : TopologyException: Input geom 1 is invalid: Self-intersection at or near point 1.1197332302192855 47.203098020153668 at 1.1197332302192855 47.203098020153668

library(raster)
library(rgeos)
library(sp)

fields <- gBuffer(fields, byid=TRUE, width=0) # Expands the given geometry to include 
the area within the specified width 

zone <- fields[fields$Type == "Zone", ]
plot <- fields[fields$Type == "Plot", ]

d <- erase(zone, plot) #issue here
spplot(d, "Rx")

# I tried using rgeos::gBuffer to avoid RGEOSBinTopology Exception but it did not worked out. Any guidance in this area would be really helpful.

zone <- gBuffer(zone, byid=TRUE, width=0)
plot <- gBuffer(plot, byid=TRUE, width=0)
string
  • 787
  • 10
  • 39
  • The error message is pretty clear. Make a plot of that location. Then fix it... If you do not know use, `crop` to select just that area so that you can share the data. – Robert Hijmans Mar 16 '20 at 04:03
  • Thank you Robert for the reply .... You can download the test shapefile from here - http://www.filedropper.com/trialmap1 – string Mar 16 '20 at 08:43
  • @RobertHijmans Did you get the chance to look into the same .... any possible solution... – string Mar 17 '20 at 20:39
  • I do not have access to the file. Preferably you would not use a file, but share the data as code (see e.g. `geom`). But if that is too difficult perhaps use google drive? Filedropper does not work for me --- just advertisements – Robert Hijmans Mar 19 '20 at 00:34

1 Answers1

0

The plot polygons are really messy. Even though they are valid. The problems arises in erase becomes it aggregates ("dissolves") argument y ("fields") but that generates an invalid topology. I need to look into how to catch that best, but, for now (and perhaps for a while), here is a work-around using cleangeo, but you might want to use a tool that cleans the fields polygons (snap vertices) --- perhaps with Rmapshaper?

library(raster)
library(rgeos)

f <- "WUERZ-I-WW _ Plot _TrialMap.shp"
fields = shapefile(f)
zone <- fields[fields$Type == "Zone", ]
plot <- fields[fields$Type == "Plot", ]

aplot <- aggregate(plot)
rgeos::gIsValid(aplot)
#[1] FALSE
#Warning message:
#In RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid") :
#  Self-intersection at or near point 13.033265979999999 51.190776640000003

library(cleangeo)
cplot <- clgeo_Clean(aplot)
rgeos::gIsValid(cplot)
# TRUE

d <- erase(zone, cplot)
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • Thank you for suggesting workaround.... it certainly helps, as I need to automate this feature using R Script so could not opt for a tool perhaps I will go through rmapshaper r package in detail. I am planning to post another question related converting shape file to ISO XML format by first converting the shape file to GeoJSON format - Could you please guide me if R Script / R Packages related Shape files could be useful in converting shape file or GeoJSON to ISO XML format - https://www.agrarheute.com/technik/traktoren/iso-xml-Datei-erklaert-563240 – string Mar 21 '20 at 19:11