Elías Cisneros

Code: Nightlight (DMSP-OLS) download and preparation

Download and pre-processing of nightlight data DMSP-OLS versions (1992-2013).

Setup

library(here)
library(raster)
library(clipr)
library(ggplot2)
library(rgeos)
system("mkdir dload")
system("mkdir data")

Download

This will download the nightlihgt data in the background and show the progress on your R console. I do prefer to copy the curl commands (stored as a) and paste them in the terminal, where it is more convinient to observe the download process. Alternativly you can also just download the files manually on the United States National Oceanic and Atmospheric Administration (NOAA) webpage: https://ngdc.noaa.gov/eog/data/web_data/v4composites. Note: The download takes some time, therefore the code checks for already existing files.

#available nightlight versions
l<-c(
  paste0("F10",1992:1994),
  paste0("F12",1994:1999),
  paste0("F14",1997:2003),
  paste0("F15",2000:2007), #    paste0("F15",2008), # a hidden 2008 file exists
  paste0("F16",2004:2009),
  paste0("F18",2010:2013)
)
#choose nightlight versions
l<-l[grep("2007|2008",l)]
#exlude already existing files
f<-list.files("dload",".tif|.tar")
f<-sub("\\.v4\\.tar","",f)
f<-sub(".v4[b|d|c]_web.stable_lights.avg_vis.tif","",f)
f<-unique(f)
#download list
l<-l[!(l %in% f)]
dlist<-paste0("https://ngdc.noaa.gov/eog/data/web_data/v4composites/",l,".v4.tar")
##download with curl - available also for Windows
a<-lapply(
  dlist,
  function(h){
    # h<-dlist[2]
    fname<-substring(h,max(c(gregexpr("/",h)[[1]]))+1)
    if (file.exists(fname)!=T){
      cmd<-paste0("curl ",h," --output ",here("dload",fname))
    } else (
      cmd<-''
    )
  }
)
a
write_clip(unlist(a),allow_non_interactive=T) #and paste it to the terminal or run:
# system(unlist(a),wait=F)

Unzip nightlight data

The downloaded tar files contain several nightlight composits with different corrections. Which you can see when you unzip the F1?YYY.v4.tar file. The tar file contains several differetn gunziped files such as: F1?YYYY_v4?_stable_lights.avg_vis.tif.gz. which is one version of the nightlight products, which already corrected for clouds and fires (for more detail read the README_V?.txt of your downloads). This file has to be unziped again with gunzip (also available for Windows).

list.files("dload","F.*\\.tar")
flist<-list.files(here("dload"),".v4.tar$",full.names=T)
# flist

for (f in flist){
  # f<-flist[1]
  print(f)
  #list files in tar file
  x<-system(paste0("tar -tf ",f),intern=T)
  #choose file to unzip
  s<-x[grep("_web.stable_lights.avg_vis.tif",x)]
  #unzip tar file
  untar(tarfile=f,exdir=here("dload"),files=s)
  #unzip gz
  system(paste0("gunzip -d -f ",file.path("dload",s)))
  #remove unziped gz file
  # not necessary as gunzip authomatically removes the original file, otherwise uncomment the following line:
  # file.remove(file.path("dload",s))
}

Compress file size

The raster tif files are already in Byte format with an average file size of about 700 MB. Using the python based GDAL libraries, the file size can be further reduced to only 23 MB (with an LZW compression) or even 5 MP (with an JPEG compression). JPEG compression formats are usally only used for rasters with a Byte data type (values 0-255), which works in this case. This smoothes a bit the raster values, marginally but noticable. I therefore prefer the LZW compression method. To install GDAL, visit the hompage here.

flist<-list.files("dload","avg_vis\\.tif$")
for (f in flist){
  # f<-flist[1]
  print(f)
  cmd<-paste0(
    c(
      "gdalwarp",
      "-co COMPRESS=LZW",
      here("dload",f),
      here("data",sub(".v4b_web.stable_lights.avg_vis.tif",".tif",f))
    ),
    collapse=" "
  )
  cmd
  system(cmd) ##run the command from within R - it does not take too much time
}

Remove downloads and intermediary products

# file.remove(here("dload"))

Plot output

#nightlight raster
r<-raster("data/F152007.tif")
# plot(r)

#lower resolution raster
ncell(r)

## [1] 725820001

# plot(r)
size.sample<-2*10^5
rd<-sampleRegular(r,size=size.sample, asRaster=T, useGDAL=T)
# table(rd[])
# plot(rd,col="red")

#plot theme
theme.nightlight<-
      ggplot() +
      theme_bw() +
      theme(
        legend.title=element_blank(),
        legend.position="right",
        axis.title=element_blank(),
        axis.text=element_blank()
      ) #+

#plot
p<-
  theme.nightlight +
  geom_tile(data=data.frame(rasterToPoints(rd)),aes(x=x,y=y,fill=F152007),color="black") +
  coord_fixed(
    # xlim = c(96, 140),
    # ylim = c(-10.5, 5.8),
    ratio = 1.3
  ) +
  scale_fill_gradient(low="black", high='white',trans="sqrt")
p


ggsave("nightlight-dmsp-download.jpg",device="jpg",height=15,width=30.8,units="cm")

You may also enjoy