Rnightlights: Satellite Nightlight Data Extraction in R

Edit: 2017-11-26 Updated example
Edit 2: 2018-10-22 New version + Updated example

Note: For any comments/assistance with the Rnightlights package it may be easier to post on the package page https://github.com/chrisvwn/Rnightlights/issues

Satellite imagery is increasingly being used in data science to complement or take the place of other data sources. The fact that satellite data are unbiased, have a high spatial resolution, are collected frequently and can be obtained at low cost make them attractive. Nightlights are satellite data that are collected by light sensitive instruments mounted on satellites orbiting around the earth at night to capture the artificial light generated on the surface of the earth i.e. light that is not from the sun or other natural sources.

The nightlights are collected and made available by USA’s NASA (National Aeronautics and Space Administration) and are available at different time resolutions from annual, monthly and most recently daily. Nightlights are available annually from 1992 to 2013 from the DMSP/OLS (Defense Meteorological Satellite Program/Operational Line Scan) system and monthly from April 2012 to present (as of this writing) from the SNPP/VIIRS (Suomi National Polar-Orbiting Partnership) satellite system in the Day/Night Band (DNB). The DMSP/OLS nightlights are presented in one GeoTIFF tile at a resolution of 0.55km x 0.55km at the equator while SNPP/VIIRS comes in 6 GeoTIFF tiles (each approximately 3 GB uncompressed) at a resolution of approximately 0.45kmx0.45km at the equator.

The Rnightlights package aims to provide access to satellite nightlights data for researchers and anyone else without them having to learn all the complexities and details that are required to obtain and process the data. To demonstrate the difference let us look at the process of obtaining monthly nightlights radiance sums at the county level in Kenya for the year 2014.

The manual method
First, the software required. You would need a browser to access and download the tiles, GIS software e.g. QGIS, ArcGIS, etc and statistical software e.g. R.

The procedure to obtain the stats is as follows:
  1. Identify and download administrative level country boundaries for Kenya. The polygon administrative boundaries can be obtained from http://www.gadm.org
  2. Identify the tiles required for Kenya. This can be obtained by getting the geospatial extents of Kenya and finding which VIIRS tiles intersect them.
  3. Locate and download the required tiles. Research would show they can be downloaded from
    https://www.ngdc.noaa.gov/eog/viirs/download_dnb_composites_iframe.html
  4. For each month in 2014:
    • load the tiles and Kenya county polygons
    • Calculate the sum of nightlights at the county level. In QGIS you can use the “zonal statistics” function
    • Save the output as a CSV
  5. Load the CSV in R and perform your analysis
    • kenyaCounties <- read.csv(“kenyaCounties.csv”) 
    • ... etc

Compare this to the Rnightlights procedure:

The Rnightlights method
  1. Open R
  2. Install the Rnightlights package
    • install.packages(“Rnightlights”)
  3. Load the Rnightlights package
    • library(Rnightlights)
  4. Download the country polygons and tiles and calculate the radiance sums at the lowest country level (ward)
    • kenyaWards <- getCtryNlData(ctryCode = "KEN", nlPeriods = nlRange("201401", "201412"), nlType = "VIIRS", nlStats = list("sum",na.rm=TRUE), ignoreMissing="FALSE")
  5. Aggregate the data to the county admin level
    • kenyaWardsMelted <-  reshape2::melt(kenyaWards, value.name="sum")
    • kenyaCounties <- aggregate(kenyaWardsMelted$sum, by=list(kenyaWardsMelted$county), FUN=sum, na.rm=T)
  6. To speed up Rnightlights you can install aria2 which allows for simultaneous downloads of file segments and GDAL which speeds up the processing (raster clipping and zonal statistics) of the tiles. Once installed you can enable them by setting some package methods as below:
    • pkgOptions(downloadMethod = "aria", cropMaskMethod = "gdal", extractMethod = "gdal", deleteTiles = TRUE)

A full code example with visualization using ggplot2 and plotly:

  1. Change the "ctry" variable to run the example on another country
  2. This will not work for countries without admin levels below the country level e.g. ATA (Antarctica)
  3. This calculates total radiances per region and so may be biased by area. Normalize by area to see regions with higher radiances per unit area e.g. to estimate areas with higher economic activity
a) In Rnightlights version 0.2.0 and later calculate directly at the county level with one line:
#install.packages(“Rnightlights”)
#install.packages("lubridate")
#install.packages("reshape2")
#install.packages("ggplot2")
#install.packages("plotly")

library(Rnightlights)
library(lubridate)
library(reshape2)

#(Optional performance enhancement if you have aria2c and gdal installed)
#pkgOptions(downloadMethod = "aria", cropMaskMethod = "gdal", extractMethod = "gdal", deleteTiles = TRUE)

#Optional performance enhancement. If extractMethod="rast" you can specify the number of
#CPU cores to use in parallel
#pkgOptions(extractMethod = "rast", numCores=4)

ctry <- "KEN" #replace to run for any other country

#download and process monthly VIIRS stats at the highest admin level
highestAdmLevelStats <- getCtryNlData(ctryCode = ctry, 
                                     admLevel = "highest",
                                     nlType = "VIIRS.M", 
                                     nlPeriods = nlRange("201401", "201412"), 
                                     nlStats = list("sum",na.rm=TRUE),
                                     ignoreMissing=FALSE)

#Optionally plot the data
library(ggplot2)
library(plotly)

#melt the stats into key-value format for easy multi-line plotting with ggplot2
highestAdmLevelStats <- melt(highestAdmLevelStats,
                              id.vars = grep("NL_", names(highestAdmLevelStats), 
                                             invert=TRUE), 
                              variable.name = "nlPeriod", 
                              value.name = "radiancesum")

#extract date from the NL col names
highestAdmLevelStats$nlPeriod <- substr(highestAdmLevelStats$nlPeriod, 12, 17)

#format period as date
highestAdmLevelStats$nlPeriod <- ymd(paste0(substr(highestAdmLevelStats$nlPeriod, 1,4), 
                                               "-",substr(highestAdmLevelStats$nlPeriod, 5,6), "-01"))

#plot 2nd admin level sums for the year
g <- ggplot(data = highestAdmLevelStats, 
            aes(x=nlPeriod, y=radiancesum, 
                color=highestAdmLevelStats[[2]])) +
  scale_x_date(date_breaks = "1 month", date_labels = "%Y-%m")+
  geom_line()+geom_point() + labs(color = names(highestAdmLevelStats)[2]) + 
  xlab("Month") + 
  ylab("Sum of Radiances") +
  ggtitle(paste0(unique(names(highestAdmLevelStats)[2]), " sum of radiances for ", ctry))

print(g)

#quick conversion to interactive map with plotly
ggplotly(g)

b) Browse the cached data with the internal Shiny app
 
exploreData()
 
One of the benefits of Rnightlights is the exploreData() function which is a Shiny app that allows for some visual analysis of the data that has been downloaded and extracted. This is still in the early stages and hopefully, this will be expanded to pack in more features. A few screenshots of the Shiny interface are below:

Country Sum

Mean Y Log-scale


Mean

View Means

View Sums


View data on the map

View data on the map

Some statistical and clustering views

Some statistical and clustering views

Some statistical and clustering views

View the raw data

Various options available

Comments

  1. Hello, Thank you very much for the example. I tried replicating exactly what you have here but it is not possible. Can you please help me with it and also for the codes in creating a shiny app

    ReplyDelete
    Replies
    1. Hi Ampofo Akwasi,

      Could you send me the error you're encountering? Either copy/paste or send a screenshot. I also realize the code is not very well formatted, so I will try to reformat it in a more friendly way.

      Also, which code are you requesting? If it is the package code, all the code is available in the package so you can study it.

      Delete
    2. Please find below the code and error message I had in running the codes. I would be happy to contact you via email should you provide that to me via akawfo@gmail.com. Thank you.


      The code in full:

      > library(Rnightlights)
      > pkgOptions(downloadMethod = "aria", cropMaskMethod = "gdal", extractMethod = "gdal", deleteTiles = TRUE)
      > kenyaWards <- getCtryNlData(ctryCode = "KEN", nlPeriods = nlRange("201401", "201412"), nlType = "VIIRS", stats = "sum")

      No data found for KEN in 201401:sum, 201402:sum, 201403:sum, 201404:sum, 201405:sum, 201406:sum, 201407:sum, 201408:sum, 201409:sum, 201410:sum, 201411:sum, 201412:sum. Returning NULL.
      Note: Set ignoreMissing=TRUE toreturn only data found or
      ignoreMissing=FALSE to download and extract missing data

      > kenyaWardsMelted <- reshape2::melt(kenyaWards, value.name="sum")
      Error in names(object) <- nm :
      'names' attribute [1] must be the same length as the vector [0]
      > kenyaCounties <- aggregate(kenyaWardsMelted$sum, by=list(kenyaWardsMelted$county), FUN=sum, na.rm=T))
      Error: unexpected ')' in "kenyaCounties <- aggregate(kenyaWardsMelted$sum, by=list(kenyaWardsMelted$county), FUN=sum, na.rm=T))"
      >


      A fuller code example with visualization using ggplot2 and plotly:

      library(Rnightlights)
      > library(reshape2)
      > library (lubridate)

      Attaching package: ‘lubridate’

      The following object is masked from ‘package:base’:

      date

      > ctry <- "KEN"
      > pkgOptions(downloadMethod = "aria", cropMaskMethod = "gdal", extractMethod = "gdal", deleteTiles = TRUE)
      > lowestAdmLevelStats <- getCtryNlData(ctryCode = ctry,
      + nlPeriods = nlRange("201401", "201412"),
      + nlType = "VIIRS",
      + stats = "sum",
      + ignoreMissing=FALSE)

      Processing missing data: KEN in 201401:sum, 201402:sum, 201403:sum, 201404:sum, 201405:sum, 201406:sum, 201407:sum, 201408:sum, 201409:sum, 201410:sum, 201411:sum, 201412:sum. This may take a while.
      Note: Set 'ignoreMissing=TRUE' to return only data found or
      'ignoreMissing=NULL' to return NULL if not all the data is found
      Error in utils::unzip(polyFnameZip, exdir = polyFnamePath) :
      'exdir' does not exist

      > lowestAdmLevelStatsMelted <- melt(lowestAdmLevelStats,
      + id.vars = grep("NL_", names(lowestAdmLevelStats),
      + invert=TRUE),
      + variable.name = "nlPeriod",
      + value.name = "radiancesum")
      Error in melt(lowestAdmLevelStats, id.vars = grep("NL_", names(lowestAdmLevelStats), :
      object 'lowestAdmLevelStats' not found

      > lowestAdmLevelStatsMelted$nlPeriod <- substr(lowestAdmLevelStatsMelted$nlPeriod, 10, 15)
      Error in substr(lowestAdmLevelStatsMelted$nlPeriod, 10, 15) :
      object 'lowestAdmLevelStatsMelted' not found
      >
      > #aggregate the data to 2nd country admin level
      > highestAdmLevelStatsAgg <- setNames(aggregate(lowestAdmLevelStatsMelted$radiancesum,
      + by=list(lowestAdmLevelStatsMelted[[2]],
      + lowestAdmLevelStatsMelted$nlPeriod),
      + FUN=sum, na.rm=T),
      + c(names(lowestAdmLevelStatsMelted)[2], "nlperiod", "sumradiancesums"))

      Error in setNames(aggregate(lowestAdmLevelStatsMelted$radiancesum, by = list(lowestAdmLevelStatsMelted[[2]], :
      object 'lowestAdmLevelStatsMelted' not found
      >

      Delete
    3. Hi,

      Please note the pkgOptions line is only viable if you have aria2 and gdal installed. Please comment it back out. I had also left out an option ignoreMissing=FALSE in the first example. Please copy again and run. I have pasted it below for your convenience:

      install.packages(“Rnightlights”)
      library(Rnightlights)

      #optional performance enhancements

      #pkgOptions(downloadMethod = "aria", cropMaskMethod = "gdal", extractMethod = "gdal", deleteTiles = TRUE, ignoreMissing="FALSE")

      kenyaWards <- getCtryNlData(ctryCode = "KEN", nlPeriods = nlRange("201401", "201412"), nlType = "VIIRS", stats = "sum")
      kenyaWardsMelted <- reshape2::melt(kenyaWards, value.name="sum")
      kenyaCounties <- aggregate(kenyaWardsMelted$sum, by=list(kenyaWardsMelted$county), FUN=sum, na.rm=T))

      Delete
    4. Thank you for getting back to me.

      I get error messages when I try downloading the packages

      > install.packages("gdal")
      --- Please select a CRAN mirror for use in this session ---
      Warning message:
      package ‘gdal’ is not available (for R version 3.4.1)
      > install.packages("aria2")
      Warning message:
      package ‘aria2’ is not available (for R version 3.4.1)
      >


      and even on the latest version of R(version 3.4.3)
      > install.packages("gdal")
      Warning message:
      package ‘gdal’ is not available (for R version 3.4.3)
      > install.packages("aria2")
      Warning message:
      package ‘aria2’ is not available (for R version 3.4.3)

      I would be glad if we can correspond via email. Thank you

      Delete
    5. So aria2 (https://aria2.github.io/) and gdal (http://www.gdal.org/) need to be installed outside of R. But they are not essential. You can skip/comment out the pkgOptions(...) line and run the rest of the code and it should work fine. Do try that and let me know if you still encounter an error.

      Delete
    6. You can create an issue on the package page at https://github.com/chrisvwn/Rnightlights/issues or reach me on chris.njuguna@gmail.com

      Delete
    7. Thank you very much, I really appreciate the assistance.

      Delete
  2. Thank You for the great example. However, I have problem while trying to used the VIIRS data before 2014.
    A call to getAllNlPeriods("VIIRS") shows a valid range which starts from January 2014.

    While trying to use the data for 2013, it alerts an error in nlRange as follows:
    Error in nlRange("201301", "201302") : Invalid start/end nlPeriod

    Is it that I am missing something or there is some way out to work with VIIRS data before 2014.

    Thank You.

    ReplyDelete
  3. Hi Bidur,

    What version of the Rnightlights package are you using? There is one version that was affected by changes in the upstream repository where rasters before 201401 were missing. This was rectified and is updated in later versions.

    Updating to the latest version should fix that. Please try that and let me know if it works.

    Thanks

    ReplyDelete
  4. Thank You very much. I got it working now.

    ReplyDelete
  5. This comment has been removed by the author.

    ReplyDelete
  6. How will I do this if I only limit to regional scale and not for the whole country?

    ReplyDelete
    Replies
    1. Hi. Currently Rnightlights only processes the whole country at the lowest level. The upcoming version will allow you to select the administration level at which you want to perform your calculations. Still, you will not be able to choose a particular region in the selected admin levels.

      Delete
    2. This may be enabled in a future release where you can specify an arbitrary region based either on a subset/aggregation of sub-polygons, by providing an alternative polygon or by specifying a rectangular area of interest. Would this be of interest to you?

      Delete
    3. Oh that's sound so unfortunate. Yea because I'm thinking that I might use this for my analysis assessment for my thesis. :(

      Delete
    4. By the way do you know any alternative where I can get the sum of lights in my region?

      Delete
    5. Why don't you get the sums at the lowest level of your country and then subset to your region and aggregate? Which country and region do you want to calculate sum for?

      Delete
    6. My problem is I don't know how to :( I'm from the Philippines by the way and I'm currently working for the data in my region in Caraga. I just wanna finish my thesis so that I can graduate 😁😁

      Delete
    7. 😁😁
      Alright. So if you look at the example you can see how to run it for the Philippines. Replace the ctryCode so that you have ctryCode <- "PHL".

      I see the admin levels in PHL are "Province_(Lalawigan|Probinsya)", "Municipality_(Bayan|Munisipyo)", "Village_(Barangay)" and Caraga is a collection of 5 provinces, right? If so, you need to first aggregate to the Province level which is level 2 or the top admin level after country. The example already does this for you.

      So after obtaining highestAdmLevelStatsAgg, you can filter the results to keep the Provinces in Caraga. Finally, sum those to get the radiance sum totals for Caraga. I hope this helps.

      Delete
    8. I hope this will work. Thankyou so much for the reply. Can I contact you through email ? So that I can send you the screenshots while I process?

      Delete
    9. You're welcome. You can reach me on chris.njuguna@isharadata.com.

      Delete
  7. The VIIRS DNB images seem to have radiance values as Negative and Positive. When we call the function < getCtryNlData(ctryCode = "SLK", nlPeriods = nlRange("201704", "201705"), nlType = "VIIRS", nlStats = "sum") >, it returns Nightlight Sum for the given ctryCode and nlStats in the given nlPeriods. Do the returned Sum include all the negative and positive values or just the Positive values? Because, the negative values are generally taken as error and discarded during analysis.

    Also, is there any way to specify the threshold value while calling the getCtryNlData() function?

    Thank You.

    ReplyDelete
    Replies
    1. Hi Bidur,

      All negative values are taken to be errors and replaced with NAs so all calculations ignore negative values. This is the generally accepted interpretation of negative radiance values as negative radiance does not have meaning and is usually attributed to instrument errors. getCtryNlData() does not have a threshold parameter. Is this something you require? How would it be used?

      Chris

      Delete
    2. Hi Chris,
      Thank You very much for your clarification. I could not search the correct part of the documentation which states such information.

      Also, I was checking if getCtryNlData() does have a threshold parameter because it seems logical to ignore some values sometimes e.g. very High radiance values.

      The VIIRS DNB monthly data that we get from ngdc.noaa.gov still contain errors ( e.g. lights from fire, boats or aurora). I was just wondering whether Rnightlights library filters such errors.

      Thank You very much for your kind support.

      Delete
    3. Hi,

      No, the package does not remove lights from known light sources. I will look into adding this functionality in a future release. For now, if you need this you can apply this correction using QGIS/ArcGIS to the tile after download and then re-running getCtryNlData. An example of how to do this can be found here: http://darrylmcleod.com/wp-content/uploads/2016/06/Night-Lights-and-ArcGIS-A-Brief-Guide.pdf

      I see what you mean about the threshold. This is not implemented at this point. If you need to, note you can write your own function and pass it into getCtryNlData() as the parameter for nlStats. nlStats can include any function that takes in a numeric vector and returns a scalar e.g.

      mySum <- function(x)
      {
      thresh <- 10000
      x <- x[x < thresh]
      sumx <- sum(x, na.rm=TRUE)
      return(sumx)
      }

      ctryNlData <- getCtryNlData(..., nlStats=mySum)

      If you feel this functionality should be built into the package, do start a discussion/feature request on the package page https://github.com/chrisvwn/Rnightlights/issues

      Chris

      Delete
    4. Hi Chris,

      I will follow your constructive advice.

      Thank You.

      Delete
    5. You're welcome. Just a note. After running getCtryNlData() you can get the path to the folder with the downloaded raster tiles by running:

      getNlDir("dirNlTiles")

      and the cropped country rasters folder by running:

      getNlDir("dirRasterOutput")

      The names should be self-explanatory.

      All the best and do not hesitate to reach out if you need any assistance.

      Delete
    6. Yes Chris,
      Thank You for your kind support.

      Delete
  8. This comment has been removed by the author.

    ReplyDelete
  9. Hi Chris,
    I have some more queries:
    1. I manually calculated radiance statistics like mean, max, min, average for VIIRS 201204 tif image for Macau . Then I compared the same statistics by using RNightlights library:

    lowStats <- getCtryNlData(ctryCode = "MAC", nlPeriods = nlRange("201204", "201204"), nlType = "VIIRS", nlStats =c("min","max","sum","mean","median"), ignoreMissing=FALSE)

    I used the same tif file and same adm shapefile file of Macau for both cases. I got same values for max and min in both cases. However, I got slightly different values for mean and sum for overall Macau. The results for both are in this file https://www.dropbox.com/s/1i7wctl4qe6z4o3/macau2014.xls?dl=0

    2. Also, is it possible to get the statistics like lighted area ( or pixels) , dark area. OR is it possible to calculate the area with specific radiance values.

    Thank You.

    ReplyDelete
    Replies
    1. Code used to manually calculated radiance statistics https://www.dropbox.com/s/68u67ymt1uz03rc/ntlStatistics4Macau.r?dl=0

      Delete
    2. Hi Bidur,

      1. From your code only one issue jumps out i.e. the treatment of negative values.

      values(ntlCountry)[values(ntlCountry)<0] = 0 # Replace Raster Values Less than 0 to 0 ( Thresholding values by ZERO)

      Since these values are unknown I recommend setting them to NA as well since zero is a valid measurement and using it will affect most aggregate calculations. Setting NA should exclude them from all calculations when you use na.rm=TRUE.

      I believe was you wanted to use the threshold for high values? I would do this separately from negative values since they are on opposite ends of your valid data so set it to work on some greater-than-zero value. However, while this *should* resolve for calculations that are affected by NAs, this should not affect sums so there may be something else going on there. I will run your code later tonight and get back to you .

      2. Currently this is not possible but is something I am thinking about. For now all values are provided as a vector so all spatial detail within the sub-polygon is lost. This is because I was initially only interested in simple aggregate stats. I am looking to implement this in a later release so that spatial-aware functions can be applied. This will also be combined with the ability to return structured output instead of single values per sub-polygon as it is now. I will add this as a feature request and tag you just so you can add to/flesh out the idea.

      Thanks!

      Delete
    3. This comment has been removed by the author.

      Delete
    4. Hi Chris,
      Thank You for your kind response.

      As a temporary solution to my problem (2) , is it possible to just count the number of pixels involved in the calculations and write the data to some files by altering the script file ctrynldata.R ( github.com/chrisvwn/Rnightlights/blob/b02c66316e5e635a03f9fd8bccb1b1a50bfbf272/R/ctrynldata.R)
      Once, we know the number of pixels then we can know the area involved in the calculation.

      Thank You.

      Delete
    5. Hi Bidur,

      You're welcome.

      If you simply want to know the count of pixels with a certain value or within a range you can write a function for that. Just filter for the pixel values you want and count. Theoretically, the area of the pixels is constant and you can find the area by multiplying by the pixel area. Unless you want a precise value for the area this should work otherwise you may have to look at a little more complex computation to get the exact area especially the farther from the equator the pixels are since the pixels are less "square" as you move towards the poles.

      Hope this helps.

      Chris

      Delete
    6. Hi Bidur,

      Apologies for taking a while to get back. I have run your code and have seen te differences in results. So I see 2 things happening here. The first is a limitation of the aggregation of aggregated values. Rnightlights currently calculates sub-polygon aggregates at the lowest country level. Getting the mean of this at the country level will not be the same as calculating the mean directly at the country level. The same is true for other aggregate stats except the sum (which is why the example only does the sum). This is addressed in the next version of Rnightlights, out in the next week or so, where one can specify the admin level at which to perform aggregations.

      The second issue is more elusive. Using your calculations I end up with a vector of 180 pixels while Rnightlights gives 168. I think something is happening here as regards pixel inclusion - maybe different strategies are being employed. I am still trying to figure out if this is a difference due to the raster::extract function or some other part of code. Will get back to you on this.

      Delete
    7. Hi Chris,

      Thank You for examining the problem and providing insights.
      I am waiting to see the next version of Rnightlights. I will explore more and communicate with you further.

      Thank You.

      Delete
    8. Hi Bidur,

      I have transferred this discussion to the package issues page for a consolidation of the information here https://github.com/chrisvwn/Rnightlights/issues/5 and it is easier to attach documents.

      Delete
    9. Hi Chris,

      Thank You very much for the consideration.

      Delete
  10. Hi, I tried install "Rnightlights" package in R, but without success. I got the following error message with the installation. Any help?

    ***********
    ERROR: installing vignettes failed
    * removing '\\storage.slu.se/Home$/jute0001/My Documents/R/win-library/3.3/stringr'
    '\\storage.slu.se\Home$\jute0001\My Documents'
    CMD.EXE was started with the above path as the current directory.
    UNC paths are not supported. Defaulting to Windows directory.
    ERROR: dependencies 'gdalUtils', 'stringr' are not available for package 'Rnightlights'
    * removing '\\storage.slu.se/Home$/jute0001/My Documents/R/win-library/3.3/Rnightlights'

    The downloaded source packages are in
    ‘C:\Users\jute0001\AppData\Local\Temp\Rtmpod5v3Q\downloaded_packages’
    Warning messages:
    1: running command '"C:/PROGRA~1/R/R-33~1.2/bin/x64/R" CMD INSTALL -l "\\storage.slu.se\Home$\jute0001\My Documents\R\win-library\3.3" C:\Users\jute0001\AppData\Local\Temp\Rtmpod5v3Q/downloaded_packages/pillar_1.2.3.tar.gz' had status 1
    2: In install.packages(NULL, .libPaths()[1L], dependencies = NA, type = type) :
    installation of package ‘pillar’ had non-zero exit status
    3: running command '"C:/PROGRA~1/R/R-33~1.2/bin/x64/R" CMD INSTALL -l "\\storage.slu.se\Home$\jute0001\My Documents\R\win-library\3.3" C:\Users\jute0001\AppData\Local\Temp\Rtmpod5v3Q/downloaded_packages/gdalUtils_2.0.1.14.tar.gz' had status 1
    4: In install.packages(NULL, .libPaths()[1L], dependencies = NA, type = type) :
    installation of package ‘gdalUtils’ had non-zero exit status
    5: running command '"C:/PROGRA~1/R/R-33~1.2/bin/x64/R" CMD INSTALL -l "\\storage.slu.se\Home$\jute0001\My Documents\R\win-library\3.3" C:\Users\jute0001\AppData\Local\Temp\Rtmpod5v3Q/downloaded_packages/stringr_1.3.1.tar.gz' had status 1
    6: In install.packages(NULL, .libPaths()[1L], dependencies = NA, type = type) :
    installation of package ‘stringr’ had non-zero exit status
    7: running command '"C:/PROGRA~1/R/R-33~1.2/bin/x64/R" CMD INSTALL -l "\\storage.slu.se\Home$\jute0001\My Documents\R\win-library\3.3" C:\Users\jute0001\AppData\Local\Temp\Rtmpod5v3Q/downloaded_packages/Rnightlights_0.2.1.tar.gz' had status 1
    8: In install.packages(NULL, .libPaths()[1L], dependencies = NA, type = type) :
    installation of package ‘Rnightlights’ had non-zero exit status

    ReplyDelete
    Replies
    1. Hi Jumens,

      Apologies. I am just seeing this comment. I hope you were able to get the package installed. If not, please get back to me. It may be easier to raise the issue on the package page at https://github.com/chrisvwn/Rnightlights/issues

      Regards,
      Chris

      Delete
  11. Hi Chris,
    Thanks for the detailed guidance ! I am getting this error "Invalid admLevels: IND:IND_adm0
    Error in getCtryNlData(ctryCode = "IND", admLevel = "IND_adm0", nlType = "VIIRS.M", :
    Invalid admLevels detected" while extracting data for India. I tried to search admLevels using "searchAdmLevel" but could not trace it.

    Can you please help me !

    ReplyDelete
    Replies
    1. Hi Nabin,

      Apologies. I am just seeing this comment. I hope you were able to solve this error. If not, please get back to me. It may be easier to raise the issue on the package page at https://github.com/chrisvwn/Rnightlights/issues

      Regards,
      Chris

      Delete
  12. Hi,

    Thank you for the awesome package.
    I was able to get country raster and split the area into administrative units, but I never get the corrects nlStats. I tried running the code on three different countries, using both VAARS and OLS data in random years between 2002 and 2016, but all the numbers always come out to NA. Do you know what could be the problem here? The tiff files look fine, and I was able open them in QGIS to see the range of values (which was between -0.76 and 192 for one of them iirc, don't know if that means anything), but I'm new to this kind of software so I don't know how to apply GDAL shapes and compute everything manually.

    Here's some sample code

    > highestAdmLevelStats <- getCtryNlData(ctryCode = ctry,
    + admLevel = "highest",
    + nlType = "VIIRS.M",
    + nlPeriods = nlRange("201401", "201401")
    + nlStats = "sum",
    + ignoreMissing=FALSE)

    And the resulting DF

    highestAdmLevelStats
    country county area_sq_km NL_VIIRS.M_201401_SUM
    1 KEN Baringo 10842.4553 NA
    2 KEN Bomet 2378.4699 NA
    3 KEN Bungoma 3573.1478 NA

    I would really appreciate it if you could point me in the right direction.

    ReplyDelete
    Replies
    1. Hi,

      Apologies for this. I need to update the documentation and blog. The latest version of Rnightlights has a new format for the nlStats parameter. Each nlStat is a list consisting of the name of the function and the parameters to pass to the function. Note that na.rm=TRUE is not longer assumed. So you can try nlStats=list("sum",na.rm=TRUE).

      Delete
    2. Thank you so much, works like a charm!
      One more quick question about how radiance is converted into the tiff. It's much brighter in some years and darker in others. Opening the image in QGIS yields an almost all-black picture for some countries, unless the grayscale is adjusted (the vast majority of the values are muh smaller than the maximum on the picture). Would adjusting the grayscale range to the same values for all rasters give me a better graphic representation of the changes over a certain area?

      Once again, thanks for the library & the support

      Delete
    3. Hi,

      You're welcome. I am glad you are finding the package of use.

      I am not sure what the best tool for this would be. For visual representation, you may try playing around with the ranges in the QGIS style under raster properties. Alternatively, the raster calculator might create a bunch of rasters with the ranges you require. You might also try doing this in R using the raster package.

      I hope this helps.

      Delete
  13. Hi Ishara,

    Thank you for your post! It is really helpful.

    I am trying to replicate your code, and I have a following issue from R:

    > highestAdmLevelStats <- getCtryNlData(ctryCode = ctry,admLevel = "highest",nlType = "VIIRS.M", nlPeriods = nlRange("201401", "201402"),nlStats = list("sum",na.rm=TRUE), ignoreMissing=FALSE)

    (...)

    Error in .local(.Object, ...) :
    TIFFReadDirectory:Failed to read directory at offset 2073600008

    In addition: There were 50 or more warnings (use warnings() to see the first 50)
    Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer", :
    Cannot create a RasterLayer object from this file.

    Do you have any ideas how to fix it? Again thank you for your help!

    Regards,
    Anh

    ReplyDelete
    Replies
    1. Dear Anh,

      Could you post the full log. I think something might have gone wrong before this error occurred. If possible also you can open an issue at https://github.com/chrisvwn/Rnightlights/issues for easier issue handling.

      Regards,

      Chris

      Delete
    2. Hi Chris,

      Thank you for your reply! I have just found out that the issue occurred because my disk space was full. But the process took a large amount of space, I estimated that each month (from 201401 to 201412) might take approximate 4,5 Gb.

      Now I working on an analysis which requires the light data at the lowest administrative level across countries. Do you think if I want to export the data from R to excel (not the graph and other related files), the process will be quicker and take lesser space? If yes, I would much appreciate if you could show me an example of code using R.

      Thank you for your time!

      Regards,
      Anh

      Delete
    3. Hi Anh,

      4.5GB per month sounds a bit high. An option though is to delete tiles after use. You can run pkgOptions(deleteNlTiles=TRUE) before running getCtryNlData to dispose of the large tiles once processing is complete.

      It is not difficult to export the data to Excel if you want to since it is just a data.frame saved as csv at the moment. You can also save directly to csv or xls(x) if you prefer using the openxlsx package or other similar package.

      To get the stats calculated at the lowest level set admLevel="lowest" either in getCtryNlData or processNlData.

      Regards,

      Chris

      Regards,

      Chris

      Delete
    4. Hi Christ,

      Thanks again for your time!

      As I am beginner in R, so sorry for asking "naive" questions. Is it possible if you could have a look at the code below and let me know if I am on the right track?

      >library(Rnightlights)
      # Selection 0
      > library(lubridate)
      > library(reshape2)
      > pkgOptions(deleteNlTiles=TRUE)
      # Above is the code you suggested last time. Unfortunately, it did not work. Here is the error:
      Warning message:
      In RNIGHTLIGHTSOPTIONS(...) :
      Ignoring options not defined in manager: deleteNlTiles

      > ctry <- "KEN"
      > highestAdmLevelStats <- getCtryNlData(ctryCode = ctry,
      + admLevel = "highest",
      + nlType = "VIIRS.M",
      + nlPeriods = nlRange("201401", "201401"),
      + nlStats = list("sum",na.rm=TRUE),
      + ignoreMissing=FALSE)
      # Please note that I only take an example of 201401. I think that will save more space and time.

      Regards,
      Anh

      Delete
    5. Hi Anh,

      I made a mistake with the tiles option. It should be pkgOptions(deleteTiles=TRUE).

      The code otherwise looks okay. Could you post the full output of the getCtryNlData?

      Regards,

      Chris

      Delete
  14. Hi Chris,

    The new code pkgOptions(deleteTiles=TRUE) works well. Now I can replicate your code. The tiles are omitted which saves a lot of space in the hard disk. It took 1 hour to retrieve the data from 201401 to 201402 in Kenya, but that's not a problem.

    About other countries, may you please let me know how I can find the country code (such as 'KEN' for Kenya)? Also, I was wondering about the measure of light data. From my understanding, the light index ranges from 0 to 63, but I am not sure if VIIRS uses the same scale. I notice that for the case of Kenya, there are many tiny numbers (for the mean) and large numbers (for the sum).

    Thank you so much!

    Regards,
    Anh

    ReplyDelete
    Replies
    1. Hi Anh,

      I'm glad this worked for you. 1 hour is too long for 2 months or did you mean 201401 to 201412? Unless the download speed is slow.

      You can find country codes using the function searchCountry or ctryNameToCode.

      VIIRS uses a different scale from OLS. It theoretically has a max of 65536. They are actually much smaller values but have been multiplied to make it easier for calculations. For info on viirs nightlights see: https://maps.ngdc.noaa.gov/viewers/VIIRS_DNB_monthly/

      Regards,

      Chris

      Delete
  15. Hi Chris,
    I want to get data on nightlight for 201205 to 201304 for Indian districts. I am not being able to understand from this portion -
    the stats into keyvalue format for easy multi-line plotting with ggplot2
    highestAdmLevelStats <- melt(highestAdmLevelStats,
    id.vars = grep("NL_", names(highestAdmLevelStats),
    invert=TRUE),
    variable.name = "nlPeriod",
    value.name = "radiancesum").
    Can you please explain to me what these mean.
    Regards

    ReplyDelete
    Replies
    1. Hi Soma,

      Essentially we are reshaping the data.frame into a form we can pass into ggplot2 for easy plotting of multiple lines at once. So it changes from a 1 row = 1 admLevel + multiple stats to 1 row = 1 admLevel+ 1 nlStat. The easiest way to understand this is to run it and compare the data.frame before and after.

      e.g. this row:
      country county area_sq_km NL_OLS.Y_1992_SUM NL_OLS.Y_1993_SUM NL_OLS.Y_1994_SUM NL_OLS.Y_1995_SUM
      1 KEN Baringo 10842.4553 29 63 63 224

      is converted into:
      country county area_sq_km nlPeriod radiancesum
      1 KEN Baringo 10842.4553 NL_OLS.Y_1992_SUM 29
      48 KEN Baringo 10842.4553 NL_OLS.Y_1993_SUM 63
      95 KEN Baringo 10842.4553 NL_OLS.Y_1994_SUM 63
      142 KEN Baringo 10842.4553 NL_OLS.Y_1995_SUM 224

      So all the values for one nlPeriod will become a "series" as it were (notice the row numbers).

      Hope this helps.

      Regards,

      Chris

      Delete
    2. Hi Chris,
      Thank you for your information. I started running the program as is directed here. However I am facing some problem while running the codes. I am attaching the entire code and the response from R-
      > highestAdmLevelStats <- getCtryNlData(ctryCode = ctry,
      + admLevel = "highest",
      + nlType = "VIIRS.M",
      + nlPeriods = nlRange("201205", "201304"),
      + nlStats = list("sum",na.rm=TRUE),
      + ignoreMissing=FALSE)
      2018-11-02 02:09:27: NlRange autodetected nlType: VIIRS.M
      Processing missing data: IND:VIIRS.M:201205:sum, VIIRS.M:201206:sum, VIIRS.M:201207:sum, VIIRS.M:201208:sum, VIIRS.M:201209:sum, VIIRS.M:201210:sum, VIIRS.M:201211:sum, VIIRS.M:201212:sum, VIIRS.M:201301:sum, VIIRS.M:201302:sum, VIIRS.M:201303:sum, VIIRS.M:201304:sum. This may take a while.
      Note: Set 'ignoreMissing=TRUE' to return only data found or
      'ignoreMissing=NULL' to return NULL if not all the data is found
      2018-11-02 02:09:30: Downloading country polygons ...
      2018-11-02 02:09:30: Downloading polygon: IND
      2018-11-02 02:09:30: Downloading ctry poly: IND
      2018-11-02 02:09:30: Polygon dir for IND:3.6 already exists
      2018-11-02 02:09:30: Downloading country polygons ... DONE
      2018-11-02 02:09:31: **** PROCESSING nlType:VIIRS.M nlPeriod:201205****
      2018-11-02 02:09:31: Checking tiles required for VIIRS.M 201205
      2018-11-02 02:09:31: IND: Stats missing. Adding tiles
      2018-11-02 02:09:331 Required tiles: 75N060E
      2018-11-02 02:09:33: Downloading tile: 2012053
      2018-11-02 02:09:33: File exists, set Overwrite = TRUE to overwrite
      2018-11-02 02:09:33: Extracting C:/Users/Soma Bhaduri/Documents/.Rnightlights/tiles/NL_TILE_VIIRS.M_201205_75N060E.tgz
      2018-11-02 02:09:33: TIF file found
      2018-11-02 02:09:33: **processNLCountry: IND gadm36_IND_1 VIIRS.M 201205
      2018-11-02 02:09:33: Check for existing data file
      2018-11-02 02:09:33: Data file not found. Creating ...
      2018-11-02 02:09:39: Data file not found. Creating ... DONE
      2018-11-02 02:09:39: Load polygon layer for crop
      2018-11-02 02:09:40: Begin processing 201205
      2018-11-02 02:09:40: Reading in the raster tiles
      2018-11-02 02:09:42: Cropping the raster tiles
      Error: cannot allocate vector of size 369.2 Mb


      Can you help me regarding what to do next?

      Regards

      Delete
    3. Hi Soma,

      I suspect the tile did not complete downloading in a previous run. Try deleting C:/Users/Soma Bhaduri/Documents/.Rnightlights/tiles/NL_TILE_VIIRS.M_201205_75N060E.tgz and try again.

      Delete
  16. Hi Chris,
    I was trying to extract the annual VIIRS composite date for 2016 and I have the code pasted below:
    ctry <- "THA"

    highestAdmLevelStats <- getCtryNlData(ctryCode = ctry,
    admLevel = "highest",
    nlType = "VIIRS.Y",
    nlPeriods = "2016",
    ignoreMissing=FALSE)

    However I get the error:
    Error in utils::download.file(url = ntLtsIndexUrlVIIRS, destfile = ntLtsPageLocalName, :
    cannot open URL 'https://www.ngdc.noaa.gov/eog/viirs/download_dnb_composites_iframe.html'
    In addition: Warning message:
    In utils::download.file(url = ntLtsIndexUrlVIIRS, destfile = ntLtsPageLocalName, :
    InternetOpenUrl failed: 'A connection with the server could not be established'

    Do I have to specify another download URL?

    ReplyDelete
    Replies
    1. Hi Ugyen,

      You do not need to specify another download URL. For some reason the package cannot access the download page. Can you access the url https://www.ngdc.noaa.gov/eog/viirs/download_dnb_composites_iframe.html in your browser?

      Regards,

      Chris

      Delete
  17. Hi Chris,
    I want to get data on nightlight for 201205 to 201304 for Indian districts. But while I am running the program as is directed here I am facing some problem while running the codes. I am attaching the response from R- Error in melt(highestAdmLevelStats, id.vars = grep("NL_", names(highestAdmLevelStats), :
    object 'highestAdmLevelStats' not found.

    Can you help me regarding how to resolve the issue?

    ReplyDelete
  18. Hi Ishara,

    I am trying to use your code for Turkey. When I run "highestAdmLevelStats" part I am getting an error as following:

    2019-06-13 09:59:46: NlRange autodetected nlType: VIIRS.M
    Processing missing data: TUR:VIIRS.M:201501:sum, VIIRS.M:201502:sum, VIIRS.M:201503:sum, VIIRS.M:201504:sum, VIIRS.M:201505:sum, VIIRS.M:201506:sum, VIIRS.M:201507:sum, VIIRS.M:201508:sum, VIIRS.M:201509:sum, VIIRS.M:201510:sum, VIIRS.M:201511:sum, VIIRS.M:201512:sum. This may take a while.
    Note: Set 'ignoreMissing=TRUE' to return only data found or
    'ignoreMissing=NULL' to return NULL if not all the data is found
    2019-06-13 09:59:49: Downloading country polygons ...
    2019-06-13 09:59:49: Downloading polygon: TUR
    2019-06-13 09:59:50: Downloading ctry poly: TUR
    2019-06-13 09:59:50: Polygon dir for TUR:3.6 already exists
    2019-06-13 09:59:50: Downloading country polygons ... DONE
    2019-06-13 09:59:53: **** PROCESSING nlType:VIIRS.M nlPeriod:201501****
    2019-06-13 09:59:53: Checking tiles required for VIIRS.M 201501
    2019-06-13 09:59:53: TUR: Stats missing. Adding tiles
    2019-06-13 09:59:551 Required tiles: 75N060W
    2019-06-13 09:59:55: Downloading tile: 2015012
    2019-06-13 09:59:56: An error occurred downloading
    2019-06-13 09:59:56: **processNLCountry: TUR gadm36_TUR_1 VIIRS.M 201501
    2019-06-13 09:59:56: Check for existing data file
    2019-06-13 09:59:56: Data file not found. Creating ...
    2019-06-13 10:00:01: Data file not found. Creating ... DONE
    2019-06-13 10:00:01: Load polygon layer for crop
    2019-06-13 10:00:02: Begin processing 201501
    2019-06-13 10:00:02: Reading in the raster tiles
    Error in .local(.Object, ...) :

    In addition: Warning message:
    running command 'aria2c -c -x2 --show-console-readout=false --summary-interval=10 https://data.ngdc.noaa.gov/instruments/remote-sensing/passive/spectrometers-radiometers/imaging/viirs/dnb_composites/v10//201501/vcmcfg/SVDNB_npp_20150101-20150131_75N060W_vcmcfg_v10_c201505111709.tgz -d C:\Users\KULLAN~1\AppData\Local\Temp\RtmpSKXAGW/.Rnightlights/tiles -o
    NL_TILE_VIIRS.M_201501_75N060W.tgz' had status 127
    Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer", :
    Cannot create a RasterLayer object from this file. (file does not exist)

    Can you help me how to solve this problem?

    Regards

    Gızem

    ReplyDelete
  19. Hi Gizem,

    The only thing I can tell is that the download failed. Could you try running the download manually from the commandline i.e.

    aria2c -c -x2 --show-console-readout=false --summary-interval=10 https://data.ngdc.noaa.gov/instruments/remote-sensing/passive/spectrometers-radiometers/imaging/viirs/dnb_composites/v10//201501/vcmcfg/SVDNB_npp_20150101-20150131_75N060W_vcmcfg_v10_c201505111709.tgz -d C:\Users\KULLAN~1\AppData\Local\Temp\RtmpSKXAGW/.Rnightlights/tiles -o NL_TILE_VIIRS.M_201501_75N060W.tgz

    And post the output from that command.

    ReplyDelete
  20. Hello Ishara,

    Thanks for the awesome effort!
    Sorry to bother you, but I'm trying this code:
    highestAdmLevelStats <- getCtryNlData(ctryCode = ctry,
    admLevel = "highest",
    nlType = "VIIRS.M",
    nlPeriods = nlRange("200501", "200512", "VIIRS.M"),
    nlStats = list("sum",na.rm=TRUE),
    ignoreMissing=FALSE)

    and facing this error when trying for Lebanon:

    "
    2019-07-26 16:49:02Invalid nlPeriods:: VIIRS.M:200501,
    Error in nlRange("200501", "200512", "VIIRS.M") :
    2019-07-26 16:49:02: Invalid nlPeriod detected for nlType VIIRS.M
    "
    could it be because there was no data for the period that I'm looking for in 2005? is there a way to check?

    Thanks a lot!
    Marie

    ReplyDelete
    Replies
    1. Hi Marie,

      Yes the problem is that there was no VIIRS data during that time. OLS is available annually from 1992-2012 while VIIRS.M is available monthly from 201204 to present.

      Delete
    2. Thanks a lot Chris! Sorry for the silly question, I was so convinced it was an error in the code I overlooked the data type. Extracting data right now. All my best, Marie

      Delete

Post a Comment