chrisvwn / Rnightlights

R package to extract data from satellite nightlights.
GNU General Public License v3.0
47 stars 14 forks source link

Invalid LayerName when ´saving country admLevel structure to CSV #22

Closed felixhaass closed 5 years ago

felixhaass commented 5 years ago

Hi Chris,

first off, thanks for the terrific package. I think it's really useful and I'm sure it'll find wide application. It's definitely useful for my work & I'll recommend it to colleagues.

I have a small problem when processing data for Indonesia, however. I execute the following code, slightly adapted from your Kenya example in the README file:


ctry <- "IDN" #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.Y", 
                                      nlPeriods = nlRange("2012", "2016"), 
                                      nlStats = list("mean",na.rm=TRUE),
                                      ignoreMissing=FALSE)

This results in the following error:

2019-02-12 14:14:38: NlRange autodetected nlType: VIIRS.Y
Processing missing data: IDN:VIIRS.Y:2012:mean, VIIRS.Y:2013:mean, VIIRS.Y:2014:mean, VIIRS.Y:2015:mean, VIIRS.Y:2016:mean. 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-02-12 14:14:38: Downloading country polygons ...
2019-02-12 14:14:38: Downloading polygon: IDN
2019-02-12 14:14:38: Downloading ctry poly: IDN
2019-02-12 14:14:38: Polygon dir for IDN:3.6 already exists
2019-02-12 14:14:39: Saving country admLevel structure to CSV
Error in ctryShpLyrName2Num(ctryCode = ctryCode, layerName = admLevel,  : 
  2019-02-12 14:14:41: Invalid layerName: NA

If I run the Rnightlights:::ctryShpLyrName2Num() function manually, however, it returns a valid integer, e.g.

> Rnightlights:::ctryShpLyrName2Num("IDN", "gadm36_IDN_1")
[1] 0

I'm running the latest GitHub version of the package (the error occurred both in the CRAN and the GitHub version).

Any idea what's going on here?

Thanks so much in advance!

Cheers Felix

chrisvwn commented 5 years ago

Hi @felixhaass,

Thank you for the kind words. I am glad the package is of use to you!

It seems there is a problem with the GADM v3.6 IDN shapefile. It is either missing the country (level 0) layer or it is mislabeled. The GADM site shows there are 5 layers (https://gadm.org/download_country_v3.html) but this is not reflected in the shapefile. Interestingly, the GADM v2.8 works but has only 4 layers including the country layer. If this works for you, you can add the parameter gadmVersion = '2.8' to getCtryNlData to use the older version.

I have written to GADM about this and you can try to do the same via their contact form (http://rasterra.com/contact/gadm_contact_form). Hopefully, they can rectify this and the package will work okay. In future, I plan to add the option to download the RDS maps instead of the shapefiles and add other map sources.

Thanks. I hope we find a solution that works for you.

Chris

felixhaass commented 5 years ago

Hi Chris,

thanks for the super quick reply and your help!

GADM 2.8 works for the time being (which fixes the issue & you can close it), but it's curious that their 3.0 version is missing the IDN_0 layer.

Your planned additions make sense and I'm sure they'll be useful. I'll be watching out for the updates and I'll provide feedback once new versions come online.

Thanks again and cheers Felix

chrisvwn commented 5 years ago

Hi Felix,

You're welcome. GADM has acknowledged receipt of the issue I raised with them concerning this. Hopefully, they will fix it soon. I will update you should they get back to me. Meanwhile, you might want to try v3.6 again maybe in a few weeks.

Thanks

chrisvwn commented 5 years ago

Hi,

So this is a workaround to allow download of the .RDS files where the shapefile zip is not adequate. In particular for the Indonesia (IDN) case.

It is not fully tested so may be buggy. It will be incorporated into a later version of the package.

dnldGADMCtryRDS <- function(ctryCode, gadmVersion = Rnightlights::pkgOptions("gadmVersion"))
{
  if(file.exists(Rnightlights:::getPolyFnameRDS(ctryCode = ctryCode)))
    stop(ctryCode, " ", gadmVersion, " RDS already exists. Please delete it to force")

  if(gadmVersion == "2.8")
    baseUrl <- paste0("https://biogeo.ucdavis.edu/data/gadm", gadmVersion, "/rds/")
  else
    baseUrl <- paste0("https://biogeo.ucdavis.edu/data/gadm", gsub("\\.", "", gadmVersion), "/Rsp/")

  ctryPolyList <- NULL

  res <- 0

  for (idx in 0:4)
  {
    if(gadmVersion == "2.8")
      fName <- paste0(ctryCode, "_adm", idx, ".rds")
    else
      fName <- paste0("gadm", gsub("\\.", "", gadmVersion), "_", ctryCode, "_", idx, "_sp.rds")

    message("Processing ", fName)

    dnldUrl <- paste0(baseUrl, fName)

    tempFname <- file.path(Rnightlights::getNlDir("dirNlTemp"), fName)

    if(!file.exists(tempFname))
      res <- try(download.file(url = dnldUrl, destfile = tempFname, method = "auto"), TRUE)

    if(inherits(res, "try-error"))
      break()

    ctryPoly <- readRDS(file = tempFname)

    lowestIDCol <- paste0("GID_", idx, "_IDX")

    #for GADM 3.6 polygons the attribute col GID_3 contains
    #strings which cannot be used by gdal_rasterize
    #Create an integer col in the shapefile corresponding to the unique GIDs
    if(gadmVersion == "3.6" && is.null(ctryPoly@data[[lowestIDCol]]))
    {
      message(Sys.time(), ": Creating integer zone attribute col for polygon")

      lowestIDColOrig <- gsub("_IDX", "", lowestIDCol)

      ctryPoly@data[,lowestIDCol] <- 1:length(sort(ctryPoly@data[,lowestIDColOrig]))  # Make new attribute
    }

    if(!dir.exists(Rnightlights::getPolyFnamePath(ctryCode = ctryCode, gadmVersion = gadmVersion)))
    {
      message(Sys.time(), ": Writing shapefile layer")
      rgdal::writeOGR(obj = ctryPoly,
                      dsn = Rnightlights::getPolyFnamePath(ctryCode = ctryCode,
                                             gadmVersion = gadmVersion),
                      layer = paste0("gadm_", ctryCode, "_", idx),
                      driver = "ESRI Shapefile",
                      overwrite_layer = T) # Save new version of shapefile
    }

    ctryPolyList <- append(ctryPolyList, ctryPoly)
  }

  message("saving combined RDS")

  if(!file.exists(Rnightlights:::getPolyFnameRDS(ctryCode = ctryCode, gadmVersion = gadmVersion)))
    saveRDS(ctryPolyList, Rnightlights:::getPolyFnameRDS(ctryCode = ctryCode, gadmVersion = gadmVersion))

  if(!file.exists(Rnightlights:::getCtryStructFnamePath(ctryCode, gadmVersion = gadmVersion)))
    Rnightlights:::createCtryStruct(ctryCode = ctryCode, gadmVersion = gadmVersion)
}

You can run it for example for Indonesia version 3.6 as:

dnldGADMCtryRDS(ctryCode = "IDN", gadmVersion = "3.6")
chrisvwn commented 5 years ago

An update. This has been incorporated into the package. Run getCtryNlData with the option gadmPolyType = "spRds" instead of the default gadmPolyType = "shpZip". There are other options ("gpkgZip", "kmlZip", "sfRds") but these are not implemented yet.