-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathagriculture_prep_accessibility.R
152 lines (100 loc) · 5.59 KB
/
agriculture_prep_accessibility.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
# ------------- PREPARE WORKSPACE -------------
setwd("C:/Users/grace/GIS/povertyequity/chad")
wd <- getwd()
library(ursa)
library(raster)
library(sp)
library(rgdal)
library(sf)
library(dplyr)
library(rgeos)
library(naniar) # replace_with_na()
library(tidyverse) # ggplot2()
library(tidyr) # complete()
library(units)
# devtools:::install_github("gearslaboratory/gdalUtils")
library(gdalUtils)
# ------------- MOSAIC AND CLIP LANDCOVER -------------
# Tutorial: https://tengkengvang.com/2018/11/12/mosaic-or-merge-rasters-in-r/
landcover <- c('GLCLU_2019_strata_10N_010E.tif', 'GLCLU_2019_strata_10N_020E.tif',
'GLCLU_2019_strata_20N_010E.tif', 'GLCLU_2019_strata_20N_020E.tif',
'GLCLU_2019_strata_30N_010E.tif', 'GLCLU_2019_strata_30N_020E.tif')
e = extent(13, 25, 7, 24) # Approximate extents of Chad
template = raster(e) # Create empty raster to add actual data onto
proj4string(template) <- CRS('+init=epsg:4326')
writeRaster(template, file="TCD_landcover.tif", format="GTiff", overwrite=TRUE)
mosaic_rasters(gdalfile=landcover, dst_dataset="TCD_landcover.tif", of="GTiff")
gdalinfo("TCD_landcover.tif") # Check contents.
# Load raster and AOI extent as objects for clip.
chad = readOGR("TCD_bbox.shp")
landcover = raster("TCD_landcover.tif")
plot(landcover, main = "Original mosaic")
plot(chad, add=TRUE)
landcover <- crop(landcover, chad)
plot(landcover, main = "Cropped to Chad bounding box")
plot(chad, add=TRUE)
writeRaster(landcover, file="TCD_landcover.tif", format="GTiff", overwrite=TRUE)
# ------------- SUBSET TO JUST CROPLAND VALUE OF LANDCOVER BAND -------------
cropland = clamp(landcover, 17, 17, useValues=FALSE) # Set all non-cropland values to null. (values below and above 17)
cropland = reclassify(cropland, c(16,18,1)) # Reclassify all values between 16-18 to 1.
writeRaster(cropland, file="TCD_cropland.tif", format="GTiff", overwrite=TRUE)
# ------------- GET CROPLAND CELL COUNT AT AGRO'S RESOLUTION (TO GET AGRO VALUE PER CROP CELL) -------------
# Get details of each dataset to:
# Determine by what factor to aggregate cropland.
# Determine which has the larger extent. Use extend() if agro is larger, or crop() if cropland is larger.
gdalinfo("TCD_spam17sum.tif")
gdalinfo("TCD_cropland.tif")
agro = raster("TCD_spam17sum.tif")
cropland = raster("TCD_cropland.tif")
# Add up # cells at same resolution as agro (9km / 30m = 333.33)
cropland_9km = raster::aggregate(cropland, fact=333, fun=sum, na.rm=TRUE)
# Give the same extent by adding NA rows/cols wherever needed (if extending) or removing rows/cols (if cropping)
cropland_9km = raster::extend(cropland_9km, agro)
plot(cropland_9km)
plot(agro, add=TRUE)
agro_aligned = resample(agro, cropland_9km, method="bilinear")
plot(agro_aligned, add=TRUE)
# The input Raster* objects should have the same extent, origin and resolution.
# If only the extent differs, the computation will continue for the intersection of the Raster objects.
agro_per_cropcell = agro_aligned / cropland_9km
plot(agro_per_cropcell)
# ------------- GET PROPORTIONAL AGRO AMOUNT AT CROPLAND'S RESOLUTION (TO ASSIGN AGRO VALUE TO CROP CELL) -------------
# Split agro cells to resolution of cropland, keeping same value as its parent cell without arithmetic.
agro_30m = raster::disaggregate(agro_per_cropcell, fact=333, method='', na.rm=TRUE)
agro_30m = raster::extend(agro_30m, cropland)
agro_aligned = resample(agro_30m, cropland, method="bilinear")
val_at_cropcells = agro_aligned * cropland
writeRaster(val_at_cropcells, file="TCD_AgroVal_byCroplandCell.tif", format="GTiff", overwrite=TRUE)
# ------------- CROPLAND VALUES TO POINT FOR USE IN TIME-COST RASTER -------------
crop_pt = rasterToPoints(val_at_cropcells, fun=NULL, spatial=TRUE)
shapefile(crop_pt, "TCD_AgroVal_byCroplandCell.shp", overwrite=TRUE)
# ------------- AGRICULTURAL VALUES TO POLYGON -------------
agro = raster("spam2017v2r1_V_sum1.tif")
agro = crop(agro, chad)
writeRaster(agro, file="TCD_spam17sum.tif", format="GTiff", overwrite=TRUE)
agro_ply = rasterToPolygons(agro, fun=NULL, n=4, na.rm=TRUE, digits=12, dissolve=FALSE)
plot(agro_ply)
shapefile(agro_ply, "TCD_spam17sum.shp", overwrite=TRUE)
# ------------- COUNT CROP CELLS IN EACH AGRO CELL (NOW POLYGON) -------------
cropland = raster("TCD_cropland.tif")
agro = raster("TCD_spam17sum.tif")
agro_ply = readOGR("TCD_spam17sum.shp")
#agro_ply$crop_sum = raster::extract(cropland, agro_ply, fun=sum, sp=TRUE) # Add raster:: to start to avoid conflict with other packages that have extract() function.
shapefile(agro_ply, "TCD_spam17sum_cropct.shp", overwrite=TRUE)
# ------------- CROPLAND VALUES TO POLYGON -------------
crop_ply = rasterToPolygons(cropland)
#crop_pt = as.data.frame(crop_pt)
#crop_pt$index = 1:nrow(crop_pt)
#crop_pt = rename(crop_pt, crop_ID = index)
#write.csv(crop_pt, file = "TCD_cropland_pt.csv")
shapefile(crop_ply, "TCD_cropland_ply.shp", overwrite=TRUE)
# ------------- JOIN AGRO VALUES and ID'S TO CROPLAND POINTS -------------
# Was having processing trouble, so switched to QGIS to lighten the load a bit.
# In QGIS, re-projected crop_pt and agro to Africa Albers.
chad = readOGR("TCD_ADM0_10km.shp")
crop_pt = readOGR("TCD_cropland_pt_albers.shp")
agro = readOGR("TCD_spam17sum_albers.shp")
agro = st_read("TCD_spam17sum.shp") # Re-read as sf data type to be compatible with st_join
crop_pt = st_read("TCD_cropland_pt.shp")
agro_crop = st_join(crop_pt, agro, join=st_within, left=T)
st_write(agro_crop, "TCD_crop_agrovalues.shp", driver="ESRI Shapefile", append=FALSE)