contour functionality using different projections

[mattijn]

I tried to use the contour transform for other projections than identity for an array that covers the whole world. My source file is a GeoTiff containing the global annual precipitation (for year 2016 based on CFSv2 data): The geographic coordinate system is WGS:84 (epsg:4326), origin (top-left corner) is (-180.02,87.58), pixel size is (1.0,-1.0) for x and y respectively and the file has a width of 360 pixels and a height of 168 pixels. I managed to visualise this file using the identity projection, but once I try to use the mercator projection, I start to see certain issues: While I see that it applies a certain projection transformation, it gets wrongly projected.

I then did some trial and error with the data transformation and prepared 6 datafiles. 3 files had a longitude from -180 to 180, one normal, one upside down, and one upside down and vertically rolled 84 pixels (array height divide by 2`) and 3 files where longitude goes from 0 to 360.

Eventually found out that the version that goes from 0 to 360 in longitude, which is placed upside down and vertically rolled by half the height of the array is positioned correctly on the mercator projection: But somehow everything below 0 degrees latitude is cut-off and there is a stiching issue at the geographical prime meridian. That are two problems, but quite glad I got this far.

Then I tried to change the rotate parmeters. I see that the geographical lines (equator, international date line, arctic- and antartic circle) are reprojecting nicely, but the contours are giving some strange results: Where mercator still has some stable rotations (rotate0), other projections (eg. robinson) have not any rotation which show the contours somehow right.

The following specification were used (vega-editor link):

Click to expand

```json { "$schema": "https://vega.github.io/schema/vega/v4.json", "width": 500, "height": 500, "autosize": "none", "signals": [ { "name": "contour_data", "value": "-180:180", "bind": {"input": "select", "options": [ "-180:180", "-180:180 flipud", "-180:180 flipud roll vert", "0:360", "0:360 flipud", "0:360 flipud roll vert"]} }, { "name": "contour_proj", "value": "identity", "bind": {"input": "select", "options": ["identity", "basemap"]} }, { "name": "basemap_proj", "value": "mercator", "bind": {"input": "select", "options": ["mercator", "robinson"]} }, { "name": "contour_fill", "value": true, "bind": {"input": "radio", "options": [true, false]} }, { "name": "scale", "value": 78, "bind": {"input": "range", "min": 54, "max": 202, "step": 8} }, { "name": "rotate0", "value": 0, "bind": {"input": "select", "options": [-360, -270,-180, -90, 0, 90, 180, 270, 360]} }, { "name": "rotate1", "value": 0, "bind": {"input": "select", "options": [-360, -270,-180, -90, 0, 90, 180, 270, 360]} }, { "name": "rotate2", "value": 0, "bind": {"input": "select", "options": [-180, -135, -90,-45, 0, 45, 90, 135, 180]} }, { "name": "center0", "value": 0, "bind": {"input": "select", "options": [-180, -135, -90,-45, 0, 45, 90, 135, 180]} }, { "name": "center1", "value": 0, "bind": {"input": "select", "options": [-89,-45, 0, 45, 89]} }, { "name": "precipitation", "update": "if(contour_data === '-180:180', data('p-180:180')[0], if(contour_data === '-180:180 flipud', data('p-180:180-flipud')[0], if(contour_data === '-180:180 flipud roll vert', data('p-180:180-flipud-roll-vert')[0], if(contour_data === '0:360', data('p0:360')[0], if(contour_data === '0:360 flipud', data('p0:360-flipud')[0], data('p0:360-flipud-roll-vert')[0])))))" } ], "data": [ { "name": "sphere", "values": [{"type": "Sphere"} ], "transform": [{"type": "geopath", "projection": "proj_basemap" } ] }, { "name": "graticule", "transform": [{"type": "graticule"}, {"type": "geopath", "projection": "proj_basemap"}] }, { "name": "geographic_lines", "url": "https://raw.githubusercontent.com/martynafford/natural-earth-geojson/master/110m/physical/ne_110m_geographic_lines.json", "format": {"type": "json", "property": "features"}, "transform": [{"type": "geopath", "projection": "proj_basemap"}] }, { "name": "world", "url": "https://raw.githubusercontent.com/vega/datalib/master/test/data/world-110m.json", "format": {"type": "topojson", "feature": "countries"}, "transform": [{"type": "geopath", "projection": "proj_basemap"}] }, { "name": "p-180:180", "url": "https://raw.githubusercontent.com/mattijn/datasets/master/annual-precip_-180-180.json" }, { "name": "p-180:180-flipud", "url": "https://raw.githubusercontent.com/mattijn/datasets/master/annual-precip_-180-180_flip-ud.json" }, { "name": "p-180:180-flipud-roll-vert", "url": "https://raw.githubusercontent.com/mattijn/datasets/master/annual-precip_-180-180_flip-ud_roll-vert.json" }, { "name": "p0:360", "url": "https://raw.githubusercontent.com/mattijn/datasets/master/annual-precip_0-360.json" }, { "name": "p0:360-flipud", "url": "https://raw.githubusercontent.com/mattijn/datasets/master/annual-precip_0-360_flip-ud.json" }, { "name": "p0:360-flipud-roll-vert", "url": "https://raw.githubusercontent.com/mattijn/datasets/master/annual-precip_0-360_flip-ud_roll-vert.json" }, { "name": "contours", "transform": [ { "type": "contour", "values": {"signal": "precipitation.values"}, "size": [{"signal": "precipitation.width"}, {"signal": "precipitation.height"}], "smooth": true, "thresholds": {"signal": "sequence(0, 3000, 500)"} } ] } ], "projections": [ { "name": "proj_basemap", "type": {"signal": "basemap_proj"}, "scale": {"signal": "scale"}, "rotate": [{"signal": "rotate0"}, {"signal": "rotate1"}, {"signal": "rotate2"}], "center": [{"signal": "center0"}, {"signal": "center1"}], "translate": [{"signal": "width/2"}, {"signal": "height/2"}] }, { "name": "proj_identity", "type": "identity", "scale": {"signal": "width / precipitation.width"}, "translate": [0, {"signal": "height/4"}] } ], "scales": [ { "name": "color", "type": "sequential", "domain": [0, 3000], "range": {"scheme": "bluepurple"} } ], "marks": [ { "type": "path", "from": {"data": "sphere"}, "encode": { "update": { "fill": {"value": "#fefef6"}, "stroke": {"value": null}, "path": {"field": {"signal": "if(contour_proj === 'basemap', 'path', 'null')"}} } } }, { "type": "path", "from": {"data": "graticule"}, "clip": {"sphere": "proj_basemap"}, "interactive": false, "encode": { "update": { "strokeWidth": {"value": 0.5}, "stroke": {"value": "#ddd"}, "fill": {"value": null}, "path": {"field": {"signal": "if(contour_proj === 'basemap', 'path', 'null')"}} } } }, { "type": "path", "from": {"data": "world"}, "clip": {"sphere": "proj_basemap"}, "interactive": false, "encode": { "update": { "strokeWidth": {"value": 1}, "stroke": {"value": "#dadada"}, "fill": {"value": "#ccc"}, "zindex": {"value": 0}, "path": {"field": {"signal": "if(contour_proj === 'basemap', 'path', 'null')"}} } } }, { "type": "path", "from": {"data": "contours"}, "clip": {"sphere": "proj_basemap"}, "encode": { "update": { "strokeWidth": {"value": 0.5}, "stroke": {"signal": "if(contour_proj === 'identity', null, 'firebrick')"}, "fill": {"signal": "if(contour_proj === 'identity', null, if(contour_fill === false, null, scale('color', datum.value)))"}, "fillOpacity": {"value": 0.5} } }, "transform": [{"type": "geopath", "field": "datum", "projection": "proj_basemap"}] }, { "type": "path", "from": { "data": "contours"}, "encode": { "update": { "strokeWidth": {"value": 0.5}, "stroke": {"signal": "if(contour_proj === 'identity', 'firebrick', null)"}, "fill": {"signal": "if(contour_proj === 'identity', if(contour_fill === false, null, scale('color', datum.value)), null)"}, "fillOpacity": {"value": 0.5} } }, "transform": [{"type": "geopath", "field": "datum", "projection": "proj_identity"}] }, { "type": "path", "from": {"data": "geographic_lines"}, "clip": {"sphere": "proj_basemap"}, "interactive": false, "encode": { "update": { "strokeWidth": {"value": 2}, "stroke": {"value": "green"}, "fill": {"value": null}, "path": {"field": {"signal": "if(contour_proj === 'basemap', 'path', 'null')"}} } } }, { "type": "path", "from": {"data": "sphere"}, "encode": { "update": { "fill": {"value": null}, "stroke": {"value": "black"}, "strokeWidth": {"value": 0.5}, "strokeDash": {"value": [8, 3]}, "path": {"field": {"signal": "if(contour_proj === 'basemap', 'path', 'null')"}} } } } ], "legends": [ { "title": "Annual Precipitation (mm)", "orient": "top-left", "type": "symbol", "fill": "color", "format": ".0f", "clipHeight": 16, "direction":"horizontal", "fillColor": "#fefef6" } ] } ```

[jheer]

Hmm... it looks like you are computing contours first, then sending the result through a projection? I think you want to do the opposite: first project your points to (x, y) coordinates, next compute the contours on those, then use geoshape or geopath with the identity projection to render the contours.

For an example of this approach, take a look at this test specification in vega-lib: https://github.com/vega/vega-lib/blob/master/test/specs-valid/contour-map.vg.json

The example uses the geopoint transform to project the (lon, lat) data points to (x, y) coordinates, then computes contours based on those projected values. Finally, a geoshape transform with the identity projection turns the GeoJSON-formatted contours into rendered on-screen paths. In your case, I think you can simply place a projection step in the data transform pipeline prior to the contour transform, then use the identity projection for your geopath transform in the contour mark definition.

[mattijn]

While I think that it is the correct approach, it is not possible as the input array has no (lon, lat) data points.

Using information of the origin and pixel size I was able to compute the coordinate vectors using sequence transformation, but things seems to get messy in the geopoint transformation.

Click to see the Vega code block

```json { "$schema": "https://vega.github.io/schema/vega/v3.0.json", "signals": [ { "name": "gtf", "update": "data('geotransform')[0]" } ], "data": [ { "name": "geotransform", "values": { "upper_left_x": -180.02, "x_size": 1, "upper_left_y": 87.58, "y_size": -1, "width": 360, "height": 168 } }, { "name": "lat_coords", "transform": [ { "type": "sequence", "start": {"signal": "gtf.upper_left_y + (gtf.y_size / 2)"}, "stop": {"signal": "(gtf.upper_left_y + (gtf.y_size / 2)) + (gtf.y_size * gtf.height)"}, "step": {"signal": "gtf.y_size"} }, { "type": "project", "fields": ["data"], "as": ["lat_coords"] } ] }, { "name": "lon_coords", "transform": [ { "type": "sequence", "start": {"signal": "gtf.upper_left_x + (gtf.x_size / 2)"}, "stop": {"signal": "(gtf.upper_left_x + (gtf.x_size / 2)) + (gtf.x_size * gtf.width)"}, "step": {"signal": "gtf.x_size"} }, { "type": "project", "fields": ["data"], "as": ["lon_coords"] } ] }, { "name": "lonlat_to_xy", "source": ["lat_coords", "lon_coords"], "transform": [ { "type": "geopoint", "projection": "projection", "fields": ["lon_coords", "lat_coords"] } ] } ], "projections": [ { "name": "projection", "type": "mercator" } ] } ```

.

I was not able to find a way to combine these two sequence arrays into a coordinate matrix (in order to flatten the matrix and stack them together).

In numpy (python) this is something that is called meshgrid.

Click to see how its done in Python

```python import numpy as np import pandas as pd ``` ```python upper_left_x = -180.02 x_size = 1 upper_left_y = 87.58 y_size = -1 width = 360 height = 168 values = np.array([392,392,392,392,393,..,166,163,165,168,169]) ``` ```python # get a sequence (coordinate vector) of longitude and latitude values lon = np.arange(upper_left_x + (x_size / 2), (upper_left_x + (x_size / 2)) + (x_size * width), x_size) lat = np.arange(upper_left_y + (y_size / 2), (upper_left_y + (y_size / 2)) + (y_size * height), y_size) ``` ```python print('lon: {}'.format(lon[0:8]),'\nlat: {}'.format(lat[0:8])) lon: [-179.52 -178.52 -177.52 -176.52 -175.52 -174.52 -173.52 -172.52] lat: [87.08 86.08 85.08 84.08 83.08 82.08 81.08 80.08] ``` ```python # get coordinate matrices from coordinate vectors lonlon, latlat = np.meshgrid(lon, lat) # plot both coordinate matrices f, (ax1, ax2) = plt.subplots(1, 2) ax1.imshow(lonlon) ax1.set_title('lon matrix') ax2.imshow(latlat) ax2.set_title('lat matrix') plt.show() ```  ```python # stack the flattened coordinate matrices lonlat = np.vstack([lonlon.flatten(), latlat.flatten()]).T ``` ```python # convert to pandas dataframe df = pd.DataFrame(lonlat, columns=['lon', 'lat']) ``` ```python # print result with pd.option_context('display.max_rows',10): print(df) lon lat 0 -179.52 87.08 1 -178.52 87.08 2 -177.52 87.08 3 -176.52 87.08 4 -175.52 87.08 ... ... ... 60475 175.48 -79.92 60476 176.48 -79.92 60477 177.48 -79.92 60478 178.48 -79.92 60479 179.48 -79.92 [60480 rows x 2 columns] ```

[jheer]

Can't you use a single sequence transform to enumerate all entries, then use formula transformations to populate values? You could derive (x, y) matrix indices from a single sequence value given the matrix width, then calculate lon/lat from there. My snippet below illustrates this, though may not have precisely correct GeoTIFF conversions.

Click to expand Vega JSON snippet

```json "signals": [ { "name": "precip", "update": "data('precip')[0]" }, { "name": "gtf", "value": { "xStart": -180.02, "xStep": 1, "yStart": 87.58, "yStep": -1 } } ], "data": [ { "name": "precip", "url": "https://raw.githubusercontent.com/mattijn/datasets/master/annual-precip_-180-180.json" }, { "name": "coords", "transform": [ { "type": "sequence", "start": 0, "stop": {"signal": "precip.width * precip.height"}, "step": 1 }, { "type": "formula", "as": "lon", "expr": "gtf.xStart + gtf.xStep * (datum.data % precip.width)" }, { "type": "formula", "as": "lat", "expr": "gtf.yStart + gtf.yStep * floor(datum.data / precip.width)" }, { "type": "formula", "as": "value", "expr": "precip.values[datum.data]" }, { "type": "geopoint", "projection": "proj", "fields": ["lon", "lat"], "as": ["x", "y"] } ] } ], ```

[mattijn]

First of all, your code snippet is really amazing! Great way to get the (x, y) matrix indices. I'd to add a half a pixel to the gtf.xStart and gtf.yStart to get the geopointsin the center of the pixels.

Following your test specification contour-map.vg.json, I managed to get a result like this:

Using:

{
  "name": "contours",
  "source": "coords",
  "transform": [
    {
      "type": "contour",
      "x": "x",
      "y": "y",
      "size": [{"signal": "width"}, {"signal": "height"}],
      "count": 10
    }
  ]
}

And I thought I was happy, but then I realised it has created contours of the projected geopoints, which are equidistant from one another. This is more visible using the mercator projection at world scale:

Furthermore, it seems that the usage of values within the contour transform cannot work together with x and y and so the derived contours will not contain the reprojected values. So it just becomes the values and the width and height of the array:

{
  "name": "contours",
  "source": "coords",
  "transform": [
    {
      "type": "contour",
      "values": {"signal": "precip.values"},
      "size": [{"signal": "precip.width"}, {"signal": "precip.height"}],
      "smooth": true,
      "thresholds": {"signal": "sequence(0, 3000, 500)"}
    }
  ]
}  

I always end up in the need of reprojecting the contours, but that is not working, as discussed above.

Therefore I have one more humble question: is there a certain existing working approach that I have overlooked? Maybe something in the identity projection?

My minimum working example is in this gist: https://gist.github.com/mattijn/f166d090dbe9c8a5e03996217f4bd594

[jheer]

The simplest solution would be if the contour transform could incorporate the values (not just the x, y coordinates) into the density estimation process. I checked and there is an open issue for this on d3-contour: https://github.com/d3/d3-contour/issues/28. I filed a PR (https://github.com/d3/d3-contour/pull/29) to see if we can get this feature added to the underlying d3 library.

[jheer]

A new version of d3-contour with a weight parameter for density estimation has been released, as well as an update version of vega-geo. So this functionality will be included in the Vega release.

[Fil]

It is also possible to compute geoContours, staying in spherical coordinates, like we do on https://beta.observablehq.com/@fil/netcdf I reckon it's not properly packaged yet but I thought it worth mentioning here for reference.

[mattijn]

Using the new vega release (congrats!), I tried the weight functionality within contours transform. The behaviour of the contours seems different to the behaviour of the contours in the notebook of @Fil (pretty awesome notebook btw!)

Using a drag rotate event (specs).

The contours in the North Atlantic (between US and EU) seems to indicate that the weight is increasing while doing a mousemove event. That seems not right. As far as I can see, this doesn't happen in @Fil 's notebook.

I also tried to use the timer eventstream to make the earth auto rotate (specs). This is an animated gif, but it also auto rotates in the editor. But here seems another issue that I maybe understand. Since it calculates the x and y coordinates located on screen, also the coordinates on the dark side of the earth are incorporated in the contours transform as these x and y are next to the visible coordinates. (edit: no, I don't think this happens)

Besides that, the timer event worked out reasonably well, two things:

• I would like to be able to stop the timer
• My browser is eaten alive. Quickly the memory increase and the page eventually crash.

Regarding scale of the contour colors, I now use "domain": [0,16] since I checked the values of the contours using

VEGA_DEBUG.view._runtime.subcontext[0].data.contours.values.value

Giving:

(10) [{…}, {…}, {…}, {…}, {…}, {…}, {…}, {…}, {…}, {…}]
0 : {type: "MultiPolygon", value: 1.6611, coordinates: Array(1), Symbol(vega_id): 2807}
1: {type: "MultiPolygon", value: 3.284, coordinates: Array(1), Symbol(vega_id): 2808}
2: {type: "MultiPolygon", value: 4.907, coordinates: Array(3), Symbol(vega_id): 2809}
3: {type: "MultiPolygon", value: 6.530, coordinates: Array(2), Symbol(vega_id): 2810}
4: {type: "MultiPolygon", value: 8.153, coordinates: Array(3), Symbol(vega_id): 2811}
5: {type: "MultiPolygon", value: 9.777, coordinates: Array(4), Symbol(vega_id): 2812}
6: {type: "MultiPolygon", value: 11.400, coordinates: Array(3), Symbol(vega_id): 2813}
7: {type: "MultiPolygon", value: 13.023, coordinates: Array(2), Symbol(vega_id): 2814}
8: {type: "MultiPolygon", value: 14.6468, coordinates: Array(3), Symbol(vega_id): 2815}
9: {type: "MultiPolygon", value: 16.270, coordinates: Array(1), Symbol(vega_id): 2816}
length: 10

I've now idea where these numbers relates to, as the original values goes from 0 to > 3000