Broken package for geographical coordinates

[rsbivand]

An SO user contacted me with this issue:

library(raster)
g <- getData(name = "GADM", country = "GBR", level = 2)
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)
plot(i)
library(maptools)
j <- unionSpatialPolygons(SpP = i, IDs = rep(1, length(i)))
plot(j)

I replied that your assign_polygons() function forces spDists to use planar coordinates (lots of warnings), and that polygon boundaries do not seem to be correctly assigned. It also takes much longer than:

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

which does not have the unconnected polygons.

[sassalley]

@rsbivand thanks for letting me know. I had no idea that it might be used for such applications. Given that the user does not seem to need assign_polygons (which was the main reason for the package) it seems like your response is definitely the most valuable approach.

Out of interest, would you suggest enforcing non-planar coordinates (i.e. using whatever spatial projection is native to the source data) as opposed to stripping the proj4string?

[rsbivand]

Very possibly - at least trap and reinforce/block the warnings from spDists. There may be an overlap with tilegramsR, perhaps with hexmapr generating tilegrams in the format @bhaskarvk uses, so comparing notes with regard to grids on geographical coordinates might make sense. See also dggridR (here too) for a good, out-of-the-box solution for global/continental representations on the sphere (@r-barnes - getting dggridR to use sf classes as well as current output would also be nice). There is a point of contact too with a recent sf issue wrt. maximum area overlay between the created grids and underlying geometries.

[guyabel]

@sassalley I was the person who posted the SO question. For what it is worth, I do not think there is a problem with unconnected polygons. As someone commented after the SO question the stray polygons are related to their geographic positions near an estuary and bay... so it was me being dumb. I confirmed this by adding some random data (I did use assign_polygons btw)...

library(raster)
g <- getData(name = "GADM", country = "GBR", level = 2)

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)

library(tidyverse)
library(broom)
set.seed(1)
ggplot(data = tidy(i) %>% 
         tbl_df() %>%
         group_by(id) %>%
         mutate(x = rnorm(1)), 
       mapping = aes(x = long, y = lat, group = group)) +
  geom_polygon(mapping = aes(fill = x)) +
  coord_map() +
  theme_bw()

... the polygons of interest are also labelled as holes.

Perhaps I should add - my motivation for the question was the desire to add some outlines to the plot of the admin level two hexagons as 1) I wanted to use a divergent colour scheme passing through white (the same colour as the background) 2) I wanted to represent admin level one regions. I did not get into this in the SO question as I thought the stray hexs were part of the same problem + I was trying to keep things simple

My full attempt was something like this...

j <- unionSpatialPolygons(SpP = i, IDs = i@data$NAME_1)
ggplot(data = tidy(i) %>% 
         tbl_df() %>%
         group_by(id) %>%
         mutate(x = rnorm(1)), 
       mapping = aes(x = long, y = lat, group = group)) +
  geom_polygon(mapping = aes(fill = x)) +
  coord_map() +
  theme_bw() +
  scale_fill_gradient2() +
  geom_polygon(data = tidy(j), fill = "transparent", colour = "black", size = 1.5)

with some unwanted pieces of lines through Scotland.

Apologies if there is any confusion with the above; I am well out my depth here, having never used R to do any serious mapping before.

[sassalley]

aha @gjabel thanks for letting me know and thanks for trying the package! Given the original context I thought you may have just wanted to create a boundary (and so assumed that you weren't using the assignments).

This may be a bit of a hack but a less well resolved boundary for the UK might help. Some more coarse borders do not separate out all of the islands on the west coast of Scotland. @rsbivand will know more about this but gSimplify discussed here might help with this.

Apologies if this is does not solve things (it would be great to know if the issue persists).

[rsbivand]

It isn't so likely that gSimplify or similar line generalisation will help, but it might be possible to filter out polygon parts of multipolygons that fall under an area threshold.

[r-barnes]

(@rsbivand I think getting sf to work with dggridR would be pretty cool, but my attention's pretty firmly elsewhere for the next few months.)

[rsbivand]

(@r-barnes OK, thanks for commenting, hope you'll be able to get back to this).