chapter 4 - More Raster functions - Part I
sam
4/17/2022
Part I
Let’s load the workspace of the last previous exercis. This is in order to avoid obnoxious problems of ‘a file not being located’.
load('E:/Documents/R Studio is here/R Studio files projects/gis_coding/.RData')Lets overlay our counties shapefile over the elevation raster to have a good feel of which counties seem to have a high altitude range, at least visually.
library(raster)## Warning: package 'raster' was built under R version 4.1.3## Loading required package: splibrary(rgdal)## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
##
## rgdal: version: 1.5-28, (SVN revision 1158)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Users/User/Documents/R/win-library/4.1/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/User/Documents/R/win-library/4.1/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-6
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## Overwritten PROJ_LIB was C:/Users/User/Documents/R/win-library/4.1/rgdal/projlibrary(rgeos)## rgeos version: 0.5-9, (SVN revision 684)
## GEOS runtime version: 3.9.1-CAPI-1.14.2
## Please note that rgeos will be retired by the end of 2023,
## plan transition to sf functions using GEOS at your earliest convenience.
## GEOS using OverlayNG
## Linking to sp version: 1.4-6
## Polygon checking: TRUElibrary(tmap)
library(sp)
library(sf)## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUEhead(county_hs_assets@data)## ADM1_EN Shape_Leng Shape_Area ADM1_PCODE ADM1_REF ADM1ALT1EN ADM1ALT2EN
## 1 Bungoma 3.062486 0.2450581 KE039 <NA> <NA> <NA>
## 2 Kakamega 4.125304 0.2443814 KE037 <NA> <NA> <NA>
## 3 Kiambu 3.193175 0.2066285 KE022 <NA> <NA> <NA>
## 4 Kiambu 3.193175 0.2066285 KE022 <NA> <NA> <NA>
## 5 Nairobi 1.658061 0.0574956 KE047 <NA> <NA> <NA>
## 6 Nakuru 10.598716 1.2836187 KE032 <NA> <NA> <NA>
## ADM0_EN ADM0_PCODE date validOn validTo County...Sub.County
## 1 Kenya KE 2017/11/03 2019/10/31 <NA> BUNGOMA
## 2 Kenya KE 2017/11/03 2019/10/31 <NA> KAKAMEGA
## 3 Kenya KE 2017/11/03 2019/10/31 <NA> KIAMBU
## 4 Kenya KE 2017/11/03 2019/10/31 <NA> KIAMBU
## 5 Kenya KE 2017/11/03 2019/10/31 <NA> NAIROBI CITY
## 6 Kenya KE 2017/11/03 2019/10/31 <NA> NAKURU
## Conventional.Households Stand...alone.Radio Desk.Top.Computer..Laptop..Tablet
## 1 357,714 62.1 4.1
## 2 432,284 62.0 4.3
## 3 46,816 61.2 28.4
## 4 792,333 61.4 19.6
## 5 1,494,676 53.4 22.7
## 6 598,237 61.1 9.6
## Functional.Television Analogue.Television Internet Bicycle Motor.Cycle
## 1 25.4 3.6 7.2 26.2 11.9
## 2 29.2 4.3 8.6 24.4 11.5
## 3 73.8 7.0 49.7 19.2 5.0
## 4 70.0 6.3 39.1 16.3 6.9
## 5 68.7 7.5 42.1 12.5 4.3
## 6 52.5 5.7 20.9 19.1 10.9
## Refrigerator Car Truck..Lorry..Bus..Three.Wheeler.truck Tuk.Tuk names
## 1 2.8 3.1 0.6 0.4 bungoma
## 2 3.3 3.3 0.6 0.4 kakamega
## 3 32.2 21.1 2.3 0.5 kiambu
## 4 19.6 12.4 1.5 0.5 kiambu
## 5 23.5 12.9 1.2 0.4 nairobi city
## 6 8.9 7.1 1.0 0.4 nakuru#show the counties overlain over the elevation data
tm_shape(elevation_counties) + tm_raster(palette = topo.colors(7, 0.7, rev = T), legend.show = T,
legend.is.portrait = T, colorNA = NULL) +
tm_shape(contours) + tm_lines(col = 'violet', lwd = 0.3, breaks = 10) +
tm_text(text = 'level', size = 0.4, col = 'black', along.lines = T, overwrite.lines = T) +
tm_shape(county_hs_assets) + tm_fill(col = NA, alpha = 0) + tm_polygons(col = 'MAP_COLORS', alpha = 0, border.col = 'black', lwd = 0.3, lty = 'solid') + tm_text('ADM1_EN', size = 'AREA', col = 'red', legend.size.show = F, legend.col.show = F) +
tm_layout(title = 'Elevation of 5 most populous counties', legend.position = c('left', 'bottom'))## stars object downsampled to 1080 by 927 cells. See tm_shape manual (argument raster.downsample)## Warning: One tm layer group has duplicated layer types, which are omitted. To
## draw multiple layers of the same type, use multiple layer groups (i.e. specify
## tm_shape prior to each of them).
Now that we’ve been able to see the elevation variation among our counties, we would like to view from which particular threshold of altitude asl most of our counties lie in.
#get the attribute data of our elevation data
elevation_counties## class : RasterLayer
## dimensions : 9332, 10869, 101429508 (nrow, ncol, ncell)
## resolution : 0.0002777778, 0.0002777778 (x, y)
## extent : 34.34347, 37.36264, -1.442083, 1.150139 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : Kenya_SRTM30meters
## values : 1213, 4316 (min, max)#the single digit ie. 1 after the elevation values such as 500, 1000, 1500 etc shows the category
#in which the elevation range falls under
reclassified <- reclassify(elevation_counties, c(-Inf, 500, 1, 500, 1000, 2,
1000, 1500, 3, 1500, 2000, 4, 2000, 2500, 5, 2500, 3000, 6, 3000, 3500, 7, 3500, 4000, 8, 4000, 4500, 9))#create a table showing the sum of pixels under each category of elevation value
ranges <- table(as.vector(reclassified))
ranges##
## 3 4 5 6 7 8 9
## 3209998 8148594 4133717 1892033 183095 64698 16644Having the sum of pixels for each category makes little sense other than seeing a tyranny of numbers at play. We will divide the sum of pixels under each category with the total sum of pixels, multiplied by 100. This will give us the percentage of area held by each elevation category.
(ranges/sum(ranges)) * 100##
## 3 4 5 6 7 8
## 18.18821574 46.17086542 23.42211322 10.72047534 1.03743721 0.36658627
## 9
## 0.09430681According to the statistics derived above, the largest elevation range for our five counties lies within group 4, which according to the reclassified function we used, stands for areas between 1500 to 2000m asl. Having a high status is not so bad after all. You may ask what happened to categories 1 and 2. Well, their elevation ranges do not exist within our counties.
Now that we have seen the percentage of various ranges of elevation values across our entire dataset, we will go a step further to investigate the variation within the counties. ie. max min, mean and other statistical measurements for the elevation values disaggregated to each of the five county units. Makes sense? If not, just follow along.
#first install theexactextractr and thereafter load into R
library(exactextractr)## Warning: package 'exactextractr' was built under R version 4.1.3#the function exact_extract returns a value of the summary of the cells in a raster covered by
#the polygons
exact_extract(elevation_counties, county_hs_assets, c('min', 'max', 'mean', 'variance'))##
|
| | 0%
|
|============ | 17%
|
|======================= | 33%
|
|=================================== | 50%
|
|=============================================== | 67%
|
|========================================================== | 83%
|
|======================================================================| 100%## min max mean variance
## 1 1213 4316 1852.509 361509.03
## 2 1213 2113 1519.226 30534.68
## 3 1279 2727 1887.173 119355.87
## 4 1279 2727 1887.173 119355.87
## 5 1450 1951 1651.480 13358.66
## 6 1468 3082 2139.271 122546.16But we are missing the names of the counties. So lets bind this dataset to our counties. But not after getting the answer this time round. We will become more sophisticated by performing the entire operation and binding it to our county_hs_assets all in one line of code.
library(raster)
library(rgdal)
library(rgeos)
library(tmap)
library(sp)
library(sf)#cbind does the magic of combining two datasets, by columns or rows. In this case by columns
county_elevation <- cbind(county_hs_assets, exact_extract(elevation_counties, county_hs_assets, c('min', 'max', 'mean')))##
|
| | 0%
|
|============ | 17%
|
|======================= | 33%
|
|=================================== | 50%
|
|=============================================== | 67%
|
|========================================================== | 83%
|
|======================================================================| 100%county_elevation## class : SpatialPolygonsDataFrame
## features : 6
## extent : 34.34341, 37.3626, -1.442148, 1.150185 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 29
## names : ADM1_EN, Shape_Leng, Shape_Area, ADM1_PCODE, ADM1_REF, ADM1ALT1EN, ADM1ALT2EN, ADM0_EN, ADM0_PCODE, date, validOn, validTo, County...Sub.County, Conventional.Households, Stand...alone.Radio, ...
## min values : Bungoma, 1.65806105822, 0.0574955973101, KE022, NA, NA, NA, Kenya, KE, 2017/11/03, 2019/10/31, NA, BUNGOMA, 1,494,676, 53.4, ...
## max values : Nakuru, 10.5987160624, 1.28361868572, KE047, NA, NA, NA, Kenya, KE, 2017/11/03, 2019/10/31, NA, NAKURU, 792,333, 62.1, ...county_elevation@data## ADM1_EN Shape_Leng Shape_Area ADM1_PCODE ADM1_REF ADM1ALT1EN ADM1ALT2EN
## 1 Bungoma 3.062486 0.2450581 KE039 <NA> <NA> <NA>
## 2 Kakamega 4.125304 0.2443814 KE037 <NA> <NA> <NA>
## 3 Kiambu 3.193175 0.2066285 KE022 <NA> <NA> <NA>
## 4 Kiambu 3.193175 0.2066285 KE022 <NA> <NA> <NA>
## 5 Nairobi 1.658061 0.0574956 KE047 <NA> <NA> <NA>
## 6 Nakuru 10.598716 1.2836187 KE032 <NA> <NA> <NA>
## ADM0_EN ADM0_PCODE date validOn validTo County...Sub.County
## 1 Kenya KE 2017/11/03 2019/10/31 <NA> BUNGOMA
## 2 Kenya KE 2017/11/03 2019/10/31 <NA> KAKAMEGA
## 3 Kenya KE 2017/11/03 2019/10/31 <NA> KIAMBU
## 4 Kenya KE 2017/11/03 2019/10/31 <NA> KIAMBU
## 5 Kenya KE 2017/11/03 2019/10/31 <NA> NAIROBI CITY
## 6 Kenya KE 2017/11/03 2019/10/31 <NA> NAKURU
## Conventional.Households Stand...alone.Radio Desk.Top.Computer..Laptop..Tablet
## 1 357,714 62.1 4.1
## 2 432,284 62.0 4.3
## 3 46,816 61.2 28.4
## 4 792,333 61.4 19.6
## 5 1,494,676 53.4 22.7
## 6 598,237 61.1 9.6
## Functional.Television Analogue.Television Internet Bicycle Motor.Cycle
## 1 25.4 3.6 7.2 26.2 11.9
## 2 29.2 4.3 8.6 24.4 11.5
## 3 73.8 7.0 49.7 19.2 5.0
## 4 70.0 6.3 39.1 16.3 6.9
## 5 68.7 7.5 42.1 12.5 4.3
## 6 52.5 5.7 20.9 19.1 10.9
## Refrigerator Car Truck..Lorry..Bus..Three.Wheeler.truck Tuk.Tuk names
## 1 2.8 3.1 0.6 0.4 bungoma
## 2 3.3 3.3 0.6 0.4 kakamega
## 3 32.2 21.1 2.3 0.5 kiambu
## 4 19.6 12.4 1.5 0.5 kiambu
## 5 23.5 12.9 1.2 0.4 nairobi city
## 6 8.9 7.1 1.0 0.4 nakuru
## min max mean
## 1 1213 4316 1852.509
## 2 1213 2113 1519.226
## 3 1279 2727 1887.173
## 4 1279 2727 1887.173
## 5 1450 1951 1651.480
## 6 1468 3082 2139.271#remove unnecessary columns to get a clean table for analysis
county_elevation2 <- subset(county_elevation@data, select = -c(2:11))#read county_elevatio2 attributes
county_elevation2## ADM1_EN validTo County...Sub.County Conventional.Households
## 1 Bungoma <NA> BUNGOMA 357,714
## 2 Kakamega <NA> KAKAMEGA 432,284
## 3 Kiambu <NA> KIAMBU 46,816
## 4 Kiambu <NA> KIAMBU 792,333
## 5 Nairobi <NA> NAIROBI CITY 1,494,676
## 6 Nakuru <NA> NAKURU 598,237
## Stand...alone.Radio Desk.Top.Computer..Laptop..Tablet Functional.Television
## 1 62.1 4.1 25.4
## 2 62.0 4.3 29.2
## 3 61.2 28.4 73.8
## 4 61.4 19.6 70.0
## 5 53.4 22.7 68.7
## 6 61.1 9.6 52.5
## Analogue.Television Internet Bicycle Motor.Cycle Refrigerator Car
## 1 3.6 7.2 26.2 11.9 2.8 3.1
## 2 4.3 8.6 24.4 11.5 3.3 3.3
## 3 7.0 49.7 19.2 5.0 32.2 21.1
## 4 6.3 39.1 16.3 6.9 19.6 12.4
## 5 7.5 42.1 12.5 4.3 23.5 12.9
## 6 5.7 20.9 19.1 10.9 8.9 7.1
## Truck..Lorry..Bus..Three.Wheeler.truck Tuk.Tuk names min max
## 1 0.6 0.4 bungoma 1213 4316
## 2 0.6 0.4 kakamega 1213 2113
## 3 2.3 0.5 kiambu 1279 2727
## 4 1.5 0.5 kiambu 1279 2727
## 5 1.2 0.4 nairobi city 1450 1951
## 6 1.0 0.4 nakuru 1468 3082
## mean
## 1 1852.509
## 2 1519.226
## 3 1887.173
## 4 1887.173
## 5 1651.480
## 6 2139.271You may have noticed Kiambu is duplicated but this is a pesky problem with long history. So we will look over it for now so long as it hardly affects our analysis. I take the blame.
Now that we have seen the variation of altitude for each of the five counties, we can do some slightly more complex analysis whereby we extract the types of landcover in each county.
Remember the countrywide elevation raster we tried to load in chapter 3? Guess what, the landcover dataset that we will use is actually heavy duty staff and we are not willing to strain our laptop RAM again. Because of its sheer size, we will load the landcover raster and clip it straight away to fit into our ‘parameters’. That is, to the bounds and boundaries of our counties.
The landcover raster is available from here: http://geoportal.rcmrd.org/data/Kenya_Sentinel2_LULC2016.zip
#load the landcover raster into R environment
landcover <- raster('D:/DUKTO/Kenya_Sentinel2_LULC2016new/Kenya_Sentinel2_LULC2016.tif')#take a quickie look at its metadata
landcover## class : RasterLayer
## dimensions : 54719, 43171, 2362273949 (nrow, ncol, ncell)
## resolution : 0.0001851853, 0.0001851853 (x, y)
## extent : 33.91169, 41.90632, -4.702468, 5.430684 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : Kenya_Sentinel2_LULC2016.tif
## names : Kenya_Sentinel2_LULC2016
## values : 1, 10 (min, max)To avoid long loading times when plotting, we will crop the landcover raster into the extremities of our counties dataset. Thereafter, we will mask it to the boundaries of our dataset, thus humbling it to the shape of our boundaries. Think of masking like cutting a piece of paper to a potrait you love. Thereon, we can be on good coding terms.
#crop the landcover raster to the bounds of the counties
landcover_counties <- crop(landcover, adm_subset)#mask the landcover to the counties boundaries
landcover_countiesmap <- mask(landcover_counties, adm_subset)Now that took quite some time. Was it the same for you? Let’s see if the raster was indeed loaded into R.
#view raster metadata
landcover_countiesmap## class : RasterLayer
## dimensions : 13998, 16304, 228223392 (nrow, ncol, ncell)
## resolution : 0.0001851853, 0.0001851853 (x, y)
## extent : 34.34335, 37.36261, -1.442097, 1.150127 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : r_tmp_2022-04-23_191751_1200_24967.grd
## names : Kenya_Sentinel2_LULC2016
## values : 1, 10 (min, max)#nota betta. save the raster into the laptop's memory.
writeRaster(landcover_countiesmap, 'landcover_countiesmap.tiff', overwrite = T)Note the NB! My experience has shown me that R will not recognize the landcover_countiesmap as soon as you quite R. Just save it in case problems arise.
library(tmap)I tried to use the below code to plot the landcover_countiesmap.
#plot the landcover_countiesmap
tm_shape(landcover_countiesmap) + tm_raster(palette = terrain.colors(10, 0.7, rev = F), legend.show = T, legend.is.portrait = T, colorNA = NULL) +
tm_layout(title = 'Landcover types of top 5 counties by population', legend.position = c('left', 'bottom'))The below error was the result.
Warning in fetch(.x, ..., downsample = downsample) :
with RasterIO defined, argument downsample is ignored
stars_proxy object shown at 16304 by 13998 cells.
Error: cannot allocate vector of size 870.6 Mb
Error during wrapup: cannot allocate vector of size 870.6 Mb
Error: no more error handlers available (recursive errors?); invoking 'abort' restartIt’s good practice to visually check, through an image if possible, if your raster dataset has been loaded, and is perfect for that matter. Judging by the error above, on the other hand, R is also plotting that I don’t sleep today.
library(raster)To reduce the landcover_countiesmap to a memory size that allows tmap to plot it, we will reduce the number of cells in the raster using the aggregate function.
#reduce the raster by 3 cells in each direction, and replacing that with one cell of the modal #cell values in both directions
countieslandcover <- aggregate(landcover_countiesmap, fact = 3, fun = modal, expand = F, na.rm = T)#view metadata of our aggregated landcover dataset
countieslandcover## class : RasterLayer
## dimensions : 4666, 5434, 25355044 (nrow, ncol, ncell)
## resolution : 0.0005555558, 0.0005555558 (x, y)
## extent : 34.34335, 37.36224, -1.442097, 1.150127 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : r_tmp_2022-04-23_191924_1200_71753.grd
## names : Kenya_Sentinel2_LULC2016
## values : 1, 10 (min, max)#save the aggregated raster
writeRaster(countieslandcover, 'countieslandcover.tiff', overwrite = T)tm_shape(countieslandcover) + tm_raster(palette = terrain.colors(10, 0.7, rev = F), n= 10, legend.show = T, legend.is.portrait = T, colorNA = NULL) +
tm_layout(title = 'Landcover types of top 5 counties by population', legend.position = c('left', 'bottom'))## stars object downsampled to 1079 by 927 cells. See tm_shape manual (argument raster.downsample)
A bitter-sweet success after the aggregation. A bit of confusion though. I expected classes to be individually listed as 1 to 10 with each class having a distinct number. Oh well. We will work with this.
Lets try to reproject it to a local UTM datum and hopefully see if the classes are distinct from each other.
library(rgdal)#reprojeting the aggregated raster to a local datum WGS 84/UTM Zone 37S
countiesland_datum <- projectRaster(countieslandcover, crs = '+init=epsg:32737', res = 200, method = 'bilinear', over = F, filename = 'countiesland_datum.tiff', overwrite = T)countiesland_datum## class : RasterLayer
## dimensions : 1447, 1693, 2449771 (nrow, ncol, ncell)
## resolution : 200, 200 (x, y)
## extent : -19743.53, 318856.5, 9839147, 10128547 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=37 +south +datum=WGS84 +units=m +no_defs
## source : countiesland_datum.tiff
## names : countiesland_datum
## values : 1, 10 (min, max)tm_shape(countiesland_datum) + tm_raster(palette = terrain.colors(10, 0.7, rev = F), n= 10, legend.show = T, legend.is.portrait = T, colorNA = NULL) +
tm_layout(title = 'Landcover types of top 5 counties by population', legend.position = c('left', 'bottom'))## stars object downsampled to 1082 by 924 cells. See tm_shape manual (argument raster.downsample)
All the same. However, I asked the experts at StackOverflow for a solution to this peculiar problem. However, by the time of going to print, I had not seen a new working solution yet. See the response here: https://stackoverflow.com/questions/71974140/raster-classes-not-appearing-in-legend/71977381#71977381
#create a duplicate raster dataset to avoid messing with original. Good practice
countiesland_datum2 <- countiesland_datum#try to restore the levels of the new countiesland_datum with those of the original landcover #raster before aggregation
levels(countiesland_datum2) <- levels(landcover_countiesmap)tm_shape(countiesland_datum2) + tm_raster(palette = terrain.colors(10, 0.7, rev = F), n= 10, legend.show = T, legend.is.portrait = T, colorNA = NULL) +
tm_layout(title = 'Landcover types of top 5 counties by population', legend.position = c('left', 'bottom'))## stars object downsampled to 1082 by 924 cells. See tm_shape manual (argument raster.downsample)
All the same again. However, big-ups to StackOverflow for its contribution to easing the life of coders.