Conversion Hotspots
Alex Zvoleff
2024-06-08
hotspots <- read_csv("conversionHotspots.csv")## Rows: 824 Columns: 23
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): ecoregion
## dbl (22): natShortVeg00, natShortVeg10, totalArea, tree00, tree10, treetoCro...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.ecoregions_raw <- st_read("ecoregions/Ecoregions2017.shp")## Reading layer `Ecoregions2017' from data source
## `C:\Users\azvol\OneDrive - Conservation International Foundation\Code\LandDegradation\sbtn-land-hub\conversion_hotspots\ecoregions\Ecoregions2017.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 847 features and 15 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -89.89197 xmax: 180 ymax: 83.62313
## Geodetic CRS: WGS 84ecoregions_raw <- st_make_valid(ecoregions_raw)
ecoregions <- st_simplify(ecoregions_raw, dTolerance = 1000)
ecoregions <- left_join(ecoregions, hotspots, by=join_by(ECO_NAME==ecoregion))
# Check that there are no NAs (i.e. that the join worked correctly)
table(is.na(ecoregions$allVegConversion00_10))##
## FALSE TRUE
## 823 24table(is.na(ecoregions$allVegConversion10_20))##
## FALSE TRUE
## 823 24Setup indicators of top ecoregions in terms of conversion
Look at top 10% of ecoregions in terms of recently converted area (from 2010 - 2020), and also in terms of total conversion since 2000. Also add an indicator of ecoregions where conversion has accelerated (i.e. higher rates of conversion in 2010-2020 version 2000-2010).
ecoregions %>%
filter(!is.na(allVegArea00)) %>%
mutate(accelerating=ifelse(allVegConversion10_20 > (allVegConversion00_10*1.1),
'Accelerating', 'Decelerating'),
topByAreaRecent=percent_rank(allVegConversion10_20)>.9,
topByPctRecent=percent_rank(allVegConversion10_20/(allVegArea00 - allVegConversion00_10))>.9,
topByAreaOverall=percent_rank(totalAllVegConversion00_20)>.9,
topByPctOverall=percent_rank(totalAllVegConversion00_20/allVegArea00)>.9) -> ecoregionsPlot some basic stats
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 38 rows containing non-finite outside the scale range
## (`stat_bin()`).
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 38 rows containing non-finite outside the scale range
## (`stat_bin()`).
## Warning: Removed 38 rows containing non-finite outside the scale range
## (`stat_ecdf()`).
What percentage of total conversion does each method capture?
# Recent conversion
sum(ecoregions$allVegConversion10_20[ecoregions$topByPctRecent], na.rm=TRUE) / sum(ecoregions$allVegConversion10_20) * 100## [1] 34.79098sum(ecoregions$allVegConversion10_20[ecoregions$topByAreaRecent], na.rm=TRUE) / sum(ecoregions$allVegConversion10_20) * 100## [1] 69.45636# Overall conversion
sum(ecoregions$totalAllVegConversion00_20[ecoregions$topByPctOverall], na.rm=TRUE) / sum(ecoregions$totalAllVegConversion00_20) * 100## [1] 38.95094sum(ecoregions$totalAllVegConversion00_20[ecoregions$topByAreaOverall], na.rm=TRUE) / sum(ecoregions$totalAllVegConversion00_20) * 100## [1] 71.24139How many ecoregions show up in multiple methods
ecoregions$nTop <- ecoregions$topByPctRecent + ecoregions$topByAreaRecent + ecoregions$topByPctOverall + ecoregions$topByAreaOverall
table(ecoregions$nTop)##
## 0 1 2 3 4
## 636 31 82 15 21And what do maps look like
By area, last 10 years
ggplot() +
geom_sf(data=world_wintri) +
geom_sf(data=filter(ecoregions, topByAreaRecent), aes(fill=ECO_NAME)) +
geom_sf(data=grat_wintri, color = "gray90", size = 0.05/.pt) +
coord_sf(datum=st_crs("+proj=wintri +datum=WGS84 +no_defs +over")) +
theme_map() +
theme(legend.position="none")
By percent converted, last 10 years
ggplot() +
geom_sf(data=world_wintri) +
geom_sf(data=filter(ecoregions, topByPctRecent), aes(fill=ECO_NAME)) +
geom_sf(data=grat_wintri, color = "gray90", size = 0.05/.pt) +
coord_sf(datum=st_crs("+proj=wintri +datum=WGS84 +no_defs +over")) +
theme_map() +
theme(legend.position="none")
By area converted, last 20 years
ggplot() +
geom_sf(data=world_wintri) +
geom_sf(data=filter(ecoregions, topByAreaOverall), aes(fill=ECO_NAME)) +
geom_sf(data=grat_wintri, color = "gray90", size = 0.05/.pt) +
coord_sf(datum=st_crs("+proj=wintri +datum=WGS84 +no_defs +over")) +
theme_map() +
theme(legend.position="none")
By percent converted, last 20 years
ggplot() +
geom_sf(data=world_wintri) +
geom_sf(data=filter(ecoregions, topByPctOverall), aes(fill=ECO_NAME)) +
geom_sf(data=grat_wintri, color = "gray90", size = 0.05/.pt) +
coord_sf(datum=st_crs("+proj=wintri +datum=WGS84 +no_defs +over")) +
theme_map() +
theme(legend.position="none")
The 21 ecoregions that show up across methods
ecoregions %>%
filter(nTop==4) %>%
arrange(ECO_NAME) %>%
st_drop_geometry() %>%
select(ECO_NAME)## ECO_NAME
## 1 Alto Paraná Atlantic forests
## 2 Aravalli west thorn scrub forests
## 3 Canadian Aspen forests and parklands
## 4 Central Deccan Plateau dry deciduous forests
## 5 Central-Southern US mixed grasslands
## 6 Cerrado
## 7 Chhota-Nagpur dry deciduous forests
## 8 Deccan thorn scrub forests
## 9 Humid Chaco
## 10 Humid Pampas
## 11 Khathiar-Gir dry deciduous forests
## 12 Maranhão Babaçu forests
## 13 Narmada Valley dry deciduous forests
## 14 Pontic steppe
## 15 Tamaulipan mezquital
## 16 Thar desert
## 17 Tocantins/Pindare moist forests
## 18 Upper Gangetic Plains moist deciduous forests
## 19 Uruguayan savanna
## 20 Victoria Basin forest-savanna
## 21 West Sudanian savannaggplot() +
geom_sf(data=world_wintri) +
geom_sf(data=filter(ecoregions, nTop==4), aes(fill=ECO_NAME)) +
geom_sf(data=grat_wintri, color = "gray90", size = 0.05/.pt) +
coord_sf(datum=st_crs("+proj=wintri +datum=WGS84 +no_defs +over")) +
theme_map() +
theme(legend.position="none")
gadm <- st_read("C:/Users/azvoleff/OneDrive - Conservation International Foundation/Data/gadm_410-levels.gpkg", "ADM_1")## Reading layer `ADM_1' from data source
## `C:\Users\azvoleff\OneDrive - Conservation International Foundation\Data\gadm_410-levels.gpkg'
## using driver `GPKG'
## Simple feature collection with 3662 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -55.98 xmax: 180 ymax: 83.65833
## Geodetic CRS: WGS 84gadm_0 <- st_read("C:/Users/azvoleff/OneDrive - Conservation International Foundation/Data/gadm_410-levels.gpkg", "ADM_0")## Reading layer `ADM_0' from data source
## `C:\Users\azvoleff\OneDrive - Conservation International Foundation\Data\gadm_410-levels.gpkg'
## using driver `GPKG'
## Simple feature collection with 263 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.65833
## Geodetic CRS: WGS 84gadm <- st_make_valid(gadm)
ecoregions$hotspot <- ecoregions$nTop==4
tic()
gadm %>%
st_filter(y=filter(ecoregions, hotspot), .predicate = st_intersects) %>%
mutate(admin_area = st_area(.)) -> hotspot_admin
gadm_intersection <- st_intersection(hotspot_admin, filter(ecoregions, hotspot))## Warning: attribute variables are assumed to be spatially constant throughout
## all geometriestoc()## 237.48 sec elapsedfull_join(
hotspot_admin,
gadm_intersection %>%
mutate(intersect_area = st_area(.)) %>%
select(GID_1, intersect_area) %>%
st_drop_geometry() %>%
group_by(GID_1) %>%
summarise(intersect_area=sum(intersect_area))
) %>%
mutate(intersect_frac=as.numeric(intersect_area / admin_area)) -> intersection_frac## Joining with `by = join_by(GID_1)`Jurisdictions overlapping the 21 ecoregions that show up across methods
Jurisdictions by percent overlap
ggplot() +
geom_sf(data=world_wintri) +
geom_sf(data=intersection_frac, aes(fill=intersect_frac > .5)) +
geom_sf(data=grat_wintri, color = "gray90", size = 0.05/.pt) +
coord_sf(datum=st_crs("+proj=wintri +datum=WGS84 +no_defs +over")) +
theme_map()
Jurisdictions with at least 50% overlap overlaid on hotspot ecoregions
ggplot() +
geom_sf(data=world_wintri) +
geom_sf(data=filter(ecoregions, hotspot), fill='cyan') +
geom_sf(data=filter(intersection_frac, intersect_frac > .5), fill='blue') +
geom_sf(data=grat_wintri, color = "gray90", size = 0.05/.pt) +
coord_sf(datum=st_crs("+proj=wintri +datum=WGS84 +no_defs +over")) +
theme_map()