This notebook creates the hexagon database using R's sf
package.
The workflow:
sf
package.#Import packages
library(sf)
library(dplyr)
Warning message: "package 'sf' was built under R version 3.6.3" Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1 Warning message: "package 'dplyr' was built under R version 3.6.3" Attaching package: 'dplyr' The following objects are masked from 'package:stats': filter, lag The following objects are masked from 'package:base': intersect, setdiff, setequal, union
#set input & output filenames
fname_need <- '../Data/raw/Basic_Data/1_Conservation_Need.shp'
fname_opps <- '../Data/raw/Basic_Data/5_Combined_Opportunity.shp'
fname_hex <- '../Data/raw/Hexagon_files/100km2_D1_LamAzimNeed.shp'
fname_output <- '../scratch/Hexagons_from_R.shp'
#Read in the needs shapefile
fcNeed <- st_read(fname_need,quiet=TRUE)
#Read in the opportunities shapefile
fcOpps <- st_read(fname_opps,quiet=TRUE)
#Read in the hexagons shapefile
fcHex <- st_read(fname_hex,quiet=TRUE)
#Transform the needs feature class to match the coordinate system of the hexagons
## *Note*: We may be introducing some error here by
## projecting the hexagons to the needs polygon dataset.
fcHex <- st_transform(fcHex,st_crs(fcNeed))
#Create a geoseries of centroids from the hexagon polygon geoseries
hex.centroids <- st_centroid(fcHex['geometry'])
#Iterate through each record in the NEEDs feature class
for (row in 1:nrow(fcNeed)) {
#Set the output column name from the row's FID value
col.name <- paste("need",row,sep='')
#Get the shape
sel.shape <- fcNeed[row,'geometry']
#Construct a Boolean mask where the hexagon features intersect the shape
sel.mask = st_contains(sel.shape,hex.centroids,sparse=FALSE)
#Set default value in a new column to zero
fcHex[col.name] = 0
#Update masked features to a value of 1
fcHex[sel.mask,col.name] = 1
}
#Iterate through each record in the OPPORTUNITYs feature class
for (row in 1:nrow(fcOpps)) {
#Set the output column name from the row's FID value
col.name <- paste("opps",row,sep='')
#Get the shape
sel.shape <- fcNeed[row,'geometry']
#Construct a Boolean mask where the hexagon features intersect the shape
sel.mask = st_contains(sel.shape,hex.centroids,sparse=FALSE)
#Set default value in a new column to zero
fcHex[col.name] = 0
#Update masked features to a value of 1
fcHex[sel.mask,col.name] = 1
}
#Create a total needs column
fcHex$NeedTotal <- as.data.frame(fcHex) %>%
select(starts_with("need")) %>%
rowSums()
#Create a total opportunites column
fcHex$OppsTotal <- as.data.frame(fcHex) %>%
select(starts_with("opps")) %>%
rowSums()
#Write results
st_write(fcHex,fname_output,delete_layer=TRUE)
Deleting layer `Hexagons_from_R' using driver `ESRI Shapefile' Writing layer `Hexagons_from_R' to data source `../scratch/Hexagons_from_R.shp' using driver `ESRI Shapefile' Writing 24918 features with 116 fields and geometry type Polygon.