# Load installed libraries
library(raster)
library(leaflet)
library(rgdal)
# import GADM country shapefiles
ken_adm_0 <- raster::getData(name = "GADM",
country = "KEN",
level = 0)
# plot
leaflet() %>%
addProviderTiles("CartoDB.Positron") %>%
addPolygons(data = ken_adm_0,
popup = ken_adm_0$NAME_1,
label = ken_adm_0$NAME_1)
# coverage for southern country area
ken_prec_south <- raster::getData(name = "worldclim",
var = "prec",
res = .5,
lon = 40,
lat = 5)
# coverage for northern country area
ken_prec_north <- raster::getData(name = "worldclim",
var = "prec",
res = .5,
lon = 40,
lat = -5)
# objects are RasterStacks
ken_prec_south
## class : RasterStack
## dimensions : 3600, 3600, 12960000, 12 (nrow, ncol, ncell, nlayers)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : 30, 60, 0, 30 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## names : prec1_27, prec2_27, prec3_27, prec4_27, prec5_27, prec6_27, prec7_27, prec8_27, prec9_27, prec10_27, prec11_27, prec12_27
## min values : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
## max values : 126, 137, 240, 697, 289, 360, 458, 457, 340, 450, 584, 184
ken_prec_north
## class : RasterStack
## dimensions : 3600, 3600, 12960000, 12 (nrow, ncol, ncell, nlayers)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : 30, 60, -30, 0 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## names : prec1_37, prec2_37, prec3_37, prec4_37, prec5_37, prec6_37, prec7_37, prec8_37, prec9_37, prec10_37, prec11_37, prec12_37
## min values : 0, 0, 0, 5, 0, 0, 0, 0, 0, 0, 17, 5
## max values : 520, 466, 563, 746, 473, 279, 275, 218, 172, 350, 458, 464
# subset to extract appropriate layer
ken_prec_north_jan <- ken_prec_north[[1]]
ken_prec_south_jan <- ken_prec_south[[1]]
# merge layers for coverage of country shape area
ken_prec_join_jan <- raster::merge(x = ken_prec_north_jan,
y = ken_prec_south_jan)
# plot to evaluate output of join
plot(ken_prec_join_jan)
# crop precipitcation map to country shape
ken_prec_crop_unmasked <- raster::crop(x = ken_prec_join_jan,
y = ken_adm_0)
# mask to within country shape boundary
ken_prec_crop <- raster::mask(x = ken_prec_crop_unmasked,
mask = ken_adm_0)
# plot to evaluate
plot(ken_prec_crop)
# use aggregate function to change resolution
ken_prec_crop_factup <- raster::aggregate(ken_prec_crop,
# change by a factor of 2 (2 * .5 = 1). Use 'mean' as function to handle continuous values
fact=2,
fun ='mean')
# compare resolution between orginal and aggregated raster layers.
ken_prec_crop
## class : RasterLayer
## dimensions : 1173, 962, 1128426 (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : 33.90833, 41.925, -4.716667, 5.058333 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : layer
## values : 0, 109 (min, max)
ken_prec_crop_factup
## class : RasterLayer
## dimensions : 587, 481, 282347 (nrow, ncol, ncell)
## resolution : 0.01666667, 0.01666667 (x, y)
## extent : 33.90833, 41.925, -4.725, 5.058333 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : layer
## values : 0, 108.75 (min, max)
# divide output of raster by 24.5 to express numerical values as inches.
ken_prec_crop_ag <- ken_prec_crop_factup / 24.5
# create palette to map
raster_colorPal_prec_ag <- colorNumeric(palette = topo.colors(64),
domain = values(ken_prec_crop_ag),
na.color = NA)
# map with legend
leaflet() %>%
addProviderTiles("CartoDB.Positron") %>%
addPolygons(data = ken_adm_0,
popup = ken_adm_0$NAME_0,
label = ken_adm_0$NAME_0,
fillOpacity = 0,
color = "red",
weight = 3,
group = "Kenya") %>%
addRasterImage(x = ken_prec_crop_ag,
color = raster_colorPal_prec_ag) %>%
addLegend(title = "Kenya<br>January Precipitation (in)<br>(1' res)",
values = values(ken_prec_crop_ag),
pal = raster_colorPal_prec_ag, position = "bottomleft")
~fin