4

I am trying to produce an outline for a hexagonal cartogram by dissolving the inner polygons via the unionSpatialPolygons or aggregate functions. I am getting stray hexs that do not dissolve... a dummy example to show the problem:

# grab a dummy example shape file
library(raster)
g <- getData(name = "GADM", country = "GBR", level = 2)
# par(mar = rep(0,4))
# plot(g)

# create a hexagonal cartogram
# library(devtools)
# install_github("sassalley/hexmapr")    
library(hexmapr)
h <- calculate_cell_size(shape = g, seed = 1,
                         shape_details = get_shape_details(g), 
                         learning_rate = 0.03, grid_type = 'hexagonal')
i <- assign_polygons(shape = g, new_polygons = h)
par(mar = rep(0,4))
plot(i)

enter image description here

# dissolve the polygons to get coastline
library(maptools)
j <- unionSpatialPolygons(SpP = i, IDs = rep(1, length(i)))
par(mar = rep(0,4))
plot(j)

# same result with aggregate in the raster package
k <- aggregate(x = i)
par(mar = rep(0,4))
plot(k)

enter image description here

With the shapefile I am actually using (not for the UK) I get even more stray hexagons - some complete - some not.

guyabel
  • 8,014
  • 6
  • 57
  • 86
  • The two stray polygons coincide with the inland extremities of two ocean inlets, the northern one being the Irish Sea, the southern being the Bristol Channel. So it might be some glitch associated with these areas of coastline. I'm not familiar with `hexmapr` - is it possible to reduce the size of the hexagons and compare results? – Stuart Allen Dec 06 '17 at 22:02
  • @StuartAllen i think you are right. both are holes in `j`. – guyabel Dec 07 '17 at 05:56

1 Answers1

2

Suggested solution from Roger Bivand (via an email exchange):

 g1 <- spTransform(x = g, CRSobj = CRS("+init=epsg:27700"))
 # cellsize from calculate_cell_size() above
 h1 <- spsample(x = g1, type="hexagonal", cellsize=38309) 
 i2 <- HexPoints2SpatialPolygons(hex = h1) 
 j2 <- unionSpatialPolygons(SpP = i2, IDs = rep(1, length(i2)))
 plot(j2)

i.e. avoid assign_polygons() in hexmapr and utilize 1) spsample to generate shape positions and 2) HexPoints2SpatialPolygons for the hexagonal grid (both in sp package).

enter image description here

guyabel
  • 8,014
  • 6
  • 57
  • 86