Merge pull request #7 from geographika/wcs

Add WCS tutorial

Author: geographika

workshop/content/docs/outputs/wcs.md

@@ -1,50 +1,183 @@
-# Web Coverage Services
+# Web Coverage Services (WCS)
 
-TODO
+## Overview
 
-https://mapserver.org/ogc/wcs_server.html
+This tutorial demonstrates the [WCS Server](https://mapserver.org/ogc/wcs_server.html) capabilities of MapServer.
+We'll be using [WCS 2.0](https://mapserver.org/ogc/wcs_server.html#wcs-2-0) for this tutorial, and will serve a
+Cloud-Optimized GeoTIFF (COG) from the Estonian Land Board as the source dataset. The dataset is a Digital Terrain Model (DTM) with a 1 m resolution.
 
-https://mapserver.org/ogc/wcs_format.html
+<div class="map">
+  <iframe src="https://mapserver.github.io/getting-started-with-mapserver-demo/wcs.html"></iframe>
+</div>
 
-We'll be using [WCS 2.0](https://mapserver.org/ogc/wcs_server.html#wcs-2-0) for this tutorial.
+## WCS Requests
 
-# GetCapabilities
-http://localhost:7000/?map=/etc/mapserver/wcs.map&SERVICE=WCS&REQUEST=GetCapabilities
+Some sample MapServer requests for testing the WCS service are listed below. You can test these in your browser.
 
-# DescribeCoverage 2.0
-http://localhost:7000/?map=/etc/mapserver/wcs.map&SERVICE=WCS&VERSION=2.0.1&REQUEST=DescribeCoverage&COVERAGEID=dtm
+- [GetCapabilities](http://localhost:7000/?map=/etc/mapserver/wcs.map&SERVICE=WCS&REQUEST=GetCapabilities&VERSION=2.0.1)
+- [DescribeCoverage 2.0](http://localhost:7000/?map=/etc/mapserver/wcs.map&SERVICE=WCS&VERSION=2.0.1&REQUEST=DescribeCoverage&COVERAGEID=dtm)
+- [GetCoverage 2.0 image/tiff full](http://localhost:7000/?map=/etc/mapserver/wcs.map&SERVICE=WCS&VERSION=2.0.1&REQUEST=GetCoverage&COVERAGEID=dtm&FORMAT=image/tiff)
 
-# GetCoverage 2.0 image/tiff full
-http://localhost:7000/?map=/etc/mapserver/wcs.map&SERVICE=WCS&VERSION=2.0.1&REQUEST=GetCoverage&COVERAGEID=dtm&FORMAT=image/tiff
+You can also connect to the MapServer Docker container and use `mapserv` to test the requests from the command line.
 
 ```bash
-gdalinfo /etc/mapserver/data/raster/54752_dtm_1m.tif
+docker exec -it mapserver bash
+mapserv -nh "QUERY_STRING=map=/etc/mapserver/wcs.map&SERVICE=WCS&REQUEST=GetCapabilities&VERSION=2.0.1"
 ```
 
-NoData Value=-9999
-Size is 5002, 5002
+## Source Dataset
 
-curl "http://localhost:7000/?map=/etc/mapserver/wcs.map&SERVICE=WCS&REQUEST=GetCapabilities"
+Let's get some information about the source dataset using the GDAL CLI command[gdal raster info](https://gdal.org/en/latest/programs/gdal_raster_info.html)
+(the modern equivalent of `gdalinfo`).
 
-!!! tip
+```bash
+gdal raster info /etc/mapserver/data/raster/54752_dtm_1m.tif
+```
 
-    The `COVERAGEID` will be the `LAYER` `NAME`
+The truncated output is shown below.
 
+```
+Driver: GTiff/GeoTIFF
+Files: /etc/mapserver/data/raster/54752_dtm_1m.tif
+       /etc/mapserver/data/raster/54752_dtm_1m.tif.aux.xml
+Size is 5000, 5000
+Coordinate System is:
+PROJCRS["Estonian Coordinate System of 1997",
+...
+    ID["EPSG",3301]]
+Origin = (655000.000000000000000,6475000.000000000000000)
+Pixel Size = (1.000000000000000,-1.000000000000000)
+...
+Image Structure Metadata:
+  LAYOUT=COG
+...
+Center      (  657500.000, 6472500.000) ( 26d41'30.35"E, 58d21'52.44"N)
+Band 1 Block=512x512 Type=Float32, ColorInterp=Gray
+  Min=30.680 Max=83.205
+  Minimum=30.680, Maximum=83.205, Mean=60.461, StdDev=9.347
+  NoData Value=-9999
+  Overviews: 2500x2500, 1250x1250, 625x625, 312x312
+  Metadata:
+    STATISTICS_MAXIMUM=83.205001831055
+...
+```
 
+From the output we can see that the dataset is in the [EPSG:3301](https://spatialreference.org/ref/epsg/3301/) coordinate reference system, with an origin at (655000, 6475000)
+and a pixel size of 1 × 1 (with a negative Y resolution, as is typical for north-up rasters). The `LAYOUT=COG` indicates that the file is structured as a Cloud-Optimized GeoTIFF (COG).
+
+## Configuring a Mapfile for WCS
+
+The Mapfile for the WCS service is similar to a WMS Mapfile, but with some differences. The `LAYER` type is set to `RASTER`, and the `METADATA` section contains
+keywords prefixed with `wcs_` to specify the WCS parameters.
+
+```scala
+WEB
+    METADATA
+        "wcs_enable_request" "*"
+        "wcs_srs" "EPSG:4326 EPSG:3857"
+        "wcs_title" "Example WCS Mapfile"
+        "wcs_description" "Test description"
+        "wcs_onlineresource" "http://localhost:7000/"
+    END
+END
 ```
-msWCSGetCoverage20(): WCS server error. Raster size out of range, width and height of resulting coverage must be no more than MAXSIZE=4096.
+
+If the Mapfile is used for multiple services such as WMS and WCS, a single metadata item can be specified using the `ows_` prefix, for example `ows_title`.
+
+[LAYER METADATA](https://mapserver.org/ogc/wcs_server.html#layer-object-metadata) can also be used to specify additional information about the coverage, 
+but is not required for this tutorial.
+
+The [OUTPUTFORMAT](https://mapserver.org/mapfile/outputformat.html) defines the properties of the output format.
+In this case we are defining a custom output format for GeoTIFFs with a `FLOAT32` data type to match the source raster
+and ensure that the full precision of the source raster is preserved in WCS responses.
+
+```scala
+OUTPUTFORMAT
+    NAME "GEOTIFF"
+    DRIVER "GDAL/GTiff"
+    MIMETYPE "image/tiff"
+    IMAGEMODE FLOAT32
+    EXTENSION "tif"
+END
 ```
 
-Set the [MAXSIZE](https://mapserver.org/mapfile/map.html#mapfile-map-maxsize) directive on the `MAP` to a larger value. By default this is set to 4096.
+We can use the full power of GDAL to define custom output formats. For example, we could define a COG output format
+by switching to the [COG Driver](https://gdal.org/en/latest/drivers/raster/cog.htm), and add statistics to the output file by adding the `STATISTICS=YES` format option:
+
+```scala
+OUTPUTFORMAT
+    NAME "GEOTIFF_COG"
+    DRIVER "GDAL/COG"
+    MIMETYPE "image/tiff"
+    IMAGEMODE FLOAT32
+    EXTENSION "tif"
+    FORMATOPTION "STATISTICS=YES"
+END
+```
+
+## Requesting a WCS in OpenLayers
 
-[WCS and NULL Values](https://github.com/geographika/wcs-test)
+Typically WCS requests are made from client applications such as QGIS, ArcGIS Pro, or custom JS code in web applications to download the raw raster data,
+rather than to display it as a map image. However, for the purposes of this tutorial we will be using OpenLayers to make requests to the WCS and display the results.
+WCS is not natively supported in OpenLayers, but we can use the [ImageWMS](https://openlayers.org/en/latest/apidoc/module-ol_source_ImageWMS.html) source as a workaround
+by overriding the request parameters to call WCS, and display the results as an image layer on the map.
+
+!!! tip
 
+    The `COVERAGEID` corresponds to the MapServer `LAYER` `NAME`
 
+## Code
+
+??? JavaScript "wcs.js"
+
+    ```js
+    --8<-- "wcs.js"
+    ```
+
+??? Mapfile "wcs.map"
+
+    ``` scala title="wcs.map"
+    --8<-- "wcs.map"
+    ```
+
+## Exercises
+
+1. From the command line, test the WCS 2.0.1 protocol by making a `GetCoverage` request and saving the output as a GeoTIFF using the configured `OUTPUTFORMAT` (MapServer format name, not a MIME type).
+   Then use `gdal raster info` to check the output file.
+
+```
+mapserv -nh "QUERY_STRING=map=/etc/mapserver/wcs.map&SERVICE=WCS&VERSION=2.0.1&REQUEST=GetCoverage&COVERAGEID=dtm&FORMAT=GEOTIFF&SUBSETTINGCRS=http://www.opengis.net/def/crs/EPSG/0/4326&SUBSET=x(26.6507,26.7362)&SUBSET=y(58.3414,58.3879)&SCALESIZE=x(400),y(400)" \
+> output.tif
+gdal raster info output.tif
 ```
-mapserv -nh "QUERY_STRING=map=test.map&SERVICE=WCS&VERSION=2.0.0&REQUEST=GetCoverage&CoverageId=test&FORMAT=GEOTIFF_INT16&BBOX=-69.955,3.420,-69.701,3.5896&CRS=EPSG:4326&WIDTH=500&HEIGHT=500" > output2.tif
-gdalinfo output.tif
+
+2. Add the COG output format to the Mapfile and make a `GetCoverage` request to download a COG-formatted output. Check the output file with `gdal raster info` to see the difference in metadata.
+
+3. Update the JavaScript code to test the WCS 1.0.0 protocol. This requires different parameters to be passed in the requests, 
+   for example `COVERAGEID` becomes `COVERAGE`, and the CRS parameters are different. You can also remove the entire `imageLoadFunction` as WCS 1.0.0
+   more closely matches the WMS protocol, using `BBOX`,`WIDTH`, and `HEIGHT` parameters to specify the area and size of the output image.
+
+```js
+params: {
+    SERVICE: 'WCS',
+    VERSION: '1.0.0',
+    REQUEST: 'GetCoverage',
+    FORMAT: 'image/png',
+    COVERAGE: 'dtm',
+    CRS: 'EPSG:3857',
+    RESPONSE_CRS: 'EPSG:3857',
+},
+```
+
+## Possible Errors
+
 ```
+msWCSGetCoverage20(): WCS server error. Raster size out of range, width and height of resulting coverage must be no more than MAXSIZE=4096.
+```
+
+Set the [MAXSIZE](https://mapserver.org/mapfile/map.html#mapfile-map-maxsize) directive in the `MAP` to a larger value. By default, this is set to 4096.
+
+## Further Reading
 
-<!--
-https://openlayers.org/workshop/en/cog/
--->
+- [MapServer WCS Use Cases](https://mapserver.org/ogc/wcs_format.html)
+- [WCS and NULL Values](https://github.com/geographika/wcs-test)

workshop/content/mkdocs.yml

@@ -25,6 +25,7 @@ nav:
     - Outputs:
         - WMS: outputs/wms.md
         - WFS: outputs/wfs.md
+        - WCS: outputs/wcs.md
         - Tiles: outputs/tiles.md
         - Vector Tiles: outputs/vector-tiles.md
         - OGC API - Features: outputs/ogcapi-features.md
@@ -35,6 +36,7 @@ nav:
         - Clusters: advanced/clusters.md
         - SLD: advanced/sld.md
         - STAC: advanced/stac.md
+        # - WCS and non-EPSG Projections: advanced/wcs-projections.md
         # - Apache: advanced/apache.md
         # - MapScript: advanced/mapscript.md
     - Summary: summary.md

workshop/exercises/app/index.html

@@ -22,10 +22,11 @@ <h2>Inputs</h2>
         <li><a href="raster.html">Raster</a></li>
         <li><a href="stars.html">Stars (Vector)</a></li>
         <li><a href="postgis.html">PostGIS (Databases)</a></li>
-    </ul>    
+    </ul>
     <h2>Outputs</h2>
     <ul>
         <li><a href="wfs.html">WFS</a></li>
+        <li><a href="wcs.html">WCS</a></li>
         <li><a href="tiles.html">Tiles Mode</a></li>
         <li><a href="vector-tiles.html">Vector Tiles</a></li>
         <li><a href="ogcapi-features.html">OGC API - Features</a></li>
@@ -36,10 +37,12 @@ <h2>Advanced</h2>
         <li><a href="gdalg.html">GDAL Pipelines</a></li>
         <li><a href="railways.html">Vector Symbols (Railways)</a></li>
         <li><a href="clusters.html">Clusters</a></li>
+        <!-- <li><a href="direction.html">Direction</a></li> -->
         <li><a href="landuse.html">Landuse</a></li>
         <li><a href="contours.html">Contours</a></li>
         <li><a href="sld.html">SLD - Styled Layer Descriptor</a></li>
         <li><a href="stac.html">STAC</a></li>
+        <!-- <li><a href="wcs-projections.html">STAC</a></li> -->
     </ul>
     <h2>Miscellaneous</h2>
     <ul>

workshop/exercises/app/js/wcs.js

@@ -0,0 +1,66 @@
+import '../css/style.css';
+import Map from 'ol/Map';
+import View from 'ol/View';
+import ImageLayer from 'ol/layer/Image';
+import ImageWMS from 'ol/source/ImageWMS';
+
+const mapserverUrl = import.meta.env.VITE_MAPSERVER_BASE_URL;
+const mapfilesPath = import.meta.env.VITE_MAPFILES_PATH;
+const url = mapserverUrl + mapfilesPath + 'wcs.map';
+
+const wcsSource = new ImageWMS({
+    url,
+    params: {
+        SERVICE: 'WCS',
+        VERSION: '2.0.1',
+        REQUEST: 'GetCoverage',
+        FORMAT: 'image/png',
+        COVERAGEID: 'dtm',
+        SUBSETTINGCRS: 'http://www.opengis.net/def/crs/EPSG/0/3857',
+        OUTPUTCRS: 'http://www.opengis.net/def/crs/EPSG/0/3857',
+    },
+    projection: 'EPSG:3857',
+    imageLoadFunction: (image, src) => {
+        const srcUrl = new URL(src);
+        const params = srcUrl.searchParams;
+
+        // Get the ImageWMS params
+        const bbox = params.get('BBOX').split(',');
+        const width = params.get('WIDTH');
+        const height = params.get('HEIGHT');
+
+        // Replace with WCS 2.0.1 equivalents
+        params.append('SUBSET', `x(${bbox[0]},${bbox[2]})`);
+        params.append('SUBSET', `y(${bbox[1]},${bbox[3]})`);
+        params.set('SCALESIZE', `x(${width}),y(${height})`);
+
+        // Remove the WMS params
+        params.delete('BBOX');
+        params.delete('WIDTH');
+        params.delete('HEIGHT');
+        params.delete('CRS');
+
+        image.getImage().src = srcUrl.toString();
+    },
+    ratio: 1,
+});
+
+const map = new Map({
+    target: 'map',
+    layers: [
+        new ImageLayer({
+            source: wcsSource
+        }),
+    ],
+    view: new View({
+        projection: 'EPSG:3857',
+        center: [2975862, 8046369],
+        zoom: 14,
+    }),
+});
+
+// apply CSS filters to enhance the contrast in the DTM
+map.once('rendercomplete', () => {
+    document.querySelector('.ol-layer canvas').style.filter =
+        'brightness(2.2) contrast(2) sepia(1) hue-rotate(90deg) saturate(3)';
+});

workshop/exercises/app/wcs.html

@@ -0,0 +1,14 @@
+<!DOCTYPE html>
+<html lang="en">
+<head>
+    <meta charset="UTF-8" />
+    <link rel="icon" type="image/x-icon" href="https://openlayers.org/favicon.ico" />
+    <link rel="stylesheet" href="node_modules/ol/ol.css">
+    <meta name="viewport" content="width=device-width, initial-scale=1.0" />
+    <title>WCS</title>
+</head>
+<body>
+    <div id="map"></div>
+    <script type="module" src="./js/wcs.js"></script>
+</body>
+</html>