5

I wonder how it's possible to add this topographic information (elevation) to the map of Santa Cruz Island and all of the Galápagos Island that I have here.

library(rgdal)
library(rgeos)
library(ggplot2)
library(ggthemes)
library(ggsn)
myGPS = data.frame(lat=c( -0.6850556,-0.6854722,  -0.6857778  ),lon=c(-90.22275,-90.22261,  -90.22272)) 
coord.deg = myGPS

class(coord.deg)
## "data.frame"
coordinates(coord.deg)<-~lon+lat
class(coord.deg)
## "SpatialPointsDataFrame"
## attr(,"package")
## "sp"

# does it have a projection/coordinate system assigned?
proj4string(coord.deg) # nope
## NA

# Manually tell R what the coordinate system is
proj4string(coord.deg)<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")

# now we can use the spTransform function to project. We will project
# the mapdata and for coordinate reference system (CRS) we will
# assign the projection from counties

coord.deg<-spTransform(coord.deg, CRS(proj4string(gal)))
# double check that they match
identical(proj4string(coord.deg),proj4string(gal))
## [1] TRUE
my_pts <- SpatialPointsDataFrame(coords = coord.deg, data=data.frame(id=1:length(coord.deg)))

URL <- "https://osm2.cartodb.com/api/v2/sql?filename=public.galapagos_islands&q=select+*+from+public.galapagos_islands&format=geojson&bounds=&api_key="
fil <- "gal.json"
if (!file.exists(fil)) download.file(URL, fil)
gal <- readOGR(fil, "OGRGeoJSON")
gal <- gSimplify(gUnaryUnion(spTransform(gal, CRS("+init=epsg:31983")), id=NULL), tol=0.001)
gal_map <- fortify(gal)
rand_pts <- spTransform(my_pts, CRS("+init=epsg:31983"))

# ggplot can't deal with a SpatialPointsDataFrame so we can convert back to a data.frame
my_pts <- data.frame(my_pts)
my_pts.final = my_pts[,2:3]
# we're not dealing with lat/long but with x/y
# this is not necessary but for clarity change variable names
names(my_pts.final)[names(my_pts.final)=="lat"]<-"y"
names(my_pts.final)[names(my_pts.final)=="lon"]<-"x"

gg <- ggplot()
gg <- gg + geom_map(map=gal_map, data=gal_map,
                    aes(x=long, y=lat, map_id=id),
                    color="black", fill="#FFFFFF", size=.5)
gg <- gg + coord_equal()
gg <- gg + ggsn:::scalebar(gal_map, dist=50, location="bottomleft", st.size=5)
gg <- gg + theme_map()
gg <- gg + geom_point(data=my_pts.final, aes(x=x, y=y), color="red")
gg

Do you have any idea? Do you have a package in R that does this?

I added this:

topo <- readOGR("public.geojson.json", layer ="OGRGeoJSON")
topo <- spTransform(topo, CRS("+init=epsg:31983"))
gg + geom_map(map = topo_map, data = topo_map, aes(map_id=id))

But it creates a map with a black shape in it... enter image description here

  • Have you tried zooming in? At this extent, anything is going to look like garbage, except for larger boundaries. – MaryBeth Nov 20 '15 at 20:02
4

Your topomap is a SpatialLinesDataframe. geom_map is used for polygons. I suggest you use geom_path as below. It connects observations in original order. (geom_line would order by x value, which you also don't want).

topo <- readOGR("public.geojson.json", layer ="OGRGeoJSON") 
topo <- spTransform(topo, CRS("+init=epsg:31983"))
class(topo)
#[1] "SpatialLinesDataFrame"
topo_map <- fortify(topo)

# add the attributes back in
library(plyr)
topo_map <- join(topo_map, topo@data, by="id")

names(topo_map)
#[1] "long"       "lat"        "order"      "piece"      "group"     
#[6] "id"         "elev"       "cartodb_id" "created_at" "updated_at"

ggplot() + 
  geom_path(data = topo_map, aes(long, lat, group=group, colour=elev), size=0.1) +
  scale_colour_gradient(low="green", high="brown") +
  coord_equal()

enter image description here

| improve this answer | |
  • Do you have a way to add colours corresponding to the elevation? geom_path(data = topo_map, aes(long, lat, group=group),size=0.1,colour=elevation) ? – M. Beausoleil Nov 20 '15 at 23:21
  • 1
    fortify removes the attributes, so you have to join them back with the dataframe. See my additions above. – cengel Nov 20 '15 at 23:49

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.