Prepared and maintained by: Leandro Parente (OpenGeoHub), Ognjen Antonijevic (GiLAB), Tom Hengl (OpenGeoHub), Luka Glusica (GiLAB), Martin Landa (CVUT Prague)
This tutorial explains (a) how to access Geo-harmonizer geospatial layers, (b) how to use them as a database i.e. to query, subset, and download, (c) how to visualize data in QGIS or similar, and (d) how to add and register new layers. Minimum software requirements to follow this tutorial:
Connected repository: https://gitlab.com/Geo-harmonizer_inea/spatial-layers
Accessing and viewing data
Raster or gridded data
Geo-harmonizer serves a number of raster or gridded layers, usually prepared as Cloud-Optimized GeoTIFFs at 30-m or coarser resolution.
Vector data
Data can be imported into QGIS in a standard way through WFS service using URL:
https://geoserver.opendatascience.eu/geoserver/wfs
The loading might take some time because the dataset is large, therefore it’s best to first zoom to area or feature of interest (for example Switzerland in the image below).
QGIS adding WFS address
QGIS overlaying points and country borders
From then on, you can perform all standard operations as with any other vector data source in QGIS. For example clip data for this country:
QGIS subsetting points for Switzerland
It’s also possible to access this data as GeoJSON directly through browser:
There are various options supported for GetFeatures request: https://docs.geoserver.org/latest/en/user/services/wfs/reference.html#getfeature
And this service can be used directly from GDAL library:
https://gdal.org/drivers/vector/wfs.html
To view vectors in Google Earth, visit the GeoServer page and then select Download → KML → Network link. Once you download the KML you need to zoom-in to an area of interest and then the points with labels should be visible directly in Google Earth.
Visualization of the WFS points in Google Earth
Opening data in QGIS
Geo-harmonizer raster layers can be accessed through HTTP service using QGIS (see also the COG tutorial). These are example of layers:
- landsat_ard_20180625_20180912_p50_RGB: http://s3.eu-central-1.wasabisys.com/eumap/landsat/landsat_ard_20180625_20180912_p50_RGB.tif
- s2l2a_ard_20180625_20180912_p50_RGB: http://s3.eu-central-1.wasabisys.com/eumap/sentinel/s2l2a_ard_20180625_20180912_p50_RGB.tif
Preview of satellite image for EU at 30m imported as an HTTP service to QGIS
Processing data using GDAL and R
To download data from COG repository best use GDAL utilities either gdal_translate or gdalwarp. This allows you to subset, reproject, translate data without a need to download the complete repository that could be tens of Terabytes in size.
Using GDAL utils on COGs
To demonstrate the usage of GDAL commands we will use a bounding box covering Amsterdam, projected in the ETRS89-extended / LAEA Europe coordinate system (EPSG:3035). An easy way to obtain the bounding box coordinates is by using the QGIS bounding box plugin:
- ulx (xmin): 3949019.319534085
- uly (xmax): 3274684.0278522763
- lrx (xmax): 3997655.528969183
- lry (ymin): 3247994.5591947753
The gdal_translate can be used to download data directly from a COG URL. You just need to pass the bounding box information via the “-projwin” parameter, informing the COG url prefixed by “/vsicurl/”:
gdal_translate -co COMPRESS=LZW -projwin 3949019.319534085 3274684.0278522763 3997655.528969183 3247994.5591947753 /vsicurl/http://s3.eu-central-1.wasabisys.com/eumap/landsat/landsat_ard_20180625_20180912_p50_RGB.tif amsterdam.tif |
It is also possible to reproject on-the-fly the raster output, generating an image in, for example, WGS84 projection (EPSG:4326). First it is necessary to create a VRT file using “gdalbuildvrt” command (note the different order for the bounding box coordinates in “-te”):
gdalbuildvrt -te 3949019.319534085 3247994.5591947753 3997655.528969183 3274684.0278522763 amsterdam_4326.vrt /vsicurl/http://s3.eu-central-1.wasabisys.com/eumap/landsat/landsat_ard_20180625_20180912_p50_RGB.tif |
Now just execute a regular gdalwarp informing the desired projection:
gdalwarp -overwrite -co COMPRESS=LZW -t_srs 'EPSG:4326' amsterdam_4326.vrt amsterdam_4326.tif |
The projected cropped image should look like this:
Example GeoTIFF subset
Overlaying a large set of points with COGs
To query the data for any coordinate/point of EU check this Jupyter notebook tutorial.
Clipping maps to a smaller subarea
To access the data for any Geo-harmonizer tile check this Jupyter notebook tutorial.
