Plot points on OpenStreetMap w/ option to place smaller inset map in NW corner of main map

<pre>
#plotting X,Y coordinate points on basemap downloaded from R OpenSteetMap program
#code is also included to get a small inset map plotted into the NW corner of your main map


#set resolution of document to be saved
ppi=1200
pdf("your_file_name_here.pdf")

#load necessary libraries
library(rJava)
library(rgdal)
library(sp)
library(OpenStreetMap)

#grab basemap, you will probably need to adjust number of tiles to your refered level of zoom
petrolina_map=openmap(c(-9.19,-40.61),c(-9.5,-40.24),type="esri-topo", minNumTiles=13)

#read in capture points
coords=read.csv("file_with_XY_points_to_plot_on_map.csv")

#i have columns called lat for latitude and long for longtitude in my .csv file with my X,Y points to plot
attach(coords)

#bind the lat and long columns together, but put long as the first column because the reverse way does not plot the points correctly, but I am not sure why 
lat_long=cbind(long, lat)

#convert lat_long into spatial data that is latitude and longtide CRS WGS84
xya=SpatialPoints(lat_long)
proj4string(xya)=CRS("+proj=longlat +ellps=WGS84")

#the points have to be converted from longlat to merc to match CRS of the OpenStreetmap basemap
xyb=spTransform(xya,CRS("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0
+lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"))

plot(petrolina_map, removeMargin=TRUE, raster=TRUE)
plot(xyb, add=TRUE, pch=19)

#add in text to name your points on the map
xy=coordinates(xyb)
names.xy=c("1","2","3","4","5","6","7","8","9")
text(xy[1,1],xy[1,2], names.xy[1], adj=1, pos=1, lwd=1, ce=1.2)
text(xy[2,1],xy[2,2], names.xy[2], adj=1, pos=1, lwd=1, ce=1.2)
text(xy[3,1],xy[3,2], names.xy[3], adj=1, pos=1, lwd=1, ce=1.2)
text(xy[4,1],xy[4,2], names.xy[4], adj=1, pos=3, lwd=1, ce=1.2)
text(xy[5,1],xy[5,2], names.xy[5], adj=1, pos=3, lwd=1, ce=1.2)
text(xy[6,1],xy[6,2], names.xy[6], adj=1, pos=3, lwd=1, ce=1.2)
text(xy[7,1],xy[7,2], names.xy[7], adj=1, pos=3, lwd=1, ce=1.2)
text(xy[8,1],xy[8,2], names.xy[8], adj=1, pos=3, lwd=1, ce=1.2)
text(xy[9,1],xy[9,2], names.xy[9], adj=1, pos=3, lwd=1, ce=1.2)

#code to get small inset map of Brazil
library(maps)
library(ggplot2)
library(ggmap)
library(grid)

#code to get inset map into larger map
brazilmap <- data.frame(map("world", "Brazil", plot = FALSE)[c("x", "y")])
# This should get your corner points for the box
insetrect <- data.frame(xmin = min(-40.5), xmax = max(-40.48),
ymin = min(-9.39), ymax = max(-9.3))

# The inset map, all of Brazil
a <- ggplot(brazilmap) +theme_bw(base_size = 22) +geom_path(data = brazilmap, aes(x, y),colour = "black", fill = "white") +theme(panel.border = element_rect(colour = "black",size = 1, linetype=1),panel.grid.major = element_blank(), panel.grid.minor=element_blank(),panel.background = element_rect( fill = "white"),legend.position = c(0.15,0.80), legend.key =element_blank(),axis.ticks = element_blank(), axis.text.x=element_blank(),axis.text.y=element_blank()) + labs(x = "", y ="")+geom_rect(data = insetrect, aes(xmin = xmin, xmax = xmax,ymin = ymin, ymax = ymax), alpha=0, colour="blue", size = 2, linetype=1)

#print the two maps together
vpa_ <- viewport(width = 0.37, height = 0.37, x = 0.15, y = 0.76) # for inset mapprint(a, vp = vpa_)
print(a, vp = vpa_)
dev.off()

Leave a comment