Chapter 7 Mapping Exercise

Now let's conduct an example where we download GIS data from a repository and map it. So let's download two geospatial datasets and create a map with those two layers. First, visit http://www.gis.pr.gov/descargaGeodatos/Pages/default.aspx and click on "Conservación", as shown below.

Next, click on "Critical wildlife areas" to download a shapefile of the critical wildlife areas in Puerto Rico. Save the data into a folder named "wildlife" and import the data using the following functions:

# Import our libraries to import the data (rgdal) and then map the data (leaflet) 
library(rgdal) # contains functions for importing and wrangling vector data
library(leaflet) # contains functions for mapping our data 

wildlife_dir <- "~/Desktop/wildlife/critical_wildlife_areas.shp"
wildlife <- readOGR(wildlife_dir)
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/eminefidan/Desktop/wildlife/critical_wildlife_areas.shp", layer: "critical_wildlife_areas"
## with 225 features
## It has 13 fields

Let's go ahead and download another dataset for mapping by visiting this GIS data repository. Next, find and click on the Download button and select "Shapefile" to download the dataset. Go ahead and import this dataset into R:

hospitals_dir <- "~/Downloads/Hospitals-shp/Hospitals.shp"
hospitals <- readOGR(hospitals_dir)
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/eminefidan/Downloads/Hospitals-shp/Hospitals.shp", layer: "Hospitals"
## with 7596 features
## It has 32 fields

Now that the data is imported into R, we should check the coordinate system and project it to a geographic coordinate system that leaflet prefers. The code below demonstrates this:

crs(wildlife)
## CRS arguments:
##  +proj=lcc +lat_0=17.8333333333333 +lon_0=-66.4333333333333
## +lat_1=18.4333333333333 +lat_2=18.0333333333333 +x_0=200000 +y_0=200000
## +ellps=GRS80 +units=m +no_defs
crs(hospitals)
## CRS arguments:
##  +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1
## +units=m +nadgrids=@null +wktext +no_defs
new_crs <- '+init=epsg:4326' 
wildlife_projected <- spTransform(wildlife, CRS(new_crs))
hospitals_projected <- spTransform(hospitals, CRS(new_crs))

# Double check the new projection of our data
crs(wildlife_projected)
## CRS arguments: +proj=longlat +datum=WGS84 +no_defs
crs(hospitals_projected)
## CRS arguments: +proj=longlat +datum=WGS84 +no_defs

If we aren't sure what type of vector data our GIS datasets are, then we can we the plot() function to generate a basic plot of our data. Within the plot function we provide the variable name of our dataset.

plot(wildlife_projected)

plot(hospitals_projected)

Although these plots aren't very pretty, it gives us some good information regarding the data type. From this figure, it is apparent that the critical wildlife areas dataset is a polygon vector data type and the hospitals dataset is a point vector data type. Thus, our next step is to create our leaflet map using the addPolygons() and addCircleMarkers() functions.

leaflet() %>% addTiles() %>% 
  addPolygons(data = wildlife_projected) %>% 
  addCircleMarkers(data = hospitals_projected)