2b:null 2c:["$","title",null,{"children":"Prototyping GeoParquet + GEOS in WebAssembly | Kyle Barron | Observable"}] 31:Td66,# Prototyping GeoParquet + GEOS in WebAssembly ## *A potential alternative to GeoRust + GeoArrow?* Christoph Pahmeyer
August 15, 2023 Last month Kyle Barron showed how efficient, vectorized geospatial operations can be performend in the browser in his notebook ["Prototyping GeoRust + GeoArrow in WebAssembly"](https://observablehq.com/@kylebarron/prototyping-georust-geoarrow-in-webassembly). Just like Kyle, I believe that the browser is the ideal universal compute platform, and his prototype highlights just how much potential there is! Especially the performance aspects of the GeoArrow + GeoRust pipeline are promising. As Kyle notes: > *Each of these steps [[outlined in the notebook]](https://observablehq.com/@kylebarron/prototyping-georust-geoarrow-in-webassembly) [are] likely more efficient than ever before in the browser for this number of geometries.* However, GeoRust is not a complete geospatial library yet. While its development is active and seemingly fast paced, it's still missing some of the key functionalities for geospatial analysis such as [buffering geometries](https://en.wikipedia.org/wiki/Buffer_analysis). In Rust, GeoRust can use functions from the computational geometry library [GEOS](https://libgeos.org/) as a fallback, but unfortunately this doesn't work in JS + WebAssembly [yet](https://github.com/kylebarron/geoarrow-rs/pull/108)! In the past weeks, fellow Observable user [Ko Nagase](https://observablehq.com/@sanak) and I spent some time porting GEOS to [WebAssembly](https://webassembly.org/) (Wasm) using [Emscripten](https://emscripten.org/). GEOS is a C/C++ library, and is probably best know for being a core dependency of PostGIS, QGIS, GDAL, and Shapely. You can find our library ([GEOS-WASM](https://chrispahm.github.io/geos-wasm/)) on [Github](https://github.com/chrispahm/geos-wasm) 😊 GEOS is very robust but, as Kyle has pointed out in [his blog post](https://kylebarron.dev/blog/geos-wasm#performance), getting it to perform well in a JS + Wasm context would be a quite difficult undertaking. Hence his enthusiasm for GeoArrow + GeoRust! Knowing the possible downsides, we agreed to only create a very thin wrapper around the low-level [C-API](https://libgeos.org/doxygen/geos__c_8h.html) in order to keep the work minimal. As a result, our library isn't really easy to use, and also doesn't have any vectorized operations (yet), but still I think it has its use-cases! In Kyles notebook, I read that GeoParquet was actually encoding geometries in [WKB](https://en.wikipedia.org/wiki/Well-known_text_representation_of_geometry) by default, which is a fairly good fit for GEOS-WASM since WKB can be transferred to WebAssembly efficiently without serialization. This made me curious: if there already is a format that stores geometries in WKB, at least one part of the serialization issue would vanish. And maybe that would even mean that GEOS-WASM wouldn't be *that* much slower than GeoRust eventually? In order to find out, I rebuild Kyles stunning demo using GeoParquet + GEOS-WASM in this notebook. The demo is visualizing 1 million building footprints with 7 million total coordinates for the U.S. state of Utah, and calculating their area on the fly. I also made a version in turf.js (+ GeoJSON), mostly because I love turf.js for its developer experience, still use it on a regular basis, and wanted to see how it would perform against the others.32:T855,table_perf_comparison = htl.html`

	GeoArrow + GeoRust	GeoParquet + GEOS-WASM	GeoJSON + turf.js
Downloading	4246 ms (46.2 MB)	4671 ms (50.7 MB)	4047 ms (40.5 MB, zipped)
Parsing	2856 ms	1599 ms	3413 ms
Reprojecting	3167 ms	2480 ms	-
Allocating Memory	3 ms	13 ms	-
Geometry creation (GEOS only)	-	1846 ms	-
Area calculation	309 ms	138 ms	289 ms**
Freeing Memory	0 ms	239 ms	-
Color scale calculation	138 ms	139 ms	303 ms
Converting to deck.gl binary format	9 ms	102 ms	3817 ms
Total	10728 ms (+0 ms)	11227 ms (-499 ms)	11869 ms (-1141 ms)

Performance comparison for calculating the area of 1 mio. buildings in a web browser. Table shows times to complete each task in miliseconds (lower is better). Tests were performed on a M1 Macbook, repeated 5x and averaged.
** Turf.js only calculates the approximate area of a polygon on a sphere, which eliminates the need for reprojection. The algorithm is based on the works of Chamberlain and Duquette (2007).

`33:Ta46,As shown in the table, all implementations take around 10-12 seconds to complete the task on my machine. On average, I find GeoArrow + GeoRust to perform the whole task **~499 ms** faster than GEOS-WASM (WKB), and **~1141 ms** faster than turf.js (GeoJSON). It's also important to note that turf.js only returns approximate area results, which drastically simplifies the calculations. The GeoArrow + GeoRust toolchain is still fairly new, so I guess further performance increases are likely to come in the future. I'm also guessing that the geometry reprojection and maybe also the data parsing have some room for improvement. With regards to the data parsing, I haven't figured out why parsing the GeoArrow encoded GeoParquet file takes almost twice as long when compared to WKB encoded GeoParquet. Also, I suspect that reprojecting geometries will eventually be faster in Rust (WASM) than in JS. If parsing and reprojection in GeoArrow + GeoRust was as fast as using WKB + GEOS-WASM (this notebook), using GeoArrow + GeoRust would beat GEOS-WASM by 2.5 seconds! Despite this, it's important to note that performance isn't the only advantage of using GeoArrow + GeoRust. As Kyle notes in [his blog post](https://kylebarron.dev/blog/geos-wasm#ease-of-binding-1), the benefit of not having to convert between formats, the overall ease of development (Rust vs C), and more permissive licensing (MIT vs LGPL 2.1) all speak for GeoRust + GeoArrow. However, until GeoArrow + GeoRust matures, I still see a few reasons for why one might choose GEOS-WASM over GeoRust: As stated above, GEOS works pretty well with WKB, and the GeoParquet [format specification](https://geoparquet.org/releases/v1.0.0-rc.1/) (v1.0.0-rc.1) notes that geometries must be encoded in WKB. While it's possible to encode geometries using GeoArrow as well, GeoParquet files encountered on the web are probably more likely to follow the standard WKB encoding. For these files, using GEOS-WASM is a good fit as outlined in this notebook. Also GEOS is a well tested library with lots of geospatial functions that can be used in the browser today using GEOS-WASM. Especially if multiple operations need to be done on the geometries, I'd say that the fix costs of having to convert WKB to GEOS geometries are tolerable. And for certain operations such as buffering geometries, there aren't really that many alternatives yet. Last but not least, this notebook shows that also performance wise GEOS-WASM doesn't have to hide, so I'd argue it's a viable alternative to GeoArrow + GeoRust if you need performant geospatial operations in the browser today! 34:T858,areas = { const start = performance.now(); // Create a proj4 object for the source and target projections const sourceProj = new proj4.Proj("WGS84"); const targetProj = new proj4.Proj( "+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs +type=crs" ); const reader = geos.GEOSWKBReader_create(); const areas = new Float32Array(geometryColumn.length); // we'll re-use the area pointer const areaPtr = geos.Module._malloc(8); // reproject the entire WKB feature array ("values") console.time("reprojecting"); const values = geometryColumn.data[0].values; const valueOffsets = geometryColumn.data[0].valueOffsets; const reprojectedValues = reprojectedWKBArray( values, valueOffsets, sourceProj, targetProj ); console.timeEnd("reprojecting"); console.time("allocating memory"); const geomPtrs = new Int32Array(geometryColumn.length); const values_ptr = geos.Module._malloc( reprojectedValues.length * reprojectedValues.BYTES_PER_ELEMENT ); geos.Module.HEAPU8.set(reprojectedValues, values_ptr); console.timeEnd("allocating memory"); console.time("creating geos geoms"); for (let i = 0; i < valueOffsets.length - 1; i++) { // Get the pointer and size of the current feature const feature_ptr = values_ptr + valueOffsets[i]; const feature_size = valueOffsets[i + 1] - valueOffsets[i]; // Read the feature from the buffer geomPtrs[i] = geos.GEOSWKBReader_read(reader, feature_ptr, feature_size); } console.timeEnd("creating geos geoms"); console.time("calculating area"); for (let i = 0; i < valueOffsets.length - 1; i++) { geos.GEOSArea(geomPtrs[i], areaPtr); areas[i] = geos.Module.getValue(areaPtr, "double"); } console.timeEnd("calculating area"); console.time("freeing memory"); for (let i = 0; i < geometryColumn.length; ++i) { geos.GEOSGeom_destroy(geomPtrs[i]); } geos.GEOSFree(values_ptr); geos.GEOSFree(areaPtr); geos.GEOSWKBReader_destroy(reader); console.timeEnd("freeing memory"); const end = performance.now(); console.log("total proj + area: " + ~~(end - start) + " ms"); return areas; }35:T7b3,resolved = { console.time("converting WKB array to deckgl binary"); var geomOffsets = geometryColumn.data[0].valueOffsets; var numCoordinates = countWKBcoords( geometryColumn.data[0].values, geomOffsets ); var flatCoordinateArray = new Float64Array(numCoordinates); var polygonIndices = new Int32Array(geomOffsets.length); var index = 0; var lastStartingPosIndex = 0; var wkbPos = 0; var offset = 0; var coordCount = 0; var wkbCopy = new Uint8Array(geometryColumn.data[0].values); // create a DataView to access wkbArray var dataView = new DataView(wkbCopy.buffer); // loop through polygon start indeces polygonIndices[i] = 0; for (var i = 0; i <= geomOffsets.length; i++) { // get the start index of the current geometry wkbPos = geomOffsets[i]; // read the byte order (0 for big endian, 1 for little endian) var byteOrder = 0; //wkbArray[wkbPos]; wkbPos += 5; var numRings; try { numRings = dataView.getUint32(wkbPos, byteOrder === 0); // number of rings } catch (error) { break; } wkbPos += 4; for (var j = 0; j < numRings; j++) { var numPoints = dataView.getUint32(wkbPos, byteOrder === 0); // number of points in each ring for (var k = 0; k < numPoints; k++) { // Each point has one coordinate pair starting from wkbPos + 4 + k * 16 // Read the x and y coordinates as double values from the data view flatCoordinateArray[index++] = dataView.getFloat64( wkbPos + 4 + k * 16, byteOrder === 0 ); flatCoordinateArray[index++] = dataView.getFloat64( wkbPos + 12 + k * 16, byteOrder === 0 ); coordCount += 1; } // Update the wkbPos to the next ring wkbPos += 4 + numPoints * 16; } polygonIndices[i + 1] = coordCount; } console.timeEnd("converting WKB array to deckgl binary"); return { flatCoordinateArray, polygonIndices }; }36:T441,colorScale = { // Convert from value to RGB values const lowerBound = 10; const upperBound = 5000; const linearScale = d3 .scaleLinear() .domain([parseFloat(lowerBound), parseFloat(upperBound)]) .range([0, 1]); // This is a fast vectorized approach: we store all colors in a single Float32Array instead of many small JS buffers // Indices maps from geometry to vertex index return (values, indices) => { const vertexArrayLength = indices[indices.length - 1]; const outputArray = new Float32Array(vertexArrayLength * 3); let lastIndex = 0; for (let i = 0; i < values.length; ++i) { let nextIndex = indices[i + 1]; const value = values[i]; const color = d3.color(d3.interpolateViridis(linearScale(value))); const r = color.r / 255; const g = color.g / 255; const b = color.b / 255; for (let j = lastIndex; j < nextIndex; ++j) { outputArray[j * 3] = r; outputArray[j * 3 + 1] = g; outputArray[j * 3 + 2] = b; } lastIndex = nextIndex; } return outputArray; }; }37:T6c7,function reprojectedWKBArray(wkbArray, geomOffsets, sourceProj, targetProj) { var wkbPos = 0; var wkbCopy = new Uint8Array(wkbArray); // create a DataView to access wkbArray var dataView = new DataView(wkbCopy.buffer); // loop through polygon start indeces for (var i = 0; i <= geomOffsets.length; i++) { // get the start index of the current geometry wkbPos = geomOffsets[i]; // read the byte order (0 for big endian, 1 for little endian) var byteOrder = 0; //wkbArray[wkbPos]; wkbPos += 5; var numRings; try { numRings = dataView.getUint32(wkbPos, byteOrder === 0); // number of rings } catch (error) { break; } wkbPos += 4; for (var j = 0; j < numRings; j++) { var numPoints = dataView.getUint32(wkbPos, byteOrder === 0); // number of points in each ring for (var k = 0; k < numPoints; k++) { // Each point has one coordinate pair starting from wkbPos + 4 + k * 16 // Read the x and y coordinates as double values from the data view var x = dataView.getFloat64(wkbPos + 4 + k * 16, byteOrder === 0); var y = dataView.getFloat64(wkbPos + 12 + k * 16, byteOrder === 0); // Reproject the coordinates using proj4 var reprojected = proj4(sourceProj, targetProj, [x, y]); // Write the reprojected coordinates back to the data view dataView.setFloat64( wkbPos + 4 + k * 16, reprojected[0], byteOrder === 0 ); dataView.setFloat64( wkbPos + 12 + k * 16, reprojected[1], byteOrder === 0 ); } // Update the wkbPos to the next ring wkbPos += 4 + numPoints * 16; } } return wkbCopy; }38:T2653,reprojectWKB = { // A function that reprojects a WKB feature using proj4js // Input: a Uint8Array representing a WKB feature, a source projection, and a target projection // Output: a Uint8Array representing the same feature in the target projection function reprojectWKB(wkb, sourceProj, targetProj) { // Copy the WKB without serialization var wkbCopy = new Uint8Array(wkb); // Create a single DataView to access the array buffer var dataView = new DataView(wkbCopy.buffer); // Extract the byte order from the first byte var byteOrder = 0; //wkbCopy[0]; // Extract the geometry type from the next four bytes var geometryType = dataView.getUint32(1, byteOrder === 0); // Define a lookup table to map geometry types to functions var geometryFunctions = { 1: reprojectPoint, 2: reprojectLineString, 3: reprojectPolygon, 4: reprojectMultiPoint, 5: reprojectMultiLineString, 6: reprojectMultiPolygon }; // Check if the geometry type is supported if (geometryFunctions.hasOwnProperty(geometryType)) { // Call the corresponding function to reproject the coordinates geometryFunctions[geometryType]( dataView, byteOrder, sourceProj, targetProj ); } else { // Throw an error and return null throw new Error("Unsupported geometry type: " + geometryType); return null; } // Return the reprojected WKB return wkbCopy; } // A helper function that reprojects a point coordinate pair from a DataView // Input: a DataView, a byte order (0 for little endian, 1 for big endian), // a source projection object, and a target projection object // Output: none (the data view is modified in place) function reprojectPoint(dataView, byteOrder, sourceProj, targetProj) { // A point has one coordinate pair starting from index 5 // Read the x and y coordinates as double values from the data view var x = dataView.getFloat64(5, byteOrder === 0); var y = dataView.getFloat64(13, byteOrder === 0); // Reproject the coordinates using proj4 var reprojected = proj4(sourceProj, targetProj, [x, y]); // Write the reprojected coordinates back to the data view dataView.setFloat64(5, reprojected[0], byteOrder === 0); dataView.setFloat64(13, reprojected[1], byteOrder === 0); } // A helper function that reprojects a linestring coordinate pairs from a DataView // Input: a DataView, a byte order (0 for little endian, 1 for big endian), // a source projection object, and a target projection object // Output: none (the data view is modified in place) function reprojectLineString(dataView, byteOrder, sourceProj, targetProj) { // A linestring has a number of points starting from index 9 var numPoints = dataView.getUint32(5, byteOrder === 0); for (var i = 0; i < numPoints; i++) { // Each point has one coordinate pair starting from index 9 + i * 16 // Read the x and y coordinates as double values from the data view var x = dataView.getFloat64(9 + i * 16, byteOrder === 0); var y = dataView.getFloat64(17 + i * 16, byteOrder === 0); // Reproject the coordinates using proj4 var reprojected = proj4(sourceProj, targetProj, [x, y]); // Write the reprojected coordinates back to the data view dataView.setFloat64(9 + i * 16, reprojected[0], byteOrder === 0); dataView.setFloat64(17 + i * 16, reprojected[1], byteOrder === 0); } } // A helper function that reprojects a polygon coordinate pairs from a DataView // Input: a DataView, a byte order (0 for little endian, 1 for big endian), // a source projection object, and a target projection object // Output: none (the data view is modified in place) function reprojectPolygon(dataView, byteOrder, sourceProj, targetProj) { // A polygon has a number of rings starting from index 9 var numRings = dataView.getUint32(5, byteOrder === 0); var offset = 9; for (var i = 0; i < numRings; i++) { // Each ring has a number of points starting from offset + 4 var numPoints = dataView.getUint32(offset, byteOrder === 0); for (var j = 0; j < numPoints; j++) { // Each point has one coordinate pair starting from offset + 4 + j * 16 // Read the x and y coordinates as double values from the data view var x = dataView.getFloat64(offset + 4 + j * 16, byteOrder === 0); var y = dataView.getFloat64(offset + 12 + j * 16, byteOrder === 0); // Reproject the coordinates using proj4 var reprojected = proj4(sourceProj, targetProj, [x, y]); // Write the reprojected coordinates back to the data view dataView.setFloat64( offset + 4 + j * 16, reprojected[0], byteOrder === 0 ); dataView.setFloat64( offset + 12 + j * 16, reprojected[1], byteOrder === 0 ); } // Update the offset to the next ring offset += 4 + numPoints * 16; } } // A helper function that reprojects a multipoint coordinate pairs from a DataView // Input: a DataView, a byte order (0 for little endian, 1 for big endian), // a source projection object, and a target projection object // Output: none (the data view is modified in place) function reprojectMultiPoint(dataView, byteOrder, sourceProj, targetProj) { // A multipoint has a number of points starting from index 9 var numPoints = dataView.getUint32(5, byteOrder === 0); for (var i = 0; i < numPoints; i++) { // Each point has one coordinate pair starting from index 13 + i * 21 // Read the x and y coordinates as double values from the data view var x = dataView.getFloat64(13 + i * 21, byteOrder === 0); var y = dataView.getFloat64(21 + i * 21, byteOrder === 0); // Reproject the coordinates using proj4 var reprojected = proj4(sourceProj, targetProj, [x, y]); // Write the reprojected coordinates back to the data view dataView.setFloat64(13 + i * 21, reprojected[0], byteOrder === 0); dataView.setFloat64(21 + i * 21, reprojected[1], byteOrder === 0); } } // A helper function that reprojects a multilinestring coordinate pairs from a DataView // Input: a DataView, a byte order (0 for little endian, 1 for big endian), // a source projection object, and a target projection object // Output: none (the data view is modified in place) function reprojectMultiLineString( dataView, byteOrder, sourceProj, targetProj ) { // A multilinestring has a number of linestrings starting from index 9 var numLineStrings = dataView.getUint32(5, byteOrder === 0); var offset = 9; for (var i = 0; i < numLineStrings; i++) { // Each linestring has a number of points starting from offset + 9 var numPoints = dataView.getUint32(offset + 5, byteOrder === 0); for (var j = 0; j < numPoints; j++) { // Each point has one coordinate pair starting from offset + 9 + j * 16 // Read the x and y coordinates as double values from the data view var x = dataView.getFloat64(offset + 9 + j * 16, byteOrder === 0); var y = dataView.getFloat64(offset + 17 + j * 16, byteOrder === 0); // Reproject the coordinates using proj4 var reprojected = proj4(sourceProj, targetProj, [x, y]); // Write the reprojected coordinates back to the data view dataView.setFloat64( offset + 9 + j * 16, reprojected[0], byteOrder === 0 ); dataView.setFloat64( offset + 17 + j * 16, reprojected[1], byteOrder === 0 ); } // Update the offset to the next linestring offset += 9 + numPoints * 16; } } // A helper function that reprojects a multipolygon coordinate pairs from a DataView // Input: a DataView, a byte order (0 for little endian, 1 for big endian), // a source projection object, and a target projection object // Output: none (the data view is modified in place) function reprojectMultiPolygon(dataView, byteOrder, sourceProj, targetProj) { // A multipolygon has a number of polygons starting from index 9 var numPolygons = dataView.getUint32(5, byteOrder === 0); var offset = 9; for (var i = 0; i < numPolygons; i++) { // Each polygon has a number of rings starting from offset + 9 var numRings = dataView.getUint32(offset + 5, byteOrder === 0); for (var j = 0; j < numRings; j++) { // Each ring has a number of points starting from offset + 13 var numPoints = dataView.getUint32(offset + 9, byteOrder === 0); for (var k = 0; k < numPoints; k++) { // Each point has one coordinate pair starting from offset + 13 + k * 16 // Read the x and y coordinates as double values from the data view var x = dataView.getFloat64(offset + 13 + k * 16, byteOrder === 0); var y = dataView.getFloat64(offset + 21 + k * 16, byteOrder === 0); // Reproject the coordinates using proj4 var reprojected = proj4(sourceProj, targetProj, [x, y]); // Write the reprojected coordinates back to the data view dataView.setFloat64( offset + 13 + k * 16, reprojected[0], byteOrder === 0 ); dataView.setFloat64( offset + 21 + k * 16, reprojected[1], byteOrder === 0 ); } // Update the offset to the next ring offset += 13 + numPoints * 16; } // Update the offset to the next polygon offset += 9 + numRings * 4; } } return reprojectWKB; }39:T745,function countWKBcoords(wkbArray, startIndexArray) { var wkbCopy = new Uint8Array(wkbArray); // create a DataView to access the wkbArray var dataView = new DataView(wkbCopy.buffer); // initialize current position in wkbArray and total number of coordinates var wkbPos = 0; var numCoords = 0; // loop through startIndexArray for (var i = 0; i < startIndexArray.length; i++) { // get the start index of the current geometry wkbPos = startIndexArray[i]; // read the byte order (0 for big endian, 1 for little endian) var byteOrder = 0; //wkbArray[wkbPos]; wkbPos++; // read the geometry type (see https://en.wikipedia.org/wiki/Well-known_text_representation_of_geometry#Well-known_binary) var geomType; try { geomType = dataView.getUint32(wkbPos, byteOrder === 0); } catch (error) { // the last index will be out of range continue; } wkbPos += 4; // determine the number of coordinates based on the geometry type switch (geomType) { case 1: // point numCoords += 1; break; case 2: // linestring var numPoints = dataView.getUint32(wkbPos, byteOrder === 0); // number of points wkbPos += 4; numCoords += numPoints; break; case 3: // polygon var numRings = dataView.getUint32(wkbPos, byteOrder === 0); // number of rings wkbPos += 4; for (var j = 0; j < numRings; j++) { var numPoints = dataView.getUint32(wkbPos, byteOrder === 0); // number of points in each ring wkbPos += 4; numCoords += numPoints; } break; default: // other geometry types are not supported in this example throw new Error("Unsupported geometry type: " + geomType); } } // return the total number of coordinates return numCoords * 2; }3a:[] 2f:{"id":"286090808f0cf4d3","slug":null,"title":"Prototyping GeoParquet + GEOS in WebAssembly","thumbnail":null,"access_level":"public","create_time":"$D2023-08-16T03:24:03.189Z","update_time":"$D2023-09-06T15:13:44.406Z","view_time":null,"first_public_time":"$D2023-08-16T03:24:03.189Z","has_unpublished_edits":null,"fork":{"id":"9146bff3f9a33b74","slug":"prototyping-geoparquet-geos-in-webassembly","title":"Prototyping GeoParquet + GEOS in WebAssembly","owner":{"id":"d2dc6c8d94d1ea2c","login":"chrispahm","name":"Christoph Pahmeyer"}},"stars":0,"starred":false,"owner":{"type":"team","id":"43d330d6f3e9ea97","login":"kylebarron","name":"Kyle Barron","avatar_url":"https://avatars-next.observableusercontent.com/75d37333b9e7068d0c7abf571919ecfa40a57510d4b2cf55d65c9b3bec1f2a9f","tier":"starter_2024"},"version":1089,"publish_version":null,"role":null,"authored":false,"authors":[{"type":"individual","id":"d003201354321004","login":"kylebarron","name":"Kyle Barron","avatar_url":"https://avatars-next.observableusercontent.com/0096ac4ec8894f01af5a4d303fe3463061511b44dda0e4803aafc9ede3f46507","tier":null}],"members":[{"type":"individual","id":"d003201354321004","login":"kylebarron","name":"Kyle Barron","avatar_url":"https://avatars-next.observableusercontent.com/0096ac4ec8894f01af5a4d303fe3463061511b44dda0e4803aafc9ede3f46507","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":false,"trashed":false,"lapsed":null,"subscribed":null,"forks":0,"collection":null,"type":"notebook","body":{"title":"Prototyping GeoParquet + GEOS in WebAssembly","cells":[{"id":673,"value":"$31","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":814,"value":"### Performance comparison results\n\nThe following table shows a comparison for the times needed to perform the different steps involved to calculate the area of 1 million buildings in Utah using different geospatial tools, and visualizing them using deck.gl. You can follow the links in the table headers to get to the respective notebooks!","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":827,"value":"$32","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":833,"value":"$33","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":249,"value":"mapContainer = html`

`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":938,"value":"## How it works\n\nThe notebook is borrowing most of it's structure from [Kyles original notebook](https://observablehq.com/@kylebarron/prototyping-georust-geoarrow-in-webassembly). \nAt a high level, the steps in this notebook are:\n\n- Download GeoParquet file and decode the array of WKB geometries.\n- Convert array of WKB geometries from the WGS84 geographic coordinate system to the UTM projected coordinate system so that area computations provide accurate results in meters, not degrees.\n- Copy the array of WKB geometries to WebAssembly memory.\n- Create GEOS geometries from the array of WKB geometries in WebAssembly.\n- Compute area from the UTM geometries.\n- Free WebAssembly memory.\n- Convert array of WKB geometries to deck.gl binary format\n- Render the geometries to deck.gl, coloring buildings by their area.\n\nMost cells contain timers, so check your browsers console for the timing on your machine!","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":985,"value":"### Pre-processing\n\nWe obtain the building datasets from Microsoft. Microsoft distributes the data as zipped GeoJSON, so as a pre-processing step we need to convert it to GeoParquet. We can use [`ogr2ogr`](https://gdal.org/programs/ogr2ogr.html#ogr2ogr) (version 3.5 or later) to do this.\n```sh\nwget https://usbuildingdata.blob.core.windows.net/usbuildings-v2/Utah.geojson.zip\ndocker run --rm -it -v $(pwd):/data osgeo/gdal:latest \\\n ogr2ogr \\\n /data/Utah.parquet \\\n /vsizip//data/Utah.geojson.zip \\\n -dialect SQLite \\\n -sql \"SELECT geometry FROM 'Utah.geojson'\" \\\n -lco COMPRESSION=BROTLI \\\n -lco POLYGON_ORIENTATION=COUNTERCLOCKWISE \\\n -lco ROW_GROUP_SIZE=9999999\n```\n\nNote how we're not passing the option `-lco GEOMETRY_ENCODING=GEOARROW`, as we explicitly want the geometries to be WKB encoded instead of GeoArrow.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":988,"value":"### Download GeoParquet file and decode to GeoArrow","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":951,"value":"To load the GeoParquet file into memory, we'll use parquet-wasm, Kyle's library that provides WebAssembly bindings to Rust Parquet libraries.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":799,"value":"buildingsUrl = \"https://humusklimanetz-couch.thuenen.de/datasets/Utah.parquet\"","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":25,"value":"parquetBuffer = {\n console.time(\"download\");\n const ab = await fetch(buildingsUrl).then((res) => res.arrayBuffer());\n console.timeEnd(\"download\");\n return ab;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":22,"value":"arrowTable = {\n console.log(\"started parsing\");\n console.time(\"parsing\");\n const parquetBytes = new Uint8Array(parquetBuffer);\n\n console.time(\"parse to ffi table\");\n const wasmArrowTable = parquet.readParquetFFI(parquetBytes);\n console.timeEnd(\"parse to ffi table\");\n\n const recordBatches = [];\n for (let i = 0; i < wasmArrowTable.numBatches(); i++) {\n // Note: Unless you know what you're doing, setting `true` below is recommended to _copy_\n // table data from WebAssembly into JavaScript memory. This may become the default in the\n // future.\n console.time(\"parse record batch\");\n const recordBatch = arrowJsFFI.parseRecordBatch(\n parquetMemory.buffer,\n wasmArrowTable.arrayAddr(i),\n wasmArrowTable.schemaAddr(),\n false\n );\n console.timeEnd(\"parse record batch\");\n\n recordBatches.push(recordBatch);\n }\n\n const table = new arrow.Table(recordBatches);\n\n console.timeEnd(\"parsing\");\n return table;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":960,"value":"This dataset has only a single column, the geometry column, so we'll select that.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":35,"value":"geometryColumn = arrowTable.getChildAt(0)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":973,"value":"By checking the `length` property of the geometryColumn, we can see that there is a total of ${geometryColumn.length} geometries stored in the GeoParquet file","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":970,"value":"geometryColumn.length","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":979,"value":"We can access single features by using the `get` accessor","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":38,"value":"firstWKB = geometryColumn.get(0)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":965,"value":"And access the entire array of WKB features by selecting the `data[0].values` property","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":632,"value":"geometryColumn.data[0].valueOffsets","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":983,"value":"In order to calculate the areas of the buldings, we first reproject the features from the WGS84 geographic coordinate system to the UTM projected coordinate system so that area computations provide accurate results in meters, not degrees. \n\nWe then pass the reprojected Uint8 WKB array to WebAssembly, and create GEOS geometries for each of the features. \n\nAfterwards, we calculate the area for each feature, and read the values from WebAssembly to a Float32Array in JavaScript.\n\nIn a final step, we free all memory used by GEOS in the WebAssembly context.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":43,"value":"$34","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1003,"value":"Voilà! We have results! Each value in the area array is the area (in square meters) of the geometry at that index of the polygon array.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":1005,"value":"### Render to deck.gl","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":1008,"value":"In order to render as performant as possible, we pass binary data to Deck.gl. Deck.gl requires a flat array of coordinates, as well as an array indicating where each geometry starts and ends (i.e. the delta from one offset to the next will be greater than 1 if it has holes).","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":262,"value":"$35","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1012,"value":"The following cell defines a function that applies the viridis color map on the array values. As a quick hack I normalize area values by taking the square root before applying a linear color scale. Maybe using a log scale would make a better visualization, but this works.\n\nNote that this is more complex than generating a color scale for points, like in the previous notebook. This is because when supplying binary attributes deck.gl needs the length of the attributes array to match the length of the vertex array, not the number of geometries. Since we're visualizing polygons, those are not the same. So this cell also needs as input the offsets mapping from geometry to vertices so that it knows how many vertices exist per geometry.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":339,"value":"$36","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1015,"value":"Then we use the function defined in the previous cell on the area values. ","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":341,"value":"colorAttribute = {\n // On my computer, generating the color scale from the area array takes ~130ms.\n // Check your browser console for your timing.\n console.time(\"colorScale\");\n const result = colorScale(areas, resolved.polygonIndices);\n console.timeEnd(\"colorScale\");\n return result;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1017,"value":"Create the deck.gl layer and set it on the map above:","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":247,"value":"deckglLayer = {\n // Refer to https://deck.gl/docs/api-reference/layers/solid-polygon-layer#use-binary-attributes\n const data = {\n // Number of geometries\n length: geometryColumn.length,\n // Indices into coordinateArray where each polygon starts\n startIndices: resolved.polygonIndices,\n attributes: {\n // Flat coordinates array\n getPolygon: { value: resolved.flatCoordinateArray, size: 2 },\n // Pass in the color values per coordinate vertex\n getFillColor: { value: colorAttribute, size: 3 }\n }\n };\n const layer = new deck.SolidPolygonLayer({\n // This is an Observable hack - changing the id will force the layer to refresh when the cell reevaluates\n id: `layer-${Date.now()}`,\n data,\n // Skip normalization for binary data\n _normalize: false,\n // Counter-clockwise winding order\n _windingOrder: \"CCW\"\n });\n\n deckglMap.setProps({ layers: [layer] });\n\n return layer;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":242,"value":"deckglMap = {\n // This is an Observable hack: clear previously generated content\n mapContainer.innerHTML = \"\";\n\n return new deck.DeckGL({\n // The HTML container to render into\n container: mapContainer,\n map: mapboxgl,\n mapStyle:\n \"https://basemaps.cartocdn.com/gl/positron-nolabels-gl-style/style.json\",\n\n // Viewport settings\n initialViewState: {\n longitude: -111.89479567035976,\n latitude: 40.69352652341092,\n zoom: 12.739031655851761,\n pitch: 0,\n bearing: 0\n },\n controller: true\n });\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":0,"value":"## WKB helper functions","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":581,"value":"$37","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":412,"value":"$38","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":346,"value":"$39","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":806,"value":"## Imports","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":18,"value":"proj4 = require(\"proj4\")","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":20,"value":"// Load the parquet-wasm library\nparquetModule = {\n const parquetModule = await import(\n \"https://unpkg.com/parquet-wasm@0.4.0/esm/arrow2.js\"\n );\n // Need to await the default export first to initialize the WebAssembly code\n const { memory } = await parquetModule.default();\n return [parquetModule, memory];\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1082,"value":"parquet = parquetModule[0]","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1085,"value":"parquetMemory = parquetModule[1]","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1087,"value":"arrowJsFFI = require(\"https://unpkg.com/arrow-js-ffi@0.3.0/dist/arrow-js-ffi.umd.js\")","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":30,"value":"arrow = require(\"apache-arrow\")","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":166,"value":"geos = {\n const initGeosJs = (\n await import(\"https://cdn.skypack.dev/geos-wasm@1.1.5?min\")\n ).default;\n const geos = await initGeosJs();\n return geos;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":244,"value":"deck = require.alias({\n h3: {}\n})(\"deck.gl@8.8.6/dist.min.js\")","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":240,"value":"mapboxgl = require(\"mapbox-gl@1.6.0/dist/mapbox-gl.js\")","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":757,"value":"## References","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":749,"value":"ref1 = md`Robert. G. Chamberlain and William H. Duquette, \"Some Algorithms for Polygons on a Sphere\", JPL Publication 07-03, Jet Propulsion Laboratory, Pasadena, CA, June 2007, URL: https://dataverse.jpl.nasa.gov/file.xhtml?fileId=45024&version=1.1`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W3a","stdlib":"1","titling":"auto"}} 30:"$2f:body"