How to extract vineyard land cover from CORINE in R?

26 Views Asked by At

I'm trying to obtain information about land use for vineyards using R by selecting the vineyard category from the TIFF downloaded from CORINE Land Cover.

corine <- rast(paste(dir$land, "U2018_CLC2018_V2020_20u1.tif", sep = "/"))
unique_values <- unique(corine)
print(unique_values) #15 vineyards 
vineyard <- corine==15

I'm not sure about this step. Moreover, I uploaded ESRI shapefile (Simple feature collection with 7901 features and 12 fields Geometry type: MULTIPOLYGON Dimension: XY Bounding box: xmin: 313279.3 ymin: 3933683 xmax: 1312016 ymax: 5220292 Projected CRS: WGS 84 / UTM zone 32N) with administrative boundaries and temperature derived from Copernicus in the NetCDF extension. Becuase I need to obtain temperature for each Italian municipalities weighted for the presence of vineyard to better capture the effect of climate on vineyard.

fname <- paste(dir$shape, "mapEU.shp", sep="/")
shp <- st_read(fname)

fname <- paste(dir$tm, "1979", "tm_19790101.nc", sep="/")
mask <- raster(fname)

# harmonize shape map and temperature layer
shp <-  shp %>% st_transform(crs = mask@crs)
st_crs(shp)

# Change extent in vineyard
e <- extent(mask)
vineyardArea <- crop(vineyard, e)

But I obtained Error:

[crop] extents do not overlap

probably because vineyard is SpatRaster while mask is a RasterLayer. How to fix this issue? Then I'll proceed with something similar ...

# Aggregate vineyard Spatialraster up to 0.1 deg, the resolution of the ERA5 data
vineyardArea <- resample(vineyardArea, mask, method="bilinear")

Thank you in advance for any help!

1

There are 1 best solutions below

0
Robert Hijmans On

I am going to assume that you get

[crop] extents do not overlap

Because the extents do not overlap. Have a look at

fname <- paste(dir$tm, "1979", "tm_19790101.nc", sep="/")
mask <- rast(fname)
ext(mask)
ext(vinyard)
mask
vinyard

Probably one of them has lonlat CRS and the other not. It is hard to help you more without a minimal, self-contained, reproducible example (which is the expectation when asking an R question); or at least some evidence. Based on the above you can perhaps edit your question.

For the polygons you can do

fname <- paste(dir$shape, "mapEU.shp", sep="/")
vct <- vect(fname)
vct <- project(vct, vinyards)