Thanks for the following blog, which teaches me how to extract nightlight data from .tif file in R:
The night light data can be downloaded at https://eogdata.mines.edu/products/vnl/ , choose Monthly Cloud-free DNB Composite. I just downloaded the nightlight data in January 2019, and I want to extract the Uganda night light data, because, for many developing countries, sub-regional economic data are not available sometimes, the night light data is a good proxy. In addition, 00N060W (south part) and 75N060W (north part) cover the whole of Uganda.
The problem with Uganda's night light data is that Uganda is an equatorial country, and two night light files cover its land. Later in my code part, you can find that the `merge` function could fix this. Below is a visualization of Uganda's night light data in January 2019, most are shown in yellow just because they are very close to 0, meanwhile, we can find that near the capital of Uganda, Kampala, the value is much higher.
library(raster)
library(sf)
library(haven)
library(tidyr)
library(ggplot2)
filepath <- "D:/Master/[Master Thesis]/night light/SVDNB_npp_20190101-20190131_00N060W_vcmcfg_v10_c201905201300/test2.tif"
filepath2 <- "D:/Master/[Master Thesis]/night light/SVDNB_npp_20190101-20190131_75N060W_vcmcfg_v10_c201905201300/test3.tif"
data <- raster(filepath)
data2 <- raster(filepath2)
uga_shape<- getData('GADM', country='uga', level=3)
uga_sf <- as(uga_shape, "sf")
nl_data_uganda_N <- crop(data2, uga_sf)
nl_data_uganda_S <- crop(data, uga_sf)
merged_raster <- merge(nl_data_uganda_N, nl_data_uganda_S)
extract_ugadata2 <- raster::extract(merged_raster, uga_shape, fun=mean, df=TRUE, na.rm = TRUE)
final_data2 <- cbind(st_drop_geometry(uga_sf), extract_ugadata2)
uganda_shapefile_cropped <- crop(uga_shape, merged_raster)
nl_data_df <- data.frame(lon = coordinates(uganda_shapefile_cropped)[, 1],
lat = coordinates(uganda_shapefile_cropped)[, 2],
nl_mean = drop_na(extract_ugadata2))
ggplot() +
geom_polygon(data = fortify(uganda_shapefile_cropped),
aes(x = long, y = lat, group = group), fill = "gray",color="black") +
geom_point(data = nl_data_df, aes(x = lon, y = lat, color = nl_mean.layer), size = 3) +
scale_color_gradient(low = "yellow", high = "red") +
coord_equal()
Comments