0

I have a netCDF file ("SSP119.nc") containing information on land-use change from 2015 to 2100, inclusive (source: https://luh.umd.edu/data.shtml). I also have a shapefile containing polygons of global Key Biodiversity Areas (KBA). I want to know how much land area overlaps with the KBA polygons for certain land-use types in certain years (yep, it's complicated... but I hope my code makes it easier to understand).

As an aside, the land-use data contains information for 14 land-use types. The resolution is at 0.25 x 0.25 degree grid-cells. Each grid-cell contains the fraction of each land-use type within it, e.g., if one land-use type covers 60% of the grid-cell, its value will be 0.6.

Here is my code, so far:

SSP119_r_primf <- stack("SSP119.nc", varname = "primf") # this creates a raster stack of the
# land-use "primary forested area" from the netCDF file

SSP119_r_primf
`class      : RasterStack 
dimensions : 720, 1440, 1036800, 86  (nrow, ncol, ncell, nlayers)
resolution : 0.25, 0.25  (x, y)
extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
names      : X0, X1, X2, X3, X4, X5, X6, X7, X8, X9, X10, X11, X12, X13, X14, ... 
years since 2015-01-01 0:0:0: 0 - 85 (range)
`
proj4string(SSP119_r_primf) = CRS("+init=EPSG:4326")

plot(SSP119_r_primf) # will produce 86 plots of primary forested areas,
# annually between 2015 and 2100

primf_2100 <- subset(SSP119_r_primf, 86) # creates a subset of SSP119_r_primf
# specifically for the year 2100

plot(primf_2100) # to give you a visual idea of the data
[Primary forested land in year 2100](https://i.stack.imgur.com/TsVo7.png)

KBA <- readOGR("KBAsGlobal_2022_September_02_POL.shp") # to read in the KBA
# shapefile polygon data

Any ideas on the next steps? Thanks a million, appreciate any help!

Conor Ó Beoláin
  • 1
  • 1
  • 1
    My first recommendation is to become familiar with the stars package. https://cran.r-project.org/web/packages/stars/vignettes/stars1.html https://r-spatial.github.io/stars/articles/ Also become familiar with the leaflet package. Leaflet is great for interactive plots. It works with shapefiles and netCDF files. https://rpubs.com/charlieb/387239 Finally become familiar with the sf package. sf has many spatial vector data functions. – YLC Jan 30 '23 at 20:10
  • The basic idea is to plot a set of KBA ploygons and overlay a set of points from you netCDF file. The stars package allows you to read and process netCDF files The stars package has functions to build a square around a point. Use the sf package to convert points and squares into suitable geometries. You can set the color of each point(or square) given a value of your choice. Use leaflet to overlay squares on top of your KBA map. – YLC Jan 30 '23 at 20:11

1 Answers1

0

It would have been useful if you had been specific about the file you use. I am using another file from the same webpage (this one has added tree cover fraction over time).

url <- "https://luh.umd.edu/LUH2/LUH2_v2f/added_tree_cover/added_tree_cover_input4MIPs_landState_ScenarioMIP_UofMD-IMAGE-ssp119-2-1-f_gn_2015-2100.nc"
f <- basename(url)
if (!file.exists(f)) download.file(url, f, mode="wb")

You are using obsolete R packages ("raster", "sp", "rgdal"). Let's use "terra" instead.

library(terra)
x <- rast(f)
x
#class       : SpatRaster 
#dimensions  : 720, 1440, 85  (nrow, ncol, nlyr)
#resolution  : 0.25, 0.25  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 
#source      : added_tree_cover_input4MIPs_landState_ScenarioMIP_UofMD-IMAGE-ssp119-2-1-f_gn_2015-2100.nc:added_tree_cover 
#varname     : added_tree_cover (area_fraction) 
#names       : added~ver_1, added~ver_2, added~ver_3, added~ver_4, added~ver_5, added~ver_6, ... 
#unit        :           1,           1,           1,           1,           1,           1, ... 
#time (years): 2015 to 2099

You could read your polygon data with

v <- vect("KBAsGlobal_2022_September_02_POL.shp")

But since we do not have that file, I will use polygons for the countries of the world

library(geodata)
v <- geodata::world(path=".")

# extract raster values, and compute the mean for each polygon 

e <- extract(x, v, mean, ID=FALSE)
e <- cbind(country=v$COUNTRY, e)

e[1:5, 1:5]
#      country added_tree_cover_1 added_tree_cover_2 added_tree_cover_3 added_tree_cover_4
#1       Aruba       0.0053531881       0.000000e+00       0.000000e+00       0.000000e+00
#2 Afghanistan       0.0000498043       4.980157e-05       4.979897e-05       4.979643e-05
#3      Angola       0.0011926381       1.192639e-03       1.199086e-03       1.199087e-03
#4    Anguilla       0.0000000000       0.000000e+00       0.000000e+00       0.000000e+00
#5       Åland       0.0000000000       0.000000e+00       0.000000e+00       0.000000e+00

If your polygons are very small relative to the grid cell size, that is when you only have a few cells per polygon, you may want to use argument exact=TRUE to better account for cells that are partly inside of a polygon.

For future reference, to ask a question like this it is better to use data generated by code, or with files that ship with R. See the examples in the R help file for how this is done (e.g. ?rast, ?vect, ?extract, ?zonal.

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63