In sustainablefarms/linking-data: OBSOLETE: Linking Data Project R Package

Woody Cover is in Albers Tile format. I think 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.

This document is to check the settings required for the package raster to correctly read the tiles. It tests chosen projection information using the very visually distinct location of Sydney.

It applies to an unofficially modified version of raster that contains partial bug fixes. To install this version of raster run: remotes::install_github('https://github.com/kasselhingee/raster', ref = 'ce63b218')") These fixes were not merged into the main raster package because they caused errors in other situations.

knitr::opts_knit$set(root.dir = rprojroot::find_root(rprojroot::has_file("DESCRIPTION")))
devtools::load_all(rprojroot::find_root(rprojroot::has_file("DESCRIPTION")))
knitr::opts_chunk$set(echo = TRUE)
invisible(lapply(c("raster", "maptools", "rgdal", "ncdf4", "lubridate"),
                 library, character.only = TRUE))
if (packageVersion("raster") != "3.0-7") {
    stop(paste("This demonstration uses the 'dims' argument of raster(). This argument requires an unofficial version of the raster package to work properly.",
    "To install this version of raster run:\n remotes::install_github('https://github.com/kasselhingee/raster', ref = 'ce63b218')"))
  }

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) and the ellipsoidal argument in the proj.4 string

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")
proj4string(r)

The values for proj, lat_1, lon_0, lat_0, x_0, y_0 are the same. But the value of lat_2 is different (45.5 vs -36) and the values a and rf do not even appear in CRS("+init=epsg:3577").

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 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", dims = 2:1)
plot(bsyd)
plot(add = TRUE, syd_epsg3577)

sydzoom <- crop(bsyd, 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 discrepancies of lat_2, a and rf can be ignored.

Correcting the CRS of the raster object works:

proj4string(bsyd) <- CRS("+init=epsg:3577")
sydzoom <- crop(bsyd, extent(buffer(syd_epsg3577, 5000)))
plot(sydzoom)
plot(add = TRUE, syd_epsg3577)

In projection as interpreted by raster:

This brings up errors quite quickly: ignore the CRS interpreted by the raster package. For example the location of Sydney in the CRS interpreted by raster is

syd_rasterCRS <- spTransform(syd_latlong,
            CRS("+proj=aea +lat_1=-18 +lon_0=132 +lat_0=0 +x_0=0 +y_0=0 +a=6378137 +rf=298.257222101 +lat_2=45.5"))
coordinates(syd_rasterCRS)

This means xmin_ymin is 21_-37, which doesn't even exist in the catalogue, and certainly doens't match the location of Syndey in the maps.

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){
                      tryCatch(
                        {ras <- raster(x, varname = "WCF", dims = 2:1)
                                     return(ras)
                        }
                      ,
                      warning = function(w) {
                              if (!grep("cannot process these parts of the CRS", as.character(w))){
                                  warning(paste("For", x, w))
                              } else {
                                  suppressWarnings(ras <- raster(x, varname = "WCF", dims = 2:1))
                              }
                      })
                      })
rall <- do.call(merge, r.l)
proj4string(rall) <- CRS("+init=epsg:3577")
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 flipped along the vertical axis in the 2019 compared to 2018. This represents a reversal of the latitude coordinate.

This difference is not reflected in the crs of the netCDF files.

crs2018 <- ncatt_get(nc2018, "crs")
crs2019 <- ncatt_get(nc2019, "crs")
all.equal(crs2018, crs2019)

It also isn't reflected in the attributes of WCF

all.equal(ncatt_get(nc2018, "WCF"), ncatt_get(nc2019, "WCF"))

It also isn't reflected in the 'x' and 'y' coordinate attributes.

all.equal(ncvar_get(nc2018, "x"), ncvar_get(nc2019, "x"))
all.equal(ncvar_get(nc2018, "y"), ncvar_get(nc2019, "y"))

There does appear to be a difference in the "dim" slot of the nc objects. The id component of y$dimvarid is different in nc2019 compared to nc2018 and nc2017.

all.equal(nc2018$dim, nc2019$dim)
all.equal(nc2018$dim, nc2017$dim)

There are other differences too. It is really hard to know which ones are due to the direction of the latitude coordinates.

all.equal(nc2018, nc2019)
all.equal(nc2018, nc2017)

All Combined Using brick_woodycover() Function

syd_epsg3577
syd_surround <- buffer(syd_epsg3577, 500)

sydbrick <- brick_woodycover(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)