Re: "Plotting a map with ggplot2, color by tile"

[duleise]

I've been trying to use shapefiles and plot maps for quite a while. This is the first time I really rolled my sleeves up. Your codes worked smoothly. Thank you for sharing.

I have a question about some of the details you used in your workflow. Hope you could give some advices.

Excluding inland waters You removed lakes by excluding their names, which is convenient only if one knows and the data provide the lake names. Besides, there are also rivers (canals in NL?) in land. Fortunately there are shapefiles for waters that I can find online, which, like the land maps, provide polygon info for inland waters. If I cover the waters polygons in front of land polygons, with different color, it looks nice. Thing is, the waters polygons do not include administrative information, which is understandable since a river can cross many places. But when I try to plot a map for the land and waters of an individual province, I can not exclude waters from other provinces. This makes the whole picture much larger than the specific province which I intend to plot. Do you have any suggestion on this?

Thank you.

[RMHogervorst]

Hi @duleise (Lei), I'm just beginning with this whole shapefile thing. Yeah I was very fortunate that those waters were specified.

But I believe you can modify the objects so that rivers are excluded from the main chart. excluding the information from the land/province info. I don't actually now how to do it, though. I guess the packages sp and maptools are supposed to work with this magic. The vignette of maptools https://cran.r-project.org/web/packages/maptools/vignettes/combine_maptools.pdf seems to talk about merging files. I think the difficulty is in making sure the shapefiles are in the same projection and then finding the right way to modify the shapefiles.

• load both of the files.
• 'substract' the rivers from the province file
• possibly fixing the polygons because the polygons might not be closed anymore

Let me know what you found and how you fixed it, because I would also like to do something like this!

[RMHogervorst]

I haven't looked or tried out the answers here but these are some links I found while searching for 'substract shapefile from shapefile r'

• http://gis.stackexchange.com/questions/140504/extracting-intersection-areas-in-r
• http://gis.stackexchange.com/questions/133264/subtracting-erasing-polygons-from-polygon-using-arcgis-for-desktop
• https://stat.ethz.ch/pipermail/r-sig-geo/2005-December/000696.html
• http://gis.stackexchange.com/questions/155328/merging-multiple-spatialpolygondataframes-into-1-spdf-in-r

[duleise]

Thank you for your suggestion. I have a quick and very-likely-dirty solution from raster. Let me take Iceland as an example.

library(sp)
library(maptools)
library(dplyr)
library(ggplot2)
library(raster)

IS <- readRDS("ISL_adm2.rds")
IS.water <- readShapePoly("ISL_wat/ISL_water_areas_dcw.shp")

# say i want to plot Austurland and exclude waters in this area
x <- raster::intersect(IS.water,
                       subset(IS, NAME_1 == "Austurland"))

# raster::crop and rgeos::gDifference may have same/similar function.

p <- ggplot() + 
  theme_minimal() + 
  xlab(NULL) + ylab(NULL) +
  theme(panel.grid.minor=element_blank(),
        panel.grid.major=element_blank(),
        axis.ticks = element_blank(),
        axis.text = element_blank(),
        legend.position="none") +
  geom_polygon(data = subset(IS, NAME_1 == "Austurland"),
               aes(x = long, y = lat, 
                   group = group),
               fill = "forestgreen") +
  geom_polygon(data = x,
               aes(x = long, y = lat, 
                   group = group),
               fill = "firebrick") +
  coord_map()
p

This is definitely not tidy or neat. I will dig more for a better solution, before which I think I need a lot of background knowledge about GIS. Thank you again.

[RMHogervorst]

Very interesting! I will try some of this as well.

[duleise]

https://github.com/Robinlovelace/Creating-maps-in-R/blob/master/intro-spatial.Rmd After reading this very helpful tutorial, which was cited by Martijn, the author of tmap, I think I found the silver bullet for my question. sp::over and rgeos::gIntersects The former can be called just by using square brackets. So in my quick and dirty codes, x can be made by x <- IS.water[subset(subset(IS, NAME_1 == "Austurland")), ] but before doing this, the CRS of IS and IS.water need to be consistent. proj4string(IS.water) <- proj4string(IS)

Thank you again for a bright inspiration. Happy mapping.