-1

Have the shape file which has multiple polygons with logical division of Zones and Plots. Plots are overlapping over Zones. The task is to dissolve / merge plots with Zones with no overlapping.

enter image description here Here is spplot of the shape file. Here plots are on top of the Field Zones. Also here is the shapefile with overlapped polygons (Zones & Plots): Shapefile

In QGIS, the same was achieved using Extracting the the Zones & Plots, Finding the difference and then dissolving using Union.Now need to program the same in R.

Tried below steps in R but would not able to get the right results , looking for ways how to dissolve this type of overlapping ploygons in R:

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

#Importing the shape files

field_boundary_fp <- "Database/gadenstedt2_outer_field 3 -26_0_zoned-plotm.shp"
poly_map_proj_str <- "+proj=longlat +datum=WGS84 +no_defs";
utm_target_proj   <- "+init=epsg:32632";

field_boundary_sdf <- maptools::readShapePoly(fn = field_boundary_fp,
                                          proj4string =  CRS(poly_map_proj_str),
                                          repair = T,
                                          delete_null_obj = T,
                                          verbose = T);
spplot(
field_boundary_sdf,"Rx"
)

# Extracting the Zones and Plots#

Zone_sdf <- field_boundary_sdf[field_boundary_sdf@data$Type == "Zone", ]
Plot_sdf <- field_boundary_sdf[field_boundary_sdf@data$Type == "Plot", ]
plot(Plot_sdf)
plot(Zone_sdf)

#Finding the Intersection Part between the both
test <- gIntersection(Zone_sdf, Plot_sdf,id="ZoneIdx")
plot(test)
plot(test, add = T, col = 'blue')

# Finding the difference

test2 <- gDifference(Zone_sdf,Plot_sdf,id="ZoneIdx")
plot(test2)
plot(test2, add = T, col = 'red')

#Trying for Union then
polygon3 <- gUnion(test2, Plot_sdf,id="ZoneIdx")
plot(polygon3)
plot(polygon3, add = T, col = 'yellow')
4

2 回答 2

1

Read the shapefile

library(raster)
fields <- shapefile("gadenstedt2_outer_field 3 -26_0_zoned-plotm.shp")

First separate the zones and fields. 

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

Erase plot from the zone

d <- erase(zone, plot)

Then append plot to d 

r <- bind(plot, d)

And now aggregate

rd <- aggregate(r, "Rx")
spplot(rd, "Rx")

---- Now with a reproducible example, so that others can also benefit; you should not ask questions that depend on a file that needs to be downloaded ----

Example data (two SpatialPolygonDataFrame objects)

library(raster)
p <- shapefile(system.file("external/lux.shp", package="raster"))
p <- aggregate(p, "NAME_1")
p$zone <- 10 + (1:length(p))
r <- raster(ncol=2, nrow=2, vals=1:4, ext=extent(6, 6.4, 49.75, 50), crs=crs(p))
names(r) <- "zone"
b <- as(r, 'SpatialPolygonsDataFrame')

erase and append

e <- erase(p, b)
pb <- bind(e, b)

data.frame(pb)
#        NAME_1 zone
#1     Diekirch   11
#2 Grevenmacher   12
#3   Luxembourg   13
#4         <NA>    1
#5         <NA>    2
#6         <NA>    3
#7         <NA>    4
于 2019-04-03T11:06:35.900 回答
0

To ensure Solution works on all Fields , added an extra line of code to above solution to add Buffer Geometry.

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)
spplot(d, "Rx")

r <- bind(plot, d)

rd <- aggregate(r, "Rx")

spplot(rd, "Rx")
于 2019-05-20T12:41:15.233 回答