A low-cost custom wind particle motion layer in mapbox-gl-js

[embed]

There are existing products of web-based wind visualization, notably earth.nullschool.net and Windy. For developers who work on making web maps, Mapbox offers raster-array source and raster-particle style layer that jointly can produce reliable and high-quality particle animation for wind. Additionally, Xweather MapsGL offers the option of adding style-customizable wind particles layer (and many other weather layers) to a mapbox map, and developers do not need to worry about data source issue. Xweather is subscription-based and also bills on extra usage.

For my own work though, I need to be able to accommodate custom wind data sources and overlay the wind layer with Mapbox-hosted vector tilesets, so Mapbox’s raster-particle layer is my first choice. Most of my map’s layers have an extent of a US county and need to zoom in to block level to visualize building and road (zoom level at least 12 and preferably 16). However, raster-array source generates tiles for each zoom level within the user defined range, and a zoom level above 10 (even after clippping of dataset extent) can take extremely long time to process (my process job even crashes after a few hours without finishing), consume a lot of compute units, and entail a significant cost. If I only generate tileset up to zoom level 10, the particle looks really coarse (like a blurry box with tails that slowly moves over my other vector layers) when zoom in. Therefore, I cannot afford Mapbox’s built-in raster-particle layer, given my limited budget and data characteristics (display particles at high zoom level, and input raster cell size preferably is small like 5-km instead of 1 degree).

I turn my attention to open-source solutions such as leaflet-velocity and mapbox-wind, and eventually I decide to implement my own wind particle motion layer based on some key ideas that mention in those works.

My implementation requires storing wind data in JPEG format somewhere that can be downloaded from Internet (e.g., public AWS S3 bucket with proper CORS policy, or as public assets of a static website). If the data needs to be updated daily, the data pipeline can be set up using cloud cron job service or on a cloud virtual machine.

I publish my work as a npm package, and also share my source code, which has a pipeline folder with necessary scripts to prepare data, as well as a react-demo folder that you can run on your machine.

In this article, I will give an overview of key ideas of my implementation.

1. Displacement in layer extent

Particle motion layer can be understood as a layer with points whose positions on the map change over time, instead of stationary like traditional vector layer. A change in position means displacement that has a horizontal component (u-component, or along x-axis) and a vertical component (v-component, or along y-axis).

The typical wind datasets such as GFS have bands with u-component velocity in m/s and v-component velocity in m/s. Vladimir Agafonkin proposes to use JPEG image’s R and G bands to store the normalized u- and v-component velocities (normalize means to scale velocities values to range between 0 and 255). Normalization requires minimum and maximum cell values for u- and v-components, mapbox-wind proposes to store min/max values as an extra attribute of JPEG image — EXIF, so that all we need is just a single image file.

When the JPEG image is passed to the layer render program, it is stored as a WebGL texture, and now the r and g bands become in 0–1 range when retrieving value (a “built-in” normalization from 0–255 to 0–1). The shader program can access this velocity texture via u_velocity_texture.

The minimum and maximum cell values of u-component velocity is stored as a 2-element vector of signed floating point numbers in WebGL and denoted as u_value_range_u. Similar, the min and max of v-component is denoted as u_value_range_v.

In geospatial world, a dataset always has an extent ([minX, maxY, maxX, minY] assumed to be in degrees), and a JPEG file does not store such information. This information needs to be passed as a bounds parameter when adding the layer. The 4 values are stored as a 4-element vector of signed floating point numbers in WebGL, and it is denoted as u_bounds.

I use bounds to normalize point positions so their x- and y-coordinates have a range of [0,1], and I store the normalized point positions as WebGL buffer.

With the above mentioned information, we can compute displacement as well as the new point (i.e., particle) position.

(1) For a point’s current position (denoted as a_position), we need to look up the cell/pixel value that it is located on by sampling the velocity texture.

vec4 velocity = texture(u_velocity_texture, a_position);

(2) We need to scale the velocities with normalized range 0–1 back to their original range.

float u = mix(u_value_range_u[0], u_value_range_u[1], velocity.r);
float v = mix(u_value_range_v[0], u_value_range_v[1], velocity.g);

(3) The next step is to convert the denormalized velocity unit (assumed to be mph in shader program) to degrees per hour (for consistence with layer extent unit). Note that texture stores values from top to bottom, while point position assumes y-axis increasing from bottom to top; this is a critical detail when converting point position to actual latitude (latitude is need to convert to degrees because the length of a degree of longitude is a function of latitude).

float lat = mix(u_bounds[3], u_bounds[1], 1.0 - a_position.y);
vec2 degrees = mphToDegreesPerFrame(vec2(u, v), lat);

(4) New point positions (in the next hour) is computed by current position plus displacement amount. The point positions are in [0, 1] range, so normalization of degrees to [0,1] range is required to reflect the extent of the layer and used to determine the new position. Animation effect is achieved by mapping point positions from one frame to next frame, typically there are around 60 frames per second. A common practice is to scale the velocity vector with a small fraction number (per hour now becomes per frame) before adding it to the current position; this way will take frame rate into account.

vec2 normalizedVelocity = vec2(
 degrees.x / (u_bounds[2] - u_bounds[0]),
 degrees.y / (u_bounds[1] - u_bounds[3])
);
vec2 newPos = a_position + normalizedVelocity * u_speed_factor;

2. Position reset mechanism

After initial new position of a point is computed, we may need to reset this initial position. The reasons for resetting positions can be exceeding layer boundary or preventing degeneration.

Exceeding layer boundary refers to the point’s new position exceeds [0,1] range that defines the bounds.

The concept of degeneration is mentioned in some other works on mapping wind particles. Ideally, particles should be distributed across the entire layer extent; in other words, in a sub-area, you should always see some particles present there (even if the particle count is small compared with some other areas). Over time, particles might be all trapped in specific areas, and most other parts of the extent have little to no particles at all; this phenomenon is degeneration.

There are two main causes of degeneration: circular air flow and extreme low wind speed. Circular air flow mainly is caused by Coriolis effect (Earth’s rotation); visually, particles form a circular pattern, no particle can move out from it while only new particles can move in to it. Examples of circular air flow include hurricanes, typhoons, or low-pressure systems (where you might expect precipitation) in general. Low wind speed means that particles will stuck at a region almost like stationary; examples may include the center of high pressure systems or possibly center of some low pressure systems as well.

Exceeding layer boundary and low wind speed (my shader program assumes below 1.5 mph as low) are easy to handle:

float windSpeed = length(vec2(u, v));
if (newPos.x < 0.0 || newPos.x > 1.0 || newPos.y < 0.0 || newPos.y > 1.0 || windSpeed < 1.5) {
    shouldReset = true;
}

To handle degeneration due to circular air flow, I use WebGL buffer to store particle age. Each particle’s age is initialized with a random number between 0 and 100 so that a large number of particles would not be updated at the same time.

There are two parameters related to ageThreshold (denoted as u_age_threshold in shader program) and maxAge (denoted as u_max_age in shader program). The more particle age exceeds ageThreshold, the more likely for it to be reset. Eventually if particle age exceeds maxAge, it must be reset.

float age = a_age + 1.0;

if (age > u_age_threshold) {
    float resetProbability = (age - u_age_threshold) / (u_max_age - u_age_threshold);
    if (random(a_position + vec2(u_time * 0.1, age * 0.01)) < resetProbability) {
        shouldReset = true;
    }
}

if (age > u_max_age) {
    shouldReset = true;
}

There is a setSource method that provides an optional parameter percentParticleWhenSetSource which default to 0.5. This parameter gives an extra opportunity to reset particles positions when the source url (e.g., wind pattern at a different timestamp) is updated, so a new wind particle distribution pattern can be emerged more quickly.

Lastly, the initial version of the package generates random coordinates in the range [0,1] when shouldReset is true. I do code a generateBoundaryPosition function in the source code, but it is not used since particles might easily go out of bounds in certain wind pattern setups.

3. Making particle position update GPU-based

Whether it is mapbox-gl-js library or WebGL code, it relies on GPU for various computation. To maximize the computational speed to allow a smooth user experience, the data transfers between GPU and CPU should be minimized. Therefore, my implementation relies on transform feedback to make particle position update fully GPU-based.

The WebGL code snippet that I provided in previous two sections are inside a dedicated vertex shader for position update (there is another vertex shader for actual mapping). This position update vertex shader takes current position and age buffers as inputs, and outputs updated position and age buffers. Transform feedback setup is needed to capture the vertex shader outputs and use them as inputs for the next frame.

First, two sets of buffers for positions and two sets of buffers for ages are needed (those code are inside onAdd method of Mapbox custom layer).

this.particleBufferA = createBuffer(gl, positions);
this.particleBufferB = createBuffer(gl, positions);
this.currentBuffer = this.particleBufferA;
this.nextBuffer = this.particleBufferB;

this.ageBufferA = createBuffer(gl, ages);
this.ageBufferB = createBuffer(gl, ages);
this.currentAgeBuffer = this.ageBufferA;
this.nextAgeBuffer = this.ageBufferB;

Next, transform feedback object needs to be created.

this.transformFeedback = gl.createTransformFeedback();

In Mapbox custom layer render method, output buffers need to be bind to the transform feedback.

gl.useProgram(this.updateProgram.program);

gl.bindTransformFeedback(gl.TRANSFORM_FEEDBACK, this.transformFeedback);
gl.bindBufferBase(gl.TRANSFORM_FEEDBACK_BUFFER, 0, this.nextBuffer);
gl.bindBufferBase(gl.TRANSFORM_FEEDBACK_BUFFER, 1, this.nextAgeBuffer);

It is different from how input buffers are bind:

gl.bindBuffer(gl.ARRAY_BUFFER, this.currentBuffer);
gl.enableVertexAttribArray(this.updateProgram.a_position);
gl.vertexAttribPointer(this.updateProgram.a_position, 2, gl.FLOAT, false, 0, 0);

gl.bindBuffer(gl.ARRAY_BUFFER, this.currentAgeBuffer);
gl.enableVertexAttribArray(this.updateProgram.a_age);
gl.vertexAttribPointer(this.updateProgram.a_age, 1, gl.FLOAT, false, 0, 0);

Now, it is possible to begin transform feedback, run update position vertex shader, and do any necessary cleanup

gl.beginTransformFeedback(gl.POINTS);

gl.drawArrays(gl.POINTS, 0, this.particleCount);

gl.endTransformFeedback();

gl.bindTransformFeedback(gl.TRANSFORM_FEEDBACK, null);
gl.bindBufferBase(gl.TRANSFORM_FEEDBACK_BUFFER, 0, null);
gl.bindBufferBase(gl.TRANSFORM_FEEDBACK_BUFFER, 1, null);

Lastly, because each frame’s output becomes next frame’s input, buffers need to be swapped:

[this.currentBuffer, this.nextBuffer] = [this.nextBuffer, this.currentBuffer];
[this.currentAgeBuffer, this.nextAgeBuffer] = [this.nextAgeBuffer, this.currentAgeBuffer];

Also note that here updateProgram is used. A WebGL program consists of a vertex shader and a fragment shader. The vertex shader of updateProgram is already discussed in the previous sections. The fragment shader in updateProgram does not need to output anything visual because we are just trying to capture new positions and ages; it serves as a placeholder here.

4. Mapping the positions of particles

So far I focus on introducing the updateProgram and especially its vertex shader. Now I will introduce rendering step, in which renderProgram is used. The vertex shader of renderProgram is responsible for mapping the new positions of particles.

The particle position buffer assumes a range [0, 1] relative to the layer extent, so it is necessary to convert it to actual longitude/latitude for further processing so it can eventually overlay with other Mapbox layers.

float lng = mix(u_bounds[0], u_bounds[2], currentPos.x);
float lat = mix(u_bounds[3], u_bounds[1], 1.0 - currentPos.y);

Mapbox has a matrix transformation that takes a web mercator [0, 1] normalized coordinate pair and outputs position used for WebGL rendering (see its example). This will be the “further processing” that is needed.

vec2 mercator = latLngToMercator(vec2(lng, lat));
gl_Position = u_matrix * vec4(mercator, 0, 1);

5. Coloring the particles

The fragment shader of renderProgram is responsible for assigning colors to the particles.

In fact, as a pre-requisite, a colormap needs to be created and stored as a WebGL texture. This involves using the user input parameter color which is an array like [ [Wind speed in mph, [R, G, B]] …] and min-speed,max-speed (denoted as u_speed_range in shader program) stored in EXIF string.

The colormap is created by first sorting the color array by wind speed. Then it creates 256 discrete stops within the range bounded by min-speed and max-speed. For each stop (a speed value), it identifies its corresponding interval in the color array (e.g., between wind speed 10 mph and 12 mph), and performs an interpolation on R, G, B values (e.g., stop wind speed is 11 mph, 10 mph has RGB values [0, 0, 0], 12 mph has RGB values [100, 100, 100], the interpolated RGB values for this stop will be [50, 50, 50]).

In the fragment shader, the normalized u- and v-component velocity needs to be retrieved from the velocity texture, and it needs to be scaled to original min/max range. A speed value then can be calculated.

vec4 velocity = texture(u_velocity_texture, v_position);
float u = mix(u_value_range_u[0], u_value_range_u[1], velocity.r);
float v = mix(u_value_range_v[0], u_value_range_v[1], velocity.g);
float speed = length(vec2(u, v));

Next we need to sample the speed value in color texture. Sampling value from texture will require a normalized value between 0 and 1. Therefore, we need to normalize the speed to [0, 1] range by taking min-speed and max-speed that we use to construct the colormap into account.

float normalizedSpeed = (speed - u_speed_range[0]) / (u_speed_range[1] - u_speed_range[0]);
vec4 color = texture(u_wind_color, vec2(normalizedSpeed, 0.5));

This is the general idea of how coloring of particles work. There are some additional processing in the fragment shader, such as making points’ shape from square to circular, or adding global opacity value; I will omit those details here.

6. Adding trail offset (a.k.a., tails)

To draw the main particle and its trail, I create buffer like following:

const trailOffsets = new Float32Array(this.trailLength + 1);
for (let i = 0; i <= this.trailLength; i++) {
    trailOffsets[i] = i; // 0 = main particle, 1,2,3... = trail segments
}
this.trailOffsetBuffer = createBuffer(gl, trailOffsets);

In the render method of Mapbox custom layer, each main particle will be drawn (trailLength+1) times with different offsets.

gl.drawArraysInstanced(gl.POINTS, 0, this.particleCount, this.trailLength + 1);

In the vertex shader of the renderProgram, besides the updated particle position as input, the trail offset (0, 1, 2, etc.) is also an input. If trail offset is non-zero (i.e., other than the main particle), the trail position will be computed by moving backwards from main particle along its velocity path.

if (trailOffset > 0.0) {
    currentPos = mainPos - velocity * u_speed_factor * trailOffset * 1.5;
}

The vertex shader also controls the size of the point. As the trail position increases, the size of the point decreases.

float size = trailOffset == 0.0 ? u_point_size : u_point_size * pow(u_trail_size_decay, trailOffset);
gl_PointSize = size;

That is it! I give a technical walk-through of how I implement mapbox-exif-layer package for wind visualization. The package can also add what I call smooth-raster layers that can be used to display various weather layers such as temperature.

I hope you learn something from it or or find it to be useful. Any community involvements in further refining this package will be much appreciated!