[GIS] How to overlay point co-ordinates on tmap shapefile

ggmapr

So far, I've managed to visually present data for different socio-economic factors at a district level in Pakistan with the help of the tmap package in RStudio to create the following output.

At the same time, I have used the ggmap package to plot the geographic co-ordinates of various towers for different Mobile Network Operators in Pakistan (roughly 13000 towers). The output is as follows.

My next requirement is to overlay the second plot over the first, essentially showing the towers on top of the district wise electricity scores.

Any ideas on which packages I could use (as well as sample code) in order to achieve this?

Best Answer

This can also be done with tmap. Have a look at the tmap vignette, especially the Plotting with tmap elements section.

In your case your code ought to look basically like this:

tm_shape(shapefile_with_electricity_data) + 
  tm_fill("Electricity_in_Schools_variable") + 
tm_shape(shapefile_with_MNO_data) + 
  tm_dots("MNO_variable")

Related Solutions

[GIS] ggmap: plot polygon from shapefile

as.data.frame() does not work for SpatialPolgons in geom_polygon, because the geometry gets lost. You have to use ggplot2::fortify (may be deprecated in the future, see ?fortify). The recommended way is now to use broom::tidy:

R> library("broom")
R> head(tidy(kommune))
Regions defined for each Polygons
   long   lat order  hole piece group  id
1 10.29 59.72     1 FALSE     1 153.1 153
2 10.32 59.70     2 FALSE     1 153.1 153
3 10.32 59.69     3 FALSE     1 153.1 153
4 10.31 59.68     4 FALSE     1 153.1 153
5 10.30 59.67     5 FALSE     1 153.1 153
6 10.28 59.67     6 FALSE     1 153.1 153

But another problem arises with your example. Since the polygon is larger than the map extent, ggmap does not correctly clip the polygon. ggmap sets the limits on the scale, this will throw away all data that's not inside these limits.

Here is a modified version of your code:

p <- ggmap(map, extent = "normal", maprange = FALSE) +
     geom_point(data = as.data.frame(subscr),
                aes(x = lon, y = lat, size=pop),
                colour = "darkgreen") +
     geom_polygon(data = fortify(kommune),
                  aes(long, lat, group = group),
                  fill = "orange", colour = "red", alpha = 0.2) +
     theme_bw() +
     coord_map(projection="mercator",
               xlim=c(attr(map, "bb")$ll.lon, attr(map, "bb")$ur.lon),
               ylim=c(attr(map, "bb")$ll.lat, attr(map, "bb")$ur.lat))

print(p)

R Tmap Package – How to Limit X and Y Axes

Have you read the help for tm_shape? It gives the following parameter:

bbox: bounding box. One of the following:

        • A bounding box (an ‘sf’ bbox object, see ‘st_bbox’, a 2
          by 2 matrix (used by the ‘sp’ package), or an ‘Extent’
          object used by the ‘raster’ package).

Its in map units rather than lat-long, so with the World data from the help you can do:

tm_shape(World, bbox=tmaptools::bb(matrix(c(0,0,2000000,2000000),2,2))) +tm_polygons()

Best Answer

Related Solutions

[GIS] ggmap: plot polygon from shapefile

R Tmap Package – How to Limit X and Y Axes

Related Question