In sustainablefarms/sflddata: Data Extraction Methods for the Linking Data Project

Woody Cover is in Albers Tile format. The coordinate scheme of these files is epsg:3577 (more later in this document on this). The raster data is stored as tiles in folders xmin_ymin, where xmin and ymin are in '000 000s of meters Eastings and Northings. However, the order of the dimensions does not match the netCDF 'CF-1' convention of "T, then Z, then Y, then X" [http://cfconventions.org/Data/cf-conventions/cf-conventions-1.8/cf-conventions.html#dimensions]. I suspect this ordering is due to the WCF data being generated by matlab or fortran. The root of this problem seems to be described here [https://github.com/Unidata/netcdf4-python/issues/337] and it seems likely due to raster assuming a 'CF-1' convention for the netcdf files (in raster's help files), which matches GDAL too [https://gdal.org/drivers/raster/netcdf.html].

This document is to check how to read in tiles correctly. It tests chosen projection information using the very visually distinct location of Sydney.

In past versions of this document, it applied an unofficially modified version of raster that enabled the dimensions to be reordered without saving the data in memory. However, this method broke in the months leading up to August 2021. This document now uses raster::trans and raster::flip to convert the initial raster data to match the actual coordinates. This should still be ok computationally because the tiles are typically less than 50MB each, however things are still much slower.

The document also includes an example method for muffling the discarded datum warning from RGDAL.

library(sflddata)
knitr::opts_chunk$set(echo = TRUE)
invisible(lapply(c("raster", "maptools", "rgdal", "ncdf4", "lubridate"),
                 library, character.only = TRUE))

An arbitrarily chosen tile

nc <- nc_open("[fillmismatch]http://dapds00.nci.org.au/thredds/dodsC/ub8/au/LandCover/DEA_ALC/0_-12/fc_metrics_0_-12_2000.nc")
ncatt_get(nc, "crs")$crs_wkt
nc_close(nc)

This text string cleaned up is:

    spatial_ref=PROJCS[
    "GDA94 / Australian Albers",
    GEOGCS["GDA94",
      DATUM["Geocentric_Datum_of_Australia_1994",
          SPHEROID["GRS1980",
                    6378137,
                    298.257222101,
                    AUTHORITY["EPSG","7019"]],
                    TOWGS84[0,0,0,0,0,0,0],
                    AUTHORITY["EPSG","6283"]],
      PRIMEM["Greenwich",
              0,
              AUTHORITY["EPSG","8901"]],
      UNIT["degree",0.0174532925199433,
           AUTHORITY["EPSG","9122"]],
      AUTHORITY["EPSG","4283"]
    ],
    PROJECTION["Albers_Conic_Equal_Area"],
    PARAMETER["standard_parallel_1",-18],
    PARAMETER["standard_parallel_2",-36],
    PARAMETER["latitude_of_center",0]
    ,PARAMETER["longitude_of_center",132],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","3577"]]

This suggests projection coordinate system is called GDA94 this corresponds to epsg:3577 according to [https://spatialreference.org/ref/epsg/3577/] (which is also given in last property in above list). The datum (spheriod) is GRS1980. The above fits with the projection type (Albers equal area), the ellipsoidal, the first parallel, the second parallel, and longitude in the proj.4 string of epsg:3577 shown below.

CRS("+init=epsg:3577")

However EPSG:3577 has some differences to the CRS that raster extracts from the ncdf files

r <- raster("[fillmismatch]http://dapds00.nci.org.au/thredds/dodsC/ub8/au/LandCover/DEA_ALC/0_-12/fc_metrics_0_-12_2000.nc", varname = "WCF")
crs(r)
rgdal::compare_CRS(CRS("+init=epsg:3577"), crs(r))

All values appear to be the same. RGDAL's compare function says they are equivalent, so the projections are the same up to naming stuff (my guess at what strict means). But the naming stuff appears to be identical too different.

Note that the equivalence seen here is an improvement over last version of this document, wherein raster had some discrepancies in the projection information.

Sydney

The longitude and latitude of Sydney is:

syd_latlong <- SpatialPoints(matrix(c(151.209900,  -33.865143), nrow = 1 ), proj4string = CRS("+proj=longlat +datum=WGS84"))
coordinates(syd_latlong)

Transformed into epsg:3577:

syd_epsg3577 <- spTransform(syd_latlong,
                            CRS("+init=epsg:3577"))
coordinates(syd_epsg3577)

This suggests tile 17_-39. The following reads in this data, converts the raster to the right dimension direction and order. Then plots to check it.

bsyd <- raster("[fillmismatch]http://dapds00.nci.org.au/thredds/dodsC/ub8/au/LandCover/DEA_ALC/17_-39/fc_metrics_17_-39_2000.nc", varname = "WCF")
tbsyd <- raster::t(bsyd)
tfbsyd <- raster::flip(tbsyd, direction = 'x')


plot(tfbsyd)
plot(add = TRUE, syd_epsg3577)

sydzoom <- crop(tfbsyd, extent(buffer(syd_epsg3577, 5000)))
plot(sydzoom)
plot(add = TRUE, syd_epsg3577)

The above tile and point location match the location of Sydney very well. This confirms that the CRS is epsg:3577 and that the transpose and flip of the data has worked as desired.

The tranpose and flip have not altered the objects crs:

compare_CRS(crs(bsyd), crs(tfbsyd))

rgdal::compare_CRS(CRS("+init=epsg:3577"), crs(tfbsyd))

Note that using terra doesn't work. I'm not sure why, as RGDAL can read netCDF.

terra::rast('http://dapds00.nci.org.au/thredds/dodsC/ub8/au/LandCover/DEA_ALC/17_-39/fc_metrics_17_-39_2000.nc')

Without flipping and transposing:

Without flipping and transposing the data the location of Sydney is wrong

plot(bsyd, maxpixel = 10000)
plot(add = TRUE, syd_epsg3577)

The location of Sydney is not even in the extent of this raster, so no point is shown.

Using the bespoke function

This bespoke function using ncdf4 should be working now.

ras <- raster_wcflike("[fillmismatch]http://dapds00.nci.org.au/thredds/dodsC/ub8/au/LandCover/DEA_ALC/17_-39/fc_metrics_17_-39_2000.nc", varname = "WCF")
plot(ras)
plot(add = TRUE, syd_epsg3577)

The year 2020 sometimes has issues

ras2020 <- raster_wcflike("[fillmismatch]http://dapds00.nci.org.au/thredds/dodsC/ub8/au/LandCover/DEA_ALC/17_-39/fc_metrics_17_-39_2020.nc", varname = "WCF")
plot(ras2020)
plot(add = TRUE, syd_epsg3577)

Muffling the Discarded datum warning

Below is a new warning handler defined that invokes a muffled restart if discarded datum geocentric is in a warning. It is invoked globally using globalCallingHandlers.

muffle_datumwarn <- function(w) if (any(grepl("Discarded datum [Uu]nknown", w) | grepl("Discarded datum Geocentric", w))) tryInvokeRestart("muffleWarning")

globalCallingHandlers(warning = muffle_datumwarn)

Tiles Relevant to South-West Slopes

Construct Region Desired

Our region of interest is:

sws_sites <- readRDS("../private/data/clean/sws_sites.rds")
points <- spTransform(sws_sites_2_spdf(sws_sites),
                       CRS("+init=epsg:3577"))
roi <- extent(points)

Find Tiles

Xmins and ymins for a square covering the points:

tilestep <- 100000
lxmin <- floor(roi@xmin / tilestep) * tilestep #lowest xmin
xmins <- seq(lxmin, -1 + ceiling(roi@xmax / tilestep) * tilestep,
    by = tilestep)
lymin <- floor(roi@ymin / tilestep) * tilestep #lowest ymin
ymins <- seq(lymin, -1 + ceiling(roi@ymax / tilestep) * tilestep,
             by = tilestep)

Tiles are then

xmin_v_ymin <- expand.grid(xmin = xmins, ymin = ymins)
tilesidx <- apply(xmin_v_ymin / tilestep, 1, function(x) paste(x, collapse = "_"))
print(tilesidx)

Get files corresponding to tiles

Files for the year 2000 are then:

tilefiles2000 <- paste0("[fillmismatch]http://dapds00.nci.org.au/thredds/dodsC/ub8/au/LandCover/DEA_ALC/",
      tilesidx,
      "/fc_metrics_",
      tilesidx,
      "_2000.nc")

These rasters are:

r.l <- lapply(tilefiles2000,
              function(x){
                         ras <- raster::raster(x, varname = "WCF")
                         firstMax <- raster::cellStats(ras, stat = "max", asSample = FALSE)
                         tras <- raster::t(ras)
                         tfras <- raster::flip(tras, direction = 'x')
                         finalMax <- raster::cellStats(tfras, stat = "max", asSample = FALSE)
                         if (finalMax != firstMax) {
                           warning(sprintf("finalMax (%i) not equal to firstMax (%i) for file %s", finalMax, firstMax, x))}
                         return(tfras)
                        })
rall <- do.call(merge, r.l)

plot(rall, maxpixels = 10000)
plot(add = TRUE, roi)
plot(add = TRUE, points)

Alternately with the bespoke function in sflddata:

r.l <- lapply(tilefiles2000,
              function(x){
                         ras <- raster_wcflike(x, varname = "WCF")
                         return(ras)
                        })
rall <- do.call(merge, r.l)

plot(rall, maxpixels = 10000)
plot(add = TRUE, roi)
plot(add = TRUE, points)

Dimension Orders Check Importing of Landcover Metrics, including Woody Cover

Using the tile over Sydney.

years <- 2018:2019
ncname2017 <- paste0("[fillmismatch]http://dapds00.nci.org.au/thredds/dodsC/ub8/au/LandCover/DEA_ALC/",
      "17_-39",
      "/fc_metrics_",
      "17_-39",
      "_2017.nc")
ncname2018 <- paste0("[fillmismatch]http://dapds00.nci.org.au/thredds/dodsC/ub8/au/LandCover/DEA_ALC/",
      "17_-39",
      "/fc_metrics_",
      "17_-39",
      "_2018.nc")
ncname2019 <- paste0("[fillmismatch]http://dapds00.nci.org.au/thredds/dodsC/ub8/au/LandCover/DEA_ALC/",
      "17_-39",
      "/fc_metrics_",
      "17_-39",
      "_2019.nc")
nc2017 <- nc_open(ncname2017)
nc2018 <- nc_open(ncname2018)
nc2019 <- nc_open(ncname2019)
wcf2017 <- ncvar_get(nc2017, "WCF")
wcf2018 <- ncvar_get(nc2018, "WCF")
wcf2019 <- ncvar_get(nc2019, "WCF")
par(mfrow = c(1, 3))
image(z = t(wcf2017), main = "2017")
image(z = t(wcf2018), main = "2018")
image(z = t(wcf2019), main = "2019")

It looks like the coordinates are no longer flipped for 2019 compared to 2018.

All Combined Using brick_woodycover() Function

syd_epsg3577
syd_surround <- buffer(syd_epsg3577, 500)

sydbrick <- fetch_woody_cover_brick(syd_surround, 2017:2019)
par(mfrow = c(1, 3))
plot(raster(sydbrick, layer = 1))
plot(add = TRUE, syd_epsg3577)
plot(raster(sydbrick, layer = 2))
plot(add = TRUE, syd_epsg3577)
plot(raster(sydbrick, layer = 3))
plot(add = TRUE, syd_epsg3577)
plot(sydbrick)

Another Tile That Has Discrepancies

ras2019 <- raster_wcflike("[fillmismatch]http://dapds00.nci.org.au/thredds/dodsC/ub8/au/LandCover/DEA_ALC/13_-39/fc_metrics_13_-39_2019.nc", varname = "WCF")
plot(ras2019)

ras2020 <- raster_wcflike("[fillmismatch]http://dapds00.nci.org.au/thredds/dodsC/ub8/au/LandCover/DEA_ALC/13_-39/fc_metrics_13_-39_2020.nc", varname = "WCF")
plot(ras2020)

sustainablefarms/sflddata documentation built on April 19, 2022, 11:19 a.m.