-
Notifications
You must be signed in to change notification settings - Fork 1
/
02_clean.R
199 lines (158 loc) · 9.29 KB
/
02_clean.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
# Copyright 2016 Province of British Columbia
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and limitations under the License.
# Spatial packages
library(rgdal)
library(sp)
library(rgeos)
library(raster)
library(bcmaps) # install using devtools::install_github("bcgov/bcmaps")
library(rmapshaper)
source("fun.R")
dir.create("tmp", showWarnings = FALSE)
## The CARTS database is downloadable from the Canadian Council on
## Ecological Areas here: http://www.ccea.org/carts/
carts <- readOGR("data/CARTS_Update_31122015.gdb", "CARTS_Update_31122015_WithoutQc", stringsAsFactors = FALSE)
## Read in the Conservation Lands database from: http://catalogue.data.gov.bc.ca/dataset/conservation-lands
cl_gdb <- "data/BCGW_conservation_lands_fgdb/WCL_CONSERVATION_LANDS_SP.gdb"
conservation_lands <- readOGR(cl_gdb, ogrListLayers(cl_gdb)[1], stringsAsFactors = FALSE)
rm(cl_gdb)
## Read in NGO conservation lands
fee_simple_ngo_lands <- readOGR("data", "BC_NGO_ConsDB_FeeSimple_31Dec2014_updated02Nov2015", stringsAsFactors = FALSE)
## Extract the year of SecDate (securement date) and store in prot_date
fee_simple_ngo_lands$prot_date <- as.numeric(substr(fee_simple_ngo_lands$SecDate1, 1, 4))
## Extract just BC from CARTS database
bc_carts <- carts[carts$LOC_E %in% c("British Columbia", "Offshore Pacific Marine"), ]
rm(carts)
## Subset only the Administrated Lands - Acquisitions and Leases
bc_admin_lands <- conservation_lands[conservation_lands$CONSERVATION_LAND_TYPE == "Administered Lands" &
(conservation_lands$TENURE_TYPE == "Acquisition" |
conservation_lands$TENURE_TYPE == "Lease" |
conservation_lands$TENURE_TYPE == "Transfer of Administration/Control"), ]
rm(conservation_lands)
## Transform CRS of all layers to BC Albers
bc_carts <- transform_bc_albers(bc_carts)
bc_admin_lands <- transform_bc_albers(bc_admin_lands)
fee_simple_ngo_lands <- transform_bc_albers(fee_simple_ngo_lands)
#### Simplify for testing
# bc_carts <- ms_simplify(bc_carts, 0.01, keep_shapes = TRUE)
# ecoregions <- ms_simplify(ecoregions, 0.01, keep_shapes = TRUE)
####
## Check validity of polygons of bc_carts, and fix:
bc_carts <- fix_geo_problems(bc_carts)
bc_admin_lands <- fix_geo_problems(bc_admin_lands)
fee_simple_ngo_lands <- fix_geo_problems(fee_simple_ngo_lands)
## Convert IUCN category to an ordered factor
bc_carts$IUCN_CAT <- factor_iucn_cats(bc_carts$IUCN_CAT)
## Union CARTS with itself to eliminate overlaps
bc_carts_agg <- raster::aggregate(bc_carts, by = "PROTDATE")
bc_carts_agg <- fix_geo_problems(bc_carts_agg)
bc_carts_agg_unioned <- self_union(bc_carts_agg) # This takes over a day to run
## Get the earliest year of protection for polygon segments that overlap,
## then aggregate by protected date
bc_carts_agg_unioned$prot_date <- sapply(bc_carts_agg_unioned$union_df, min, na.rm = TRUE)
bc_carts_agg_unioned <- raster::aggregate(bc_carts_agg_unioned, by = "prot_date")
## Union and get attributes of Fee Simple lands
fee_simple_ngo_lands_agg <- raster::aggregate(fee_simple_ngo_lands, by = "prot_date")
fee_simple_ngo_lands_unioned <- self_union(fee_simple_ngo_lands_agg)
fee_simple_ngo_lands_unioned$prot_date <- sapply(fee_simple_ngo_lands_unioned$union_df, min, na.rm = TRUE)
fee_simple_ngo_lands_unioned$prot_date[is.infinite(fee_simple_ngo_lands_unioned$prot_date)] <- NA
fee_simple_ngo_lands_unioned$designation_type <- "NGO Conservation Areas"
fee_simple_ngo_lands_unioned$designation <- "Fee Simple"
fee_simple_ngo_lands_unioned$prot_area <- gArea(fee_simple_ngo_lands_unioned, byid = TRUE)
fee_simple_ngo_lands_agg_unioned <- raster::aggregate(fee_simple_ngo_lands_unioned, by = "prot_date")
## Union and get attributes of BC Administered lands
bc_admin_lands_agg <- raster::aggregate(bc_admin_lands, by = "TENURE_TYPE")
bc_admin_lands_unioned <- self_union(bc_admin_lands_agg)
bc_admin_lands_unioned$prot_date <- max(c(bc_carts_agg_unioned$prot_date,
fee_simple_ngo_lands_unioned$prot_date), na.rm = TRUE)
bc_admin_lands_unioned$designation_type <- "BC Administered Conservation Lands"
bc_admin_lands_unioned$designation <- sapply(bc_admin_lands_unioned$union_df, min, na.rm = TRUE)
bc_admin_lands_unioned$prot_area <- gArea(bc_admin_lands_unioned, byid = TRUE)
bc_admin_lands_agg_unioned <- raster::aggregate(bc_admin_lands_unioned, by = "prot_date")
## Union the first two layers
admin_fee_simple_unioned <- raster::union(bc_admin_lands_agg_unioned, fee_simple_ngo_lands_agg_unioned)
admin_fee_simple_unioned$prot_date <- pmin(admin_fee_simple_unioned$prot_date.1,
admin_fee_simple_unioned$prot_date.2, na.rm = TRUE)
admin_fee_simple_unioned <- raster::aggregate(admin_fee_simple_unioned, by = "prot_date")
## Finally union the BC CARTS data with the previously unioned other layers
prot_areas_unioned <- raster::union(bc_carts_agg_unioned, admin_fee_simple_unioned)
prot_areas_unioned$prot_date <- pmin(prot_areas_unioned$prot_date.1,
prot_areas_unioned$prot_date.2, na.rm = TRUE)
prot_areas_unioned$prot_date[is.infinite(prot_areas_unioned$prot_date)] <- max(c(bc_carts_agg_unioned$prot_date,
fee_simple_ngo_lands_unioned$prot_date), na.rm = TRUE)
## Aggregate the final unioned product by prot_date, so we have a unified layer
## of protected lands and waters for the province, with overlaps removed,
## grouped by date of protection
prot_areas_agg <- raster::aggregate(prot_areas_unioned, by = "prot_date")
save(list = ls(), file = "tmp/prot_areas_clean.rda")
rm(list = ls())
# Process Ecoregions ------------------------------------------------------
## Marine ecoregions
m_ecoregions <- c("HCS", "IPS", "OPS", "SBC", "TPC")
## load ecoregions data from bcmaps package
ecoregions <- bcmaps::ecoregions()
bc_bound_hres <- bcmaps::bc_bound_hres()
## Extract the terrestrial and marine portions of GPB into separate objects
gpb_terrestrial <- ms_clip(ecoregions[ecoregions$CRGNCD == "GPB",],
bc_bound_hres)
gpb_marine <- ms_erase(ecoregions[ecoregions$CRGNCD == "GPB",],
bc_bound_hres)
## Fix it up:
gpb_terrestrial <- fix_geo_problems(gpb_terrestrial)
gpb_marine <- fix_geo_problems(gpb_marine)
## Add terrestrial portion of GPB back to terrestrial ecoregions
ecoregions_t <- rbind(ecoregions[!ecoregions$CRGNCD %in% c("GPB", m_ecoregions), ],
gpb_terrestrial[, setdiff(names(gpb_terrestrial), "rmapshaperid")])
## Add marine portion of GPB back to marine ecoregions
ecoregions_m <- rbind(ecoregions[ecoregions$CRGNCD %in% m_ecoregions, ],
gpb_marine[, setdiff(names(gpb_terrestrial), "rmapshaperid")])
## Calcualte the area of the polygons
ecoregions_t$area <- gArea(ecoregions_t, byid = TRUE)
ecoregions_m$area <- gArea(ecoregions_m, byid = TRUE)
## Create simplified versions for visualization
ecoregions_t_simp <- ms_simplify(ecoregions_t, 0.01)
ecoregions_m_simp <- ms_simplify(ecoregions_m, 0.01)
rm(list = c("gpb_marine", "gpb_terrestrial"))
save(list = ls(), file = "tmp/ecoregions_clean.rda")
rm(list = ls())
# BEC ---------------------------------------------------------------------
## First clip the bec shapefile to the terrestrial BC boundaries
## Using rmapshaper - runs out of memory
# bec <- readOGR("data/BEC_POLY", "BEC_POLY_polygon", stringsAsFactors = FALSE)
# bec_t <- rmapshaper::ms_clip(bec, bc_bound_hres)
## Using mapshaper on the command line. Requires Node installed (https://nodejs.org),
## and install mapshaper with: 'npm install -g mapshaper'
unlink(paste0("data/", c("bc_bound.geojson", "bec_clip*")))
geojsonio::geojson_write(bc_bound_hres, file = "data/bc_bound.geojson") # bc_bound_hres is from bcmaps package
system("mapshaper data/BEC_POLY/BEC_POLY_polygon.shp -clip data/bc_bound.geojson -explode -o data/bec_clip.shp")
## Check for validity of bec polygons
bec_t <- readOGR("data", "bec_clip", stringsAsFactors = FALSE)
if (any(!gIsValid(bec_t, byid = TRUE))) {
bec_t <- gBuffer(bec_t, byid = TRUE, width = 0)
any(!gIsValid(bec_t, byid = TRUE))
}
## Simplify BEC pologyons for use in display
bec_t$poly_id <- row.names(bec_t)
unlink("data/bec_t*")
writeOGR(bec_t, "data", "bec_t", "ESRI Shapefile")
system("mapshaper data/bec_t.shp -simplify 0.01 keep-shapes -o data/bec_t_simp.shp")
bec_t_simp <- readOGR("data", "bec_t_simp", stringsAsFactors = FALSE)
## Repair orphaned holes
bec_t_simp <- gBuffer(bec_t_simp, byid = TRUE, width = 0)
## Put area back in m2
bec_t$area <- gArea(bec_t, byid = TRUE)
bec_t_simp$area <- bec_t$area
## Create a map of bec zones
bec_zone_simp <- raster::aggregate(bec_t_simp, by = "ZONE")
bec_zone_simp$area <- gArea(bec_zone_simp, byid = TRUE)
# bec_t_simp <- ms_simplify(bec_t, keep = 0.01, keep_shapes = TRUE)
save("bec_t", "bec_t_simp", "bec_zone_simp", file = "tmp/bec_clean.rda")