|
| 1 | +# Working with STAC - SpatioTemporal Asset Catalogs |
| 2 | + |
| 3 | +## Overview |
| 4 | + |
| 5 | +MapServer can display data directly from a STAC server using GDAL's [STACIT driver](https://gdal.org/en/latest/drivers/raster/stacit.html). |
| 6 | + |
| 7 | +In this example, we will build on the [Download data from a STAC API using GDAL and the command line](https://stacspec.org/en/tutorials/gdal_cli/) tutorial, |
| 8 | +and demonstrate how to render data from Microsoft's Planetary Computer using a STAC connection. |
| 9 | + |
| 10 | +<div class="map"> |
| 11 | + <iframe src="https://mapserver.github.io/getting-started-with-mapserver-demo/stac.html"></iframe> |
| 12 | +</div> |
| 13 | + |
| 14 | +## Checking the Data using GDAL |
| 15 | + |
| 16 | +First, let's inspect the data using GDAL command-line tools. Since GDAL is already installed in the MapServer Docker container, |
| 17 | +we can connect to the container and run the necessary commands: |
| 18 | + |
| 19 | +```bash |
| 20 | +# open a bash shell |
| 21 | +docker exec -it mapserver /bin/bash |
| 22 | +# check the GDAL version |
| 23 | +gdal --version |
| 24 | +# check we have the STACIT driver |
| 25 | +gdal vector --formats | grep STACIT |
| 26 | +# STACIT -raster- (rovs): Spatio-Temporal Asset Catalog Items |
| 27 | +``` |
| 28 | + |
| 29 | +First, let's retrieve information about the 2021 LCMAP primary land cover classification (`lcpri`) asset using GDAL: |
| 30 | + |
| 31 | +```bash |
| 32 | +gdal raster info "STACIT:\"https://planetarycomputer.microsoft.com/api/stac/v1/search?&collections=usgs-lcmap-conus-v13&datetime=2021-01-01/2021-12-31\":asset=lcpri" --output-format text |
| 33 | +``` |
| 34 | + |
| 35 | +We can also provide a bounding box to inspect the file names and details for a specific area: |
| 36 | + |
| 37 | +```bash |
| 38 | +gdal raster info "STACIT:\"https://planetarycomputer.microsoft.com/api/stac/v1/search?&collections=usgs-lcmap-conus-v13&datetime=2021-01-01/2021-12-31&bbox=-73.94658,41.95589,-73.88281,42.00546\":asset=lcpri" --output-format text |
| 39 | +``` |
| 40 | + |
| 41 | +To get statistics for the raster data in the selected area, use the `--stats` flag: |
| 42 | + |
| 43 | +```bash |
| 44 | +gdal raster info "STACIT:\"https://planetarycomputer.microsoft.com/api/stac/v1/search?&collections=usgs-lcmap-conus-v13&datetime=2021-01-01/2021-12-31&bbox=-73.94658,41.95589,-73.88281,42.00546\":asset=lcpri" --stats |
| 45 | +``` |
| 46 | + |
| 47 | +At the end of the output, GDAL lists the bands and their data types. We will use this information later when writing our Mapfile. |
| 48 | + |
| 49 | +```json |
| 50 | +"raster:bands":[ |
| 51 | + { |
| 52 | + "data_type":"uint8", |
| 53 | + "stats":{ |
| 54 | + "minimum":1.0, |
| 55 | + "maximum":8.0, |
| 56 | + "mean":3.475, |
| 57 | + "stddev":1.443 |
| 58 | + }, |
| 59 | + "nodata":0 |
| 60 | + } |
| 61 | +], |
| 62 | +``` |
| 63 | + |
| 64 | +When inspecting the dataset with GDAL, `null` is returned for the `proj:epsg` value. Looking at the `projjson` property, the conversion method is listed as "Albers Equal Area". |
| 65 | +This means the dataset uses a custom projection without an official EPSG code. While MapServer can handle custom projections, for the purposes of this workshop we will use |
| 66 | +[EPSG:5070](https://spatialreference.org/ref/epsg/5070/) in our Mapfile. EPSG:5070 corresponds to NAD83 / Conus Albers, |
| 67 | +a standard Albers Equal Area implementation for the United States, which works well with this dataset. |
| 68 | + |
| 69 | +```json |
| 70 | +"proj:epsg":null, |
| 71 | +"proj:projjson":{ |
| 72 | + "$schema":"https://proj.org/schemas/v0.7/projjson.schema.json", |
| 73 | + "type":"ProjectedCRS", |
| 74 | + "name":"AEA WGS84", |
| 75 | + "base_crs":{ |
| 76 | + "type":"GeographicCRS", |
| 77 | + "name":"WGS 84", |
| 78 | + ... |
| 79 | + }, |
| 80 | + "id":{ |
| 81 | + "authority":"EPSG", |
| 82 | + "code":4326 |
| 83 | + } |
| 84 | + ... |
| 85 | + "conversion":{ |
| 86 | + "name":"unnamed", |
| 87 | + "method":{ |
| 88 | + "name":"Albers Equal Area", |
| 89 | + "id":{ |
| 90 | + "authority":"EPSG", |
| 91 | + "code":9822 |
| 92 | + } |
| 93 | + }, |
| 94 | +``` |
| 95 | + |
| 96 | + |
| 97 | +## The Mapfile |
| 98 | + |
| 99 | +In this Mapfile, we use the `BBOX` parameter sent by WMS clients (such as OpenLayers) to dynamically update the `DATA` clause. We accomplish this using |
| 100 | +[runtime substitution](https://mapserver.org/cgi/runsub.html) and by defining the `BBOX` parameter in the `WEB` `VALIDATION` block. |
| 101 | +A regular expression is used to ensure that the parameter is in the expected format. |
| 102 | + |
| 103 | +```scala |
| 104 | +WEB |
| 105 | + VALIDATION |
| 106 | + # check the parameter is in the form -74.134,41.871,-73.694,42.089 |
| 107 | + BBOX '^-?\d+(?:\.\d+)?,-?\d+(?:\.\d+)?,-?\d+(?:\.\d+)?,-?\d+(?:\.\d+)?$' |
| 108 | + END |
| 109 | +``` |
| 110 | + |
| 111 | +Whenever a `BBOX` parameter is passed to MapServer in a querystring, it becomes available as the dynamic variable `%BBOX%` in the Mapfile. |
| 112 | + |
| 113 | +In this example, we use `%BBOX% ` to dynamically update the STAC `bbox` parameter, allowing the Mapfile to request only the data for the area currently viewed by the WMS client. |
| 114 | + |
| 115 | +```scala |
| 116 | +LAYER |
| 117 | + NAME "lcpri" |
| 118 | + TYPE RASTER |
| 119 | + DATA 'STACIT:\"https://planetarycomputer.microsoft.com/api/stac/v1/search?&collections=usgs-lcmap-conus-v13&datetime=2021-01-01/2021-12-31&bbox=%BBOX%\":asset=lcpri' |
| 120 | +``` |
| 121 | + |
| 122 | +In the OpenLayers client, we use WMS 1.1.1 because its `BBOX` coordinate order matches the longitude/latitude order used by the STAC `bbox` parameter. |
| 123 | + |
| 124 | +By contrast, WMS 1.3.0 uses a latitude/longitude order for its BBOX when using EPSG:4326, which can cause coordinate mismatches if used directly with STAC queries. |
| 125 | + |
| 126 | +```js |
| 127 | +new ImageLayer({ |
| 128 | + source: new ImageWMS({ |
| 129 | + url: mapserverUrl + mapfilesPath + 'stac.map&', |
| 130 | + params: { 'LAYERS': 'lcpri', 'STYLES': '', VERSION: '1.1.1' } |
| 131 | + }), |
| 132 | +``` |
| 133 | + |
| 134 | +## Code |
| 135 | + |
| 136 | +!!! example |
| 137 | + |
| 138 | + - Local OpenLayers example: <http://localhost:7001/stac.html> |
| 139 | + |
| 140 | +??? JavaScript "stac.js" |
| 141 | + |
| 142 | + ``` js |
| 143 | + --8<-- "stac.js" |
| 144 | + ``` |
| 145 | + |
| 146 | +??? Mapfile "stac.map" |
| 147 | + |
| 148 | + ``` scala |
| 149 | + --8<-- "stac.map" |
| 150 | + ``` |
| 151 | + |
| 152 | +## Exercises |
| 153 | + |
| 154 | +- The [metadata](https://data.usgs.gov/datacatalog/metadata/USGS.6759c005d34edfeb8710a490.xml) provides what pixel values represent in the dataset. |
| 155 | + Add new classes to the Mapfile to represent each of the 8 landcover types. |
| 156 | + |
| 157 | + ```xml |
| 158 | + <edom> |
| 159 | + <edomv>3</edomv> |
| 160 | + <edomvd> |
| 161 | + Grass/Shrub. Land predominantly covered with shrubs and perennial or annual natural and domesticated grasses (e.g., pasture), |
| 162 | + forbs, or other forms of herbaceous vegetation. The grass and shrub cover must comprise at least 10% of the area and tree cover is less than 10% of the area. |
| 163 | + </edomvd> |
| 164 | + <edomvds>LCMAP Legend Land Cover Class Descriptions</edomvds> |
| 165 | + </edom> |
| 166 | + ``` |
| 167 | + |
| 168 | +- Comment-out the classes and uncomment the `PROCESSING` directives to create a black and white output of the land cover dataset. |
| 169 | + |
| 170 | +## Further Links |
| 171 | + |
| 172 | +- https://planetarycomputer.microsoft.com/docs/quickstarts/reading-stac/ |
0 commit comments