0

Hopefully this is a simple query. I have pieced together some code to make a topography plot using elevation data and ggplot2.

However, now I would like to overlay onto this plot the country outlines, along with points to show the location / elevation of points of interest e.g. cities and also other polygons and lines to show for example routes. However I am struggling to see how best to do this. Any help or suggestions would be appreciated.

My code is as follows:

library(dplyr)
library(ggplot2)
library(raster)
library(scales)
library(rgeos)

#get elevation data
dem1<- getData("SRTM",lat=34.5,lon=33)
dem2<- getData("SRTM",lat=34.5,lon=36)
dem <- merge(dem1,dem2)

#create polygon to crop the elevation data file   
library(sp)
coords = matrix(c(34, 31,
                  34, 32.5,
                  36, 32.5,
                  36, 31.0,
                  34, 31), 
                ncol = 2, byrow = TRUE)
P1 = Polygon(coords)
Ps1 = SpatialPolygons(list(Polygons(list(P1), ID = "a")), proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

#crop the elevation data using the polygon
dem <- crop(dem,Ps1)

#lower the reolution to enable faster plotting
dem_lower_res <- aggregate(dem, fact=10)
dem.p  <-  rasterToPoints(dem_lower_res )
df <-  data.frame(dem.p)
colnames(df) = c("lon", "lat", "alt")

#plot the data
ggplot(df, aes(lon,lat)) +
  geom_raster(aes(fill = alt))+
  scale_fill_gradientn(colours = terrain.colors(10))
1
  • Welcome to GIS.SE. May this that help you? Apr 28 '20 at 12:02
0

Just as an example, I used GADM country outline and two city locations (notional).

  • GADM dataset was downloaded as sp object, and plotted by geom_polygon()
  • A dataframe of two points (my_points) was converted to sf object, and plotted by geom_sf. Also geom_sf_label to print the label (city name) may be of interest for you.

It highly depends what kind of point/polygon/line you have... hope this helps.


library(raster)
#get elevation data
dem1<- getData("SRTM",lat=34.5,lon=33)
dem2<- getData("SRTM",lat=34.5,lon=36)
dem <- merge(dem1,dem2)

# get country border ################################################
gadm <- getData('GADM', country = "ISR", level = 0)
# use st_read() if you have a shapefile

# Add points  #######################################################
my_points = data.frame(
  city = c('Jerusalem', 'Tel Aviv-Yafo'),
  lon =  c(35.1, 34.9),
  lat =  c(31.8, 32.1)
)
library(sf)
my_points = st_as_sf(my_points, coords = c("lon", "lat"), crs = 4326)

# simplified your code a little, using extent() function ############
#create polygon to crop the elevation data file   
Ps1 = as(extent(34, 36, 31, 32.5), 'SpatialPolygons')
crs(Ps1) = "+proj=longlat +datum=WGS84 +no_defs"

#crop the elevation data using the polygon
dem = crop(dem, Ps1, snap= 'out')

#lower the reolution to enable faster plotting
dem_lower_res <- aggregate(dem, fact=10)
dem.p  <-  rasterToPoints(dem_lower_res )
df <-  data.frame(dem.p)
colnames(df) = c("lon", "lat", "alt")

#plot the data
library(ggplot2)
ggplot(df) +
  geom_raster(aes(lon, lat, fill = alt))+
  scale_fill_gradientn(colours = terrain.colors(10))+
  geom_sf(data= my_points, color= 'red')+                                         # point
  geom_sf_label(data= my_points, aes(label= city), nudge_y = 0.1,  color= 'red')+ # label
  geom_polygon(data= gadm, aes(long, lat), color= 'blue', fill= NA) +             # polygon
  coord_sf(xlim = c(34, 36), ylim = c(31, 33))+
  theme_bw()

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.