Hengistbury
library(raster)
library(mapview)
library(tmap)
library(giscourse)
library(sf)
library(tidyverse)
conn<-connect()
mapviewOptions(basemaps =c("Esri.WorldImagery","OpenStreetMap","Esri.WorldStreetMap", "CartoDB.Positron", "OpenTopoMap", "OpenStreetMap.BlackAndWhite", "OpenStreetMap.HOT","Stamen.TonerLite", "CartoDB.DarkMatter", "Esri.WorldShadedRelief", "Stamen.Terrain"))
query<-"select geom,lnr_name from local_reserves where lnr_name like 'Hengist%'"
lnr<-st_read(conn,query=query)
bound<-as(lnr, "Spatial")
mapview(lnr, alpha.regions = 0, legend=FALSE, lwd =3) %>% extras()
hh_habitat<-sprintf("select main_habit mhabitat,
s1habtype,s2habtype,s3habtype habitat,
st_intersection(h.geom,s.geom) geom from
(%s) s,
ph_v2_1 h
where st_intersects(h.geom,s.geom)",query)
habitats<-st_read(conn, query=hh_habitat)
habitats$area <- st_area(habitats)
habitats%>% st_buffer(0.001)%>% group_by(mhabitat) %>% summarise(area=sum(area))-> habitats
qtm(habitats, "mhabitat") +tm_style("classic")
mapview(habitats, zcol="mhabitat", burst=TRUE,legend=FALSE, hide = TRUE) %>% extras()
dsm<- pgGetRast(conn, "dsm2m",boundary=bound)
dtm<- pgGetRast(conn, "dtm2m",boundary=bound)
dsm[dsm<0.1]<-NA
dtm[dtm<0.1]<-NA
slope<-terrain(dtm, opt='slope', unit='degrees')
aspect<-terrain(dtm, opt='aspect',unit='degrees')
sloper<-terrain(dtm, opt='slope', unit='radians')
aspectr<-terrain(dtm, opt='aspect',unit='radians')
hillshade<-hillShade(sloper,aspectr,angle=5)
plot(hillshade, col = gray(0:100 / 100), legend = FALSE)
plot(dtm, col = terrain.colors(25), alpha = 0.3, legend = TRUE, add = TRUE)
