GIS and Mapping with R

Distance calculations

The goal is combining national park and city databases. They are both in lat-long projections, so no transformation needed.

library(rgdal)
library(sp)
## Assuming the GIS files are stored in the working directory (so first argument is the blank string '')
USparks=readOGR('','ne_10m_parks_and_protected_lands_area')  ## or readShapePoly
UScities=readOGR('','citiesx010g')                           ## or readShapePoints
Worldparks=readOGR('','INTpol_RWB')
Worldcities=readOGR('','ne_10m_populated_places')

Extract a single park using subset. Find how many Polygons it has.

isleroyale=(subset(USparks,unit_name=='Isle Royale NP'))
length(isleroyale@polygons)                      ## This should be 1 because one park was extracted

[1] 1

npoly=length(isleroyale@polygons[[1]]@Polygons)  ## For a park with a complex map, this might be > 1
npoly

[1] 1

Subset cities within a set distance of the bounding box of the park (perhaps 1 degree latitude in all directions).

closecity=subset(UScities,LONGITUDE<isleroyale@bbox['x','max']+1)
closecity=subset(closecity,LONGITUDE>isleroyale@bbox['x','min']-1)
closecity=subset(closecity,LATITUDE<isleroyale@bbox['y','max']+1)
closecity=subset(closecity,LATITUDE>isleroyale@bbox['y','min']-1)
dim(closecity@coords)

[1] 26  2

# sp::plot(closecity)          ## Map omitted
# sp::plot(dehoop,add=TRUE)

Extract coordinates for the first polygon of the park.

coord1=isleroyale@polygons[[1]]@Polygons[[1]]@coords

Find distance from polygon to cities.

alldist=sqrt(dsquare(pts=coord1,src=coordinates(closecity)))