Agricultural economics · Thünen Institute of Farm Economics / University Bonn
parquetBuffer = {
console.time("download");
const ab = await fetch(buildingsUrl).then((res) => res.arrayBuffer());
console.timeEnd("download");
return ab;
}
arrowTable = {
console.time("parsing parquet");
const parquetBytes = new Uint8Array(parquetBuffer);
const decodedBytes = parquet.readParquet(parquetBytes);
const arrowTable = arrow.tableFromIPC(decodedBytes);
console.timeEnd("parsing parquet");
return arrowTable;
}
geometryColumn = arrowTable.getChildAt(0)
geometryColumn.length
firstWKB = geometryColumn.get(0)
geometryColumn.data[0].valueOffsets
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;
}
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 };
}
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;
};
}
colorAttribute = {
// On my computer, generating the color scale from the area array takes ~130ms.
// Check your browser console for your timing.
console.time("colorScale");
const result = colorScale(areas, resolved.polygonIndices);
console.timeEnd("colorScale");
return result;
}
deckglLayer = {
// Refer to https://deck.gl/docs/api-reference/layers/solid-polygon-layer#use-binary-attributes
const data = {
// Number of geometries
length: geometryColumn.length,
// Indices into coordinateArray where each polygon starts
startIndices: resolved.polygonIndices,
attributes: {
// Flat coordinates array
getPolygon: { value: resolved.flatCoordinateArray, size: 2 },
// Pass in the color values per coordinate vertex
getFillColor: { value: colorAttribute, size: 3 }
}
};
const layer = new deck.SolidPolygonLayer({
// This is an Observable hack - changing the id will force the layer to refresh when the cell reevaluates
id: `layer-${Date.now()}`,
data,
// Skip normalization for binary data
_normalize: false,
// Counter-clockwise winding order
_windingOrder: "CCW"
});
deckglMap.setProps({ layers: [layer] });
return layer;
}
deckglMap = {
// This is an Observable hack: clear previously generated content
mapContainer.innerHTML = "";
return new deck.DeckGL({
// The HTML container to render into
container: mapContainer,
map: mapboxgl,
mapStyle:
"https://basemaps.cartocdn.com/gl/positron-nolabels-gl-style/style.json",
// Viewport settings
initialViewState: {
longitude: -111.89479567035976,
latitude: 40.69352652341092,
zoom: 12.739031655851761,
pitch: 0,
bearing: 0
},
controller: true
});
}
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;
}
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;
}
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;
}
proj4 = require("proj4")
// Load the parquet-wasm library
parquet = {
const parquetModule = await import(
"https://unpkg.com/parquet-wasm@0.4.0-beta.3/esm/arrow2.js"
);
// Need to await the default export first to initialize the WebAssembly code
await parquetModule.default();
return parquetModule;
}
arrow = require("apache-arrow")
geos = {
const initGeosJs = (
await import("https://cdn.skypack.dev/geos-wasm@1.1.5?min")
).default;
const geos = await initGeosJs();
return geos;
}
deck = require.alias({
h3: {}
})("deck.gl@8.8.6/dist.min.js")
mapboxgl = require("mapbox-gl@1.6.0/dist/mapbox-gl.js")
Purpose-built for displays of data
Observable is your go-to platform for exploring data and creating expressive data visualizations. Use reactive JavaScript notebooks for prototyping and a collaborative canvas for visual data exploration and dashboard creation.
