Spatial layers for Europe at 30-m available as Cloud Optimized GeoTiffs (corrected)

Thanks, Tom.

How can one use the .qml file you point to in an R workflow to get 
category labels and colors for the COG data you shared originally?

On 16/03/2021 12:24, Tomislav Hengl wrote:

Hi Edzer thanks for spotting this.

Cloud Optimized GeoTiffs seems to be somewhat complicated when it 
comes to adding metadata etc in the tif file directly. For example, 
the COG validator reports:

# rio cogeo validate 
/data/raster/lcv/lcv_landcover.hcl_lucas.corine.rf_p_30m_0..0cm_2019_eumap_epsg3035_v0.1.tif 

he following errors were found:
- The offset of the main IFD should be < 300. It is 1681443308 instead
- The offset of the IFD for overview of index 0 is 750, whereas it 
should be greater than the one of the main image, which is at byte 
1681443308
/data/raster/lcv/lcv_landcover.hcl_lucas.corine.rf_p_30m_0..0cm_2019_eumap_epsg3035_v0.1.tif 
is NOT a valid cloud optimized GeoTIFF

But it seems that the tif fails the check, but works fine on the end.

We are still only testing if this has any effect on speed of access 
etc. The worst case scenario is that users can not access the data 
because of some small glitch in the tif header.

To produce a legend for land cover maps, we recommend using:

https://gitlab.com/geoharmonizer_inea/spatial-layers/-/blob/master/map-style/lcv_landcover.hcl_lucas.corine.rf_p_30m_0..0cm_2000_eumap_epsg3035_v0.1.qml 


You can of course always visually check what you get via:

https://maps.opendatascience.eu

HTH,

T. Hengl
https://opengeohub.org/about

On 3/15/21 4:40 PM, Edzer Pebesma wrote:

Hi Tom, great work!

When reading this either with terra or stars, as in

in.tif = 
"/vsicurl/http://s3.eu-central-1.wasabisys.com/eumap/lcv/lcv_landcover.hcl_lucas.corine.rf_p_30m_0..0cm_2019_eumap_epsg3035_v0.1.tif" 


library(terra)
# terra version 1.1.5
plot(rast(in.tif))

#library(raster)
#plot(raster(in.tif))

library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
plot(read_stars(in.tif))

I see that both plots show a continuous raster, not a categorical. We 
spent quite a bit of time recently in stars & tmap development to get 
handling and plotting of categorical rasters right (including those 
having a color table), but this file doesn't give proper access to 
the categories.

The metadata tags do have them:

stars::gdal_metadata(in.tif)
#? [1] "1=111 - Urban fabric"
#? [2] "10=212 - Permanently irrigated arable land"
#? [3] "11=213 - Rice fields"
#? [4] "12=221 - Vineyards"
#? [5] "13=222 - Fruit trees and berry plantations"
#? [6] "14=223 - Olive groves"
#? [7] "15=231 - Pastures"
#? [8] "16=311 - Broad-leaved forest"
#? [9] "17=312 - Coniferous forest"
# [10] "18=321 - Natural grasslands"
# [11] "19=322 - Moors and heathland"
# [12] "2=122 - Road and rail networks and associated land"
# [13] "20=323 - Sclerophyllous vegetation"
# [14] "21=324 - Transitional woodland-shrub"
# [15] "22=331 - Beaches, dunes, sands"
# [16] "23=332 - Bare rocks"
# [17] "24=333 - Sparsely vegetated areas"
# [18] "25=334 - Burnt areas"
# [19] "26=335 - Glaciers and perpetual snow"
# [20] "27=411 - Inland wetlands"
# [21] "28=421 - Maritime wetlands"
# [22] "29=511 - Water courses"
# [23] "3=123 - Port areas"
# [24] "30=512 - Water bodies"
# [25] "31=521 - Coastal lagoons"
# [26] "32=522 - Estuaries"
# [27] "33=523 - Sea and ocean"
# [28] "4=124 - Airports"
# [29] "5=131 - Mineral extraction sites"
# [30] "6=132 - Dump sites"
# [31] "7=133 - Construction sites"
# [32] "8=141 - Green urban areas"
# [33] "9=211 - Non-irrigated arable land"
# [34] "AREA_OR_POINT=Area"

but GDAL doesn't give access to those in a programmatic way. I've 
tried to add a .aux.xml file with the table, this worked locally (for 
both stars - after using droplevels() - and terra) and might as well 
work over the /vsicurl connection. File attached.

Many regards,


On 02/03/2021 18:30, Tomislav Hengl wrote:

We have mapped land cover classes for the 2000-2019 period for 
continental Europe at 30-m resolution using spatiotemporal Machine
Learning (we used R and python for modeling). Explore the dynamic EU 
landscapes on your palm using the ODS-Europe viewer: 
https://maps.opendatascience.eu

To access almost 10TB of data using R you use the terra or similar
packages e.g.:

R> library(terra)
R> in.tif = 
"/vsicurl/http://s3.eu-central-1.wasabisys.com/eumap/lcv/lcv_landcover.hcl_lucas.corine.rf_p_30m_0..0cm_2019_eumap_epsg3035_v0.1.tif" 

R> tif = rast(in.tif)

?From here you can use any native operation e.g. to crop some polygon
or resample / aggregate values (there is no need to download whole
data sets). A detailed tutorial on how to work with Cloud Optimized
GeoTiffs is available here: 
https://gitlab.com/openlandmap/global-layers/-/blob/master/tutorial/OpenLandMap_COG_tutorial.md. 


Complete list of Cloud Optimized GeoTiffs we produced so far for Europe
is available here: 
https://gitlab.com/geoharmonizer_inea/eumap/-/blob/master/gh_raster_layers.csv 


If not otherwise specified, the data available on this portal is 
licensed under the Open Data Commons Open Database License 
<https://opendatacommons.org/licenses/odbl/> (ODbL) and/or Creative
Commons Attribution-ShareAlike 4.0 
<https://creativecommons.org/licenses/by-sa/4.0/legalcode> and/or 
Creative Commons Attribution 4.0 
<https://creativecommons.org/licenses/by/4.0/legalcode> International
license (CC BY).

Read more in: 
https://opengeohub.medium.com/europe-from-above-space-time-machine-learning-reveals-our-changing-environment-1b05cb7be520 


If you experience any technical problems or if you discover a bug,
please report via: 
https://gitlab.com/geoharmonizer_inea/spatial-layers/-/issues

T. Hengl
https://opengeohub.org/about