library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(raster)
## Loading required package: sp
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
##
## intersect, select, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(spData)
## Warning: package 'spData' was built under R version 4.0.3
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
library(tmap)
library(leaflet)
library(ggplot2)
library(maptools)
## Warning: package 'maptools' was built under R version 4.0.3
## Checking rgeos availability: TRUE
#reading data
ar_cnty <- sf::st_read(dsn = "data/county_nrcs_a_ar.gdb")
## Reading layer `county_nrcs_a_ar' from data source `D:\Fall 2020\Seminar\Assignments\Assignment 10\Tmap\data\county_nrcs_a_ar.gdb' using driver `OpenFileGDB'
## Simple feature collection with 75 features and 31 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 165469.4 ymin: -147183.4 xmax: 612641.6 ymax: 243444.9
## projected CRS: NAD83 / Arkansas North
dtgw <- read_sf("data/semAs10.shp")
st_crs(dtgw)
st_crs(ar_cnty)
#Changing the CRS to be the same
dtgw_prj <- st_transform(dtgw, st_crs(ar_cnty))
#checking CRS
st_crs(dtgw_prj)==st_crs(ar_cnty)
## [1] TRUE
#Clip geometry
ar_gw_intsct <- st_intersection(dtgw_prj, ar_cnty)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
plot(ar_gw_intsct)
## Warning: plotting the first 10 out of 34 attributes; use max.plot = 34 to plot
## all
ar_gw_intsct
## Simple feature collection with 151 features and 34 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 394548.7 ymin: -88770 xmax: 570389.3 ymax: 239138.2
## projected CRS: NAD83 / Arkansas North
## # A tibble: 151 x 35
## Longitude Latitude X2016dtgw_m OBJECTID FIPS_C FIPS_I FIPSST FIPSCO STPO
## * <dbl> <dbl> <dbl> <int> <chr> <dbl> <chr> <chr> <chr>
## 1 -91.0 34.3 1.59 94 05107 5107 05 107 AR
## 2 -91.0 34.3 3.9 94 05107 5107 05 107 AR
## 3 -90.7 34.5 2.74 94 05107 5107 05 107 AR
## 4 -90.9 34.5 6.9 94 05107 5107 05 107 AR
## 5 -91.0 34.4 4.8 94 05107 5107 05 107 AR
## 6 -90.8 34.5 6.3 94 05107 5107 05 107 AR
## 7 -90.9 34.5 7.5 94 05107 5107 05 107 AR
## 8 -91.0 34.5 5.24 94 05107 5107 05 107 AR
## 9 -90.7 34.5 2.57 94 05107 5107 05 107 AR
## 10 -91.0 34.6 8.4 94 05107 5107 05 107 AR
## # ... with 141 more rows, and 26 more variables: COUNTYNAME <chr>,
## # CNTYDISP <chr>, CNTYSHORT <chr>, CNTYSORT <chr>, CNTYCATEGO <chr>,
## # CNTYACTIVE <chr>, INDEPCITY <chr>, CNTYSTAND <chr>, SEATLAT <dbl>,
## # SEATLONG <dbl>, NAD83UTM <dbl>, NAD83STATE <dbl>, NAD27STATE <dbl>,
## # STATENAME <chr>, CNTYSTARTD <dttm>, CNTYENDD <dttm>, LASTCHGD <dttm>,
## # NOTE <chr>, BOTTOM <dbl>, TOP_ <dbl>, LEFT_ <dbl>, RIGHT_ <dbl>,
## # GlobalID <chr>, Shape_Length <dbl>, Shape_Area <dbl>, geometry <POINT [m]>
#add raster object for visualization
gw_depth <- raster("data/2016 Depth to water1.tif")
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum Unknown based on GRS80 ellipsoid in CRS definition
crs(gw_depth)
## CRS arguments:
## +proj=lcc +lat_0=34.3333333333333 +lon_0=-92 +lat_1=36.2333333333333
## +lat_2=34.9333333333333 +x_0=400000 +y_0=0 +ellps=GRS80 +units=m
## +no_defs
## #tmap#
# Add fill layer to ar_cnty shape
tm_shape(ar_cnty) +
tm_fill()
# Add border layer to ar_cnty shape
tm_shape(ar_cnty) +
tm_borders()
# Add fill and border layers to ar_cnty shape
tm_shape(ar_cnty) +
tm_fill() +
tm_borders()
#Map objects
map_cnty = tm_shape(ar_cnty) + tm_polygons()
class(map_cnty)
## [1] "tmap"
print(map_cnty)
map_dtgw = tm_shape(ar_gw_intsct)+tm_dots()
print(map_dtgw)
#Adding new shapes
map_ar_gw = map_cnty +
tm_shape(ar_gw_intsct) + tm_dots(alpha = 0.5)
print(map_ar_gw)
tmap_arrange(map_cnty, map_dtgw, map_ar_gw)
map_gw_depth = tm_shape(gw_depth)+tm_raster()
print(map_gw_depth)
map_ar_dtgw = map_cnty + map_gw_depth
print(map_ar_dtgw)
tmap_arrange(map_cnty, map_gw_depth, map_ar_dtgw)
#add aesthetics
ma1 = tm_shape(ar_cnty) + tm_fill(col = "black")
ma2 = tm_shape(ar_cnty) + tm_fill(col = "black")
ma3 = tm_shape(ar_cnty) + tm_borders(col = "snow4")
ma4 = tm_shape(ar_cnty) + tm_borders(lwd = 1)
ma5 = tm_shape(ar_cnty) + tm_borders(lty = 2)
ma6 = tm_shape(ar_cnty) + tm_fill(col = "black") +
tm_borders(col = "snow4", lwd = 1, lty = 2)
tmap_arrange(ma1, ma2, ma3, ma4, ma5, ma6)
plot(st_geometry(ar_cnty), col = ar_gw_intsct$X2016dtgw_m)
legend_title = expression("Depth to Water (m)")
ma7 <- ma6 + tm_shape(ar_gw_intsct, position = c("right", "center")) +
tm_dots(col = "X2016dtgw_m", title = legend_title, n = 7,
palette = "Reds", style = "equal", size = 2) +
tm_legend(legend.text.col = "black", legend.bg.color = "white",
legend.title.col = "black", position = c("LEFT", "top")) +
tm_compass(type = "4star", position = c("LEFT", "BOTTOM"), size = 0.7) +
tm_scale_bar(breaks = c(0, 50), text.size = 0.5,
position = c("RIGHT", "BOTTOM")) + tm_layout(bg.color = "lightblue",
title.color = "darkslategray1", legend.outside = TRUE,
legend.title.size = 0.9, legend.frame = "black", legend.text.size = 0.9)
print(ma7)
## #Interactive Map#
tmap_mode("view")
## tmap mode set to interactive viewing
print(ma7)
## Compass not supported in view mode.
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.