Skip to content

Commit bef67cc

Browse files
geographikageographika
authored
Merge pull request #8 from geographika/wcs-proj
Add non-EPSG projections tutorial
2 parents 5e9b545 + 155ce21 commit bef67cc

9 files changed

Lines changed: 624 additions & 18 deletions

File tree

‎workshop/content/docs/advanced/other-projections.md‎

Lines changed: 281 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,281 @@
1+
# Working with non-EPSG Coordinate Reference Systems
2+
3+
## Overview
4+
5+
Most geospatial workflows use EPSG-defined coordinate reference systems. However, there are cases where we need to work with custom or non-standard projections that are not included in the EPSG registry.
6+
7+
In this example we will be working with the [ESRI:53009](https://spatialreference.org/ref/esri/53009/) - a spherical [Mollweide projection](https://en.wikipedia.org/wiki/Mollweide_projection),
8+
and using WMS, WFS, and WCS services to access data in this projection.
9+
10+
<div class="map">
11+
<iframe src="https://mapserver.github.io/getting-started-with-mapserver-demo/other-projections.html"></iframe>
12+
</div>
13+
14+
`ESRI:53009` can be identified using different formats, as shown in the table below:
15+
16+
17+
| Format | Example |
18+
|--------|--------|
19+
| Authority code | ESRI:53009 |
20+
| OGC URN | urn:ogc:def:crs:ESRI::53009 |
21+
| OGC HTTP URI | http://www.opengis.net/def/crs/ESRI/0/53009 |
22+
23+
## The Mapfile
24+
25+
Projections are referenced in several places in the Mapfile, listed below.
26+
27+
The top-level projection of the Mapfile is set to `ESRI:53009`:
28+
29+
``` scala
30+
PROJECTION
31+
"ESRI:53009"
32+
END
33+
```
34+
35+
The source data used for the "raster" and "countries" layers are in the `EPSG:4326` projection.
36+
These need to be set explicitly in the Mapfile using the `PROJECTION` object within the `LAYER` definition, as follows:
37+
38+
``` scala
39+
LAYER
40+
NAME "raster"
41+
TYPE RASTER
42+
EXTENT -180 -90 180 90
43+
DATA "data/naturalearth/NE2_50M_SR_SMALL.tif"
44+
PROJECTION
45+
"EPSG:4326"
46+
END
47+
END
48+
```
49+
50+
The "cities" layer uses a GeoJSON file already in the `ESRI:53009` projection. Even though it is not strictly necessary to specify the projection for this layer,
51+
it is good practice to do so to ensure that MapServer can correctly identify the source CRS and reproject the data as needed.
52+
53+
```scala
54+
LAYER
55+
NAME "cities"
56+
TYPE POINT
57+
PROJECTION
58+
"ESRI:53009"
59+
END
60+
CONNECTIONTYPE OGR
61+
CONNECTION "data/naturalearth/worldcity_53009.geojson"
62+
```
63+
64+
!!! note
65+
66+
GeoJSON traditionally assumes EPSG:4326. Using other CRS relies on client/server
67+
agreement and is not part of the official GeoJSON specification.
68+
69+
Finally we want to make the projection available to all the WxS protocols - WMS, WFS, and WCS. This is done using the `ows_srs` metadata item in the `WEB` section of the Mapfile, as follows:
70+
71+
``` scala
72+
WEB
73+
METADATA
74+
"ows_srs" "ESRI:53009 EPSG:4326"
75+
END
76+
END
77+
```
78+
79+
Using the `ows_` prefix makes the projection available to all OWSs (OGC Web Services), and avoids having to have duplicate entries such as `wms_srs`, `wfs_srs`, and `wcs_srs`
80+
if we were to use the service-specific metadata items.
81+
82+
!!! note
83+
84+
The order of the layers in the Mapfile is important when requesting a map directly using the MapServer CGI interface
85+
(for example <http://localhost:7000/?map=/etc/mapserver/other-projections.map&mode=map&layer=countries&layer=cities&layer=raster>).
86+
The raster layer must be first in the Mapfile, otherwise the raster will be drawn on top of the vector layers and obscure them. The order of the layers in the request
87+
itself does not affect the rendering order, only the order in the Mapfile.
88+
89+
## OpenLayers
90+
91+
The OpenLayers client for this tutorial is set up to make requests to the WMS, WFS, and WCS services using the `ESRI:53009` projection.
92+
The client is configured to request data in this projection and display it on the map.
93+
94+
OpenLayers does not include definitions for all projections by default, so we use the [proj4js](https://proj4js.org/) library
95+
to define the `ESRI:53009` projection and register it with OpenLayers.
96+
97+
```js
98+
proj4.defs('ESRI:53009',
99+
'+proj=moll +lon_0=0 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs'
100+
);
101+
register(proj4);
102+
103+
const coverageExtent = [-18019909, -9009954, 18019909, 9009954];
104+
105+
const mollweideProjection = new Projection({
106+
code: 'ESRI:53009',
107+
units: 'm',
108+
extent: coverageExtent,
109+
worldExtent: [-180, -90, 180, 90],
110+
global: true,
111+
});
112+
113+
addProjection(mollweideProjection);
114+
```
115+
116+
We set the projection on the map object to `ESRI:53009` so that all requests sent to MapServer return data in this projection.
117+
118+
```js
119+
const map = new Map({
120+
target: 'map',
121+
layers: [imageLayer, wmsLayer, graticule, wfsLayer],
122+
view: new View({
123+
projection: 'ESRI:53009',
124+
center: [0, 0],
125+
zoom: 1,
126+
}),
127+
});
128+
```
129+
130+
To help visually confirm the reprojection is working correctly we add a [graticule layer](https://openlayers.org/en/latest/apidoc/module-ol_layer_Graticule-Graticule.html) to the map, which shows the lines of latitude and longitude.
131+
The graticule uses the `ESRI:53009` projection of the `map` object, and display lines representing degrees of latitude and longitude. We also enable labels to show the coordinates of the graticule lines.
132+
133+
```js
134+
const graticule = new Graticule({
135+
strokeStyle: new Stroke({
136+
color: 'rgba(50,50,50,0.5)',
137+
width: 1,
138+
}),
139+
showLabels: true,
140+
wrapX: false
141+
});
142+
```
143+
144+
### WFS
145+
146+
The "cities" point layer is displayed as a WFS layer in the OpenLayers client.
147+
As w"e construct the URLs to send to MapServer ourselves we need to set the `SRSNAME` to `ESRI:53009` (MapServer interprets the `BBOX` values in the CRS specified by `SRSNAME`).
148+
149+
We need to set the `dataProjection` to `ESRI:53009` in the `GeoJSON` format options to ensure the features are correctly reprojected and displayed on the map.
150+
151+
```js
152+
const wfsSource = new VectorSource({
153+
format: new GeoJSON({
154+
dataProjection: 'ESRI:53009',
155+
}),
156+
url: function (extent) {
157+
const [minx, miny, maxx, maxy] = extent;
158+
return `${url}&SERVICE=WFS&VERSION=2.0.0&REQUEST=GetFeature` +
159+
`&TYPENAMES=ms:cities` +
160+
`&OUTPUTFORMAT=geojson` +
161+
`&SRSNAME=ESRI:53009` +
162+
`&BBOX=${minx},${miny},${maxx},${maxy}`;
163+
},
164+
strategy: bbox,
165+
});
166+
```
167+
168+
### WMS
169+
170+
The "countries" polygon layer is displayed as a WMS. As the OpenLayers map itself is set to use the `ESRI:53009` projection we don't need to
171+
specify any further parameters here.
172+
173+
```js
174+
const wmsLayer = new ImageLayer({
175+
source: new ImageWMS({
176+
url: url,
177+
params: {
178+
LAYERS: 'countries',
179+
FORMAT: 'image/png',
180+
TRANSPARENT: true,
181+
VERSION: '1.3.0',
182+
},
183+
serverType: 'mapserver',
184+
}),
185+
});
186+
```
187+
188+
## WCS
189+
190+
!!! note
191+
192+
OpenLayers does not natively support WCS. In this example we use the same approach as in the [Web Coverage Services (WCS)](../outputs/wcs.md) tutorial for testing the WCS protocol.
193+
Images are requested as PNGs using the [ImageWMS](https://openlayers.org/en/latest/apidoc/module-ol_source_ImageWMS-ImageWMS.html) class and a custom
194+
[imageLoadFunction](https://openlayers.org/en/latest/apidoc/module-ol_Image.html#~LoadFunction),
195+
as displaying GeoTIFFs directly in OpenLayers is not supported without additional libraries.
196+
197+
We construct the WCS request parameters manually in the `imageLoadFunction` and ensure that the `SUBSETTINGCRS` and `OUTPUTCRS` parameters are set to `http://www.opengis.net/def/crs/ESRI/0/53009`
198+
to ensure the data is returned in the correct projection. In this example we use the OGC HTTP URI form of the CRS identifier for the CRS requests.
199+
200+
201+
```js
202+
const wcsSource = new ImageWMS({
203+
url,
204+
params: {
205+
SERVICE: 'WCS',
206+
VERSION: '2.0.1',
207+
REQUEST: 'GetCoverage',
208+
FORMAT: 'image/png',
209+
COVERAGEID: 'raster',
210+
SUBSETTINGCRS: 'http://www.opengis.net/def/crs/ESRI/0/53009',
211+
OUTPUTCRS: 'http://www.opengis.net/def/crs/ESRI/0/53009',
212+
},
213+
imageLoadFunction: (image, src) => {
214+
const srcUrl = new URL(src);
215+
const params = srcUrl.searchParams;
216+
217+
// Get the ImageWMS params
218+
const bbox = params.get('BBOX').split(',');
219+
const width = params.get('WIDTH');
220+
const height = params.get('HEIGHT');
221+
222+
// Replace with WCS 2.0.1 equivalents
223+
params.append('SUBSET', `x(${bbox[0]},${bbox[2]})`);
224+
params.append('SUBSET', `y(${bbox[1]},${bbox[3]})`);
225+
params.set('SCALESIZE', `x(${width}),y(${height})`);
226+
...
227+
```
228+
229+
!!! note
230+
231+
`SUBSETTINGCRS` is not strictly necessary here as MapServer assumes the `SUBSET` parameters are in the same CRS as the `OUTPUTCRS`,
232+
but it is good practice to include it.
233+
234+
## Code
235+
236+
!!! example
237+
238+
- MapServer request: <http://localhost:7000/?map=/etc/mapserver/other-projections.map&mode=map&layer=countries&layer=cities&layer=raster>
239+
- OpenLayers example: <http://localhost:7001/other-projections.html>
240+
241+
??? JavaScript "other-projections.js"
242+
243+
``` js
244+
--8<-- "other-projections.js"
245+
```
246+
247+
??? Mapfile "other-projections.map"
248+
249+
``` scala
250+
--8<-- "other-projections.map"
251+
```
252+
253+
## Exercises
254+
255+
1. Currently the CRS is specified using the `http://www.opengis.net/def/crs/ESRI/0/53009` format,
256+
(OGC CRS HTTP URIs) but this can also be specified using other CRS formats.
257+
258+
Update the OpenLayers WCS request parameters to use this format, and test that the requests are still working correctly.
259+
260+
```js
261+
// OGC URN CRS identifiers
262+
SUBSETTINGCRS: 'urn:ogc:def:crs:ESRI::53009',
263+
OUTPUTCRS: 'urn:ogc:def:crs:ESRI::53009',
264+
265+
// authority codes
266+
SUBSETTINGCRS: 'ESRI:53009',
267+
OUTPUTCRS: 'ESRI:53009',
268+
```
269+
270+
2. Update the Mapfile to use another non-EPSG projection, for example `ESRI:54030` (the [Robinson projection](https://en.wikipedia.org/wiki/Robinson_projection)).
271+
272+
Test everything is configured correctly by making a direct request to the MapServer CGI interface
273+
using <http://localhost:7000/?map=/etc/mapserver/other-projections.map&mode=map&layer=countries&layer=cities&layer=raster>.
274+
275+
276+
## Further Reading
277+
278+
- [Reproject Data with GDAL Tutorial](https://gdal.org/en/latest/tutorials/reprojecting_data.html)
279+
- [MapServer Projections](https://mapserver.org/mapfile/projection.html)
280+
- [Spatial Reference](https://spatialreference.org/)
281+
- [PROJ](https://proj.org/)

‎workshop/content/docs/outputs/wcs.md‎

Lines changed: 36 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -122,6 +122,26 @@ rather than to display it as a map image. However, for the purposes of this tuto
122122
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
123123
by overriding the request parameters to call WCS, and display the results as an image layer on the map.
124124

125+
Typically WCS is used to return raw raster formats such as GeoTIFFs, but viewing these in a browser requires rendering the raw
126+
pixel values using additional JavaScript libraries such as [geotiff.js](https://geotiffjs.github.io/).
127+
128+
In this tutorial we instead request a rendered image (`FORMAT: 'image/png'`) and use [CSS filters](https://developer.mozilla.org/en-US/docs/Web/CSS/Reference/Properties/filter)
129+
to visually confirm that data is being returned correctly.
130+
131+
```js
132+
const wcsSource = new ImageWMS({
133+
url,
134+
params: {
135+
SERVICE: 'WCS',
136+
VERSION: '2.0.1',
137+
REQUEST: 'GetCoverage',
138+
FORMAT: 'image/png',
139+
COVERAGEID: 'dtm',
140+
SUBSETTINGCRS: 'http://www.opengis.net/def/crs/EPSG/0/3857',
141+
OUTPUTCRS: 'http://www.opengis.net/def/crs/EPSG/0/3857',
142+
}
143+
```
144+
125145
!!! tip
126146
127147
The `COVERAGEID` corresponds to the MapServer `LAYER` `NAME`
@@ -145,29 +165,29 @@ by overriding the request parameters to call WCS, and display the results as an
145165
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).
146166
Then use `gdal raster info` to check the output file.
147167
148-
```
149-
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)" \
150-
> output.tif
151-
gdal raster info output.tif
152-
```
168+
```bash
169+
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)" \
170+
> output.tif
171+
gdal raster info output.tif
172+
```
153173
154174
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.
155175
156176
3. Update the JavaScript code to test the WCS 1.0.0 protocol. This requires different parameters to be passed in the requests,
157177
for example `COVERAGEID` becomes `COVERAGE`, and the CRS parameters are different. You can also remove the entire `imageLoadFunction` as WCS 1.0.0
158178
more closely matches the WMS protocol, using `BBOX`,`WIDTH`, and `HEIGHT` parameters to specify the area and size of the output image.
159179
160-
```js
161-
params: {
162-
SERVICE: 'WCS',
163-
VERSION: '1.0.0',
164-
REQUEST: 'GetCoverage',
165-
FORMAT: 'image/png',
166-
COVERAGE: 'dtm',
167-
CRS: 'EPSG:3857',
168-
RESPONSE_CRS: 'EPSG:3857',
169-
},
170-
```
180+
```js
181+
params: {
182+
SERVICE: 'WCS',
183+
VERSION: '1.0.0',
184+
REQUEST: 'GetCoverage',
185+
FORMAT: 'image/png',
186+
COVERAGE: 'dtm',
187+
CRS: 'EPSG:3857',
188+
RESPONSE_CRS: 'EPSG:3857',
189+
},
190+
```
171191
172192
## Possible Errors
173193

‎workshop/content/mkdocs.yml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -36,7 +36,7 @@ nav:
3636
- Clusters: advanced/clusters.md
3737
- SLD: advanced/sld.md
3838
- STAC: advanced/stac.md
39-
# - WCS and non-EPSG Projections: advanced/wcs-projections.md
39+
- Non-EPSG CRSs: advanced/other-projections.md
4040
# - Apache: advanced/apache.md
4141
# - MapScript: advanced/mapscript.md
4242
- Summary: summary.md

‎workshop/exercises/app/index.html‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -42,7 +42,7 @@ <h2>Advanced</h2>
4242
<li><a href="contours.html">Contours</a></li>
4343
<li><a href="sld.html">SLD - Styled Layer Descriptor</a></li>
4444
<li><a href="stac.html">STAC</a></li>
45-
<!-- <li><a href="wcs-projections.html">STAC</a></li> -->
45+
<li><a href="other-projections.html">Working with non-EPSG Coordinate Reference Systems (CRS)</a></li>
4646
</ul>
4747
<h2>Miscellaneous</h2>
4848
<ul>

0 commit comments

Comments
 (0)