I’m (very) new to R. I’m trying to produce a map of the world with the following features:
- ISEA3H grid (with a spacing of 6 miles between hexagons (I don’t have enough memory to generate one with 6 mile hexagons; the closest I’ve gotten is 24 mile hexagons)).
- No political borders (i.e., only “natural borders, like coastlines, lakes, rivers, etc.) I’ve got this down as well.
- Elevation/relief data. This is where I’m having trouble.
- The cell number printed inside the cells. I haven’t gotten this far yet.
Specifically, the problem I’m having is I can’t seem to “mesh” the data from my tif file with my map of the world.
I’m using the dggridR package, which produces an ISEA3H grid of the world with variable spacing. I’ve gotten pretty far.
Here’s my code thus far:
library(rgdal) library(sp) library(dggridR) library(mapdata) library(rnaturalearth) library(rnaturalearthdata) library(raster) # I've increased the spacing to 500 miles so the map doesn't take so long to render between iterations dggs <- dgconstruct(projection = "ISEA", aperture = 3, topology = "HEXAGON", precision = 7, spacing = 500, metric = FALSE) # Save a fairly large image at the end of it jpeg("mymap.jpg", width = 9000, height = 6000, quality = 95) # I'm using the physical data from https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/10m_physical.zip # and the raster data from https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/raster/NE2_HR_LC_SR_W.zip world <- readOGR("/home/mario/Downloads/data/10m_physical", "ne_10m_land") rivers <- readOGR("/home/mario/Downloads/data/10m_physical", "ne_10m_rivers_lake_centerlines") ocean <- readOGR("/home/mario/Downloads/data/10m_physical", "ne_10m_ocean") lakes <- readOGR("/home/mario/Downloads/data/10m_physical", "ne_10m_lakes") coast <- readOGR("/home/mario/Downloads/data/10m_physical", "ne_10m_coastline") rel <- raster("/home/mario/Downloads/data/NE2_HR_LC_SR_W/NE2_HR_LC_SR_W.tif") rel_spdf <- as(rel, "SpatialPixelsDataFrame") rel <- as.data.frame(rel_spdf) # use grid of the entire earth sa_grid <- dgearthgrid(dggs,frame=TRUE, wrapcells=TRUE) # actually plot the data p<- ggplot() + # what should be the relief layer geom_tile(data = rel, aes(x = "x", y = "y")) + # the world geom_polygon(data=world, aes(x=long, y=lat, group=group), fill=NA, color="black") + # lakes geom_polygon(data=lakes, aes(x = long, y = lat, group = group), fill = '#ADD8E6') + # more lakes geom_path(data=lakes, aes(x = long, y = lat, group = group), color = 'blue') + # rivers geom_path(data=rivers, aes(x=long, y=lat, group=group), alpha=0.4, color="blue") + # coastline geom_path(data=coast, aes(x = long, y = lat, group = group), color = 'blue') + # hexagonal grid geom_polygon(data=sa_grid, aes(x=long, y=lat, group=group), fill="white", alpha=0.4) + # hexagonal grid geom_path(data=sa_grid, aes(x=long, y=lat, group=group), alpha=0.4, color="grey") + # some necessary r code that I don't understand coord_equal() # change from flat map to globe p+coord_map("ortho", orientation = c(41.0082, 28.9784, 0))+ xlab('')+ylab('')+ theme(axis.ticks.x=element_blank())+ theme(axis.ticks.y=element_blank())+ theme(axis.text.x=element_blank())+ theme(axis.text.y=element_blank())+ ggtitle('World map') # finish writing to jpeg dev.off()
This is as far as I’ve gotten: https://imgur.com/y9LPqVS
The code above is currently chugging away, and has been for the past 4 hours. It’s remained within the bounds of the memory of my machine, so that’s a good sign.
What is the idiomatic way of projecting a geodesic grid over a relief map of the Earth? How could I include the relief data from the naturalearth tif file with my code thus far?