What Does Wind Look Like?

Jordi Aranda Procesamiento de datos

In our recent collaboration with United Nations Global Pulse to measure the economic impact of natural disasters we analyzed purchase behaviors when Hurricane Odile struck Baja California Sur in September 2014. As part of the study we visualized the strength of the winds and showed how the hurricane initially formed until finally making landfall.

In this article I describe how we produced this visualization by gathering wind data and animate it on the web.

Data sources

A great resource for wind data is NOAA, which offers plenty of data services, including GFS. GFS is a weather forecast model produced by the National Centers for Environmental Prediction, with dozens of atmospheric and land-soil variables available.

This forecast model covers the entire globe, forming a grid of points (at different scales) which is used by the operational forecasters to predict weather out to 16 days in the future, at different hourly intervals. We extracted historical data on three variables: the date range (September 2014) and the wind vector components (horizontal and vertical) of the desired geographic extension to cover (Mexico). Aside from that, time granularity of the weather forecasts was important, so we chose the maximum avaialable, 3-hourly intervals. With these variables, we were good to go.

Data wrangling

Data extraction

GFS offers historical data and weather forecast data in a specific file format called GRIB, which is a binary file format widely adopted in meteorology to store historical data and forecast weather data. In order to process these datasets we had to install wgrib2, a command-line tool used to query and extract data from grib files.

We extracted the following data:

  •  U-component (horizontal, UGRD) of wind vector (10m above ground)
  •  V-component (vertical, VGRD) of wind vector (10m above ground)
  • Temperature (TMP) (100m above ground)

We used the following variables to query the datasets:

  •  Dates range between 10 September 2014 and 18 September 2014 (hurricane Odile made landfall on 15th September 2014)
    • Bounding box of covered region (in longitude/latitude coordinates): *
      • Bottom: 14.5
      • Left: -118.5 *
      • Top: 32.5
      • Right: -86.5
    • 3-hourly interval forecasts (maximum time granularity available)
    • 0.5 grid scale (measures taken every 0.5 degrees in longitude/latitude space)

Here you can find the script we used to extract the required data in csv format.

This is what the extracted dataset looked like:


$ head /tmp/result.csv
"t0","t1","variable","filter","lon","lat","value"
"2014-09-10 00:00:00","2014-09-10 03:00:00","UGRD","10 m above ground",-118.5,14.5,-3.48
"2014-09-10 00:00:00","2014-09-10 03:00:00","UGRD","10 m above ground",-118,14.5,-3.55
"2014-09-10 00:00:00","2014-09-10 03:00:00","UGRD","10 m above ground",-117.5,14.5,-4.96
"2014-09-10 00:00:00","2014-09-10 03:00:00","UGRD","10 m above ground",-117,14.5,-6.28
"2014-09-10 00:00:00","2014-09-10 03:00:00","UGRD","10 m above ground",-116.5,14.5,-6.82
"2014-09-10 00:00:00","2014-09-10 03:00:00","UGRD","10 m above ground",-116,14.5,-6.23
"2014-09-10 00:00:00","2014-09-10 03:00:00","UGRD","10 m above ground",-115.5,14.5,-6.88
"2014-09-10 00:00:00","2014-09-10 03:00:00","UGRD","10 m above ground",-115,14.5,-7.53
"2014-09-10 00:00:00","2014-09-10 03:00:00","UGRD","10 m above ground",-114.5,14.5,-6.69

The extracted fields are fairly self-explanatory:

  • t0: Timestamp of when the forecast was made
  • t1: Timestamp of the forecast
  • variable: name of measure extracted
  • filter: height of measure
  • lon: longitude of measure
  • lat: latitude of measure
  • value: value of measure

Data processing

At this stage we had a grid of points with wind and temperature measures available every 3 hours for the specified dates range and bounding box. Since hurricanes evolve at a very high pace, after some tests we noticed this time granularity was not enough to reproduce the hurricane trajectory properly. On the other hand, the grid was not dense enough either (measures were only available every 0.5 degrees in longitude/latitude coordinates) so we end up interpolating values for both space and time, so that we could rely on more points over time to draw the wind.

We used mostly pandas to carry out this grid interpolation. This script does pretty much all the necessary work and generates a new json file (data.json) containing grid interpolated measures of wind data (horizontal/vertical components, wind speed) in an hourly basis, which was meant to be consumed later on by the visualization.

Data visualization

The input for the data visualization was a list of grids, every one corresponding to the wind vectors for an specific hour in a day. Every grid could actually be seen as a vector field.

Vector fields are often used to model speed and direction of moving fluids, or the strength and direction of some force, for example. In our case, our whole purpose was to model the wind speed and direction.

In order to visualize how the wind moves over time we used the JavaScript canvas API, simulating a particle system driven by the wind vector field at some point in time t. Every some number of frames, the grid of wind vectors changes dynamically and so do the particles, achieving the effect of how the hurricane formed and moved over time.

Since our grid of wind measures was finite, we had to interpolate particles positions using their closest grid points as well.

Building the particle system

Particles can be seen as a list of objects with the following properties:

  •  x: current x canvas position
  • y: current y canvas position
  • startX: initial x canvas position
  • startY: initial y canvas position
  •  lastX: last x canvas position
  • lastY: last y canvas position
  • age: the particle age, used to reset particles’ positions after some time

The following piece of code initializes the canvas and the particle system:


const canvas = document.getElementById('canvas');
const ctx = canvas.getContext('2d');
const width = canvas.width;
const height = canvas.height;

// Particles related
const N_PARTICLES = 2000;
const PARTICLE_MAX_AGE = 30;
let particles = _.map(_.range(N_PARTICLES), () => {
  let x = Math.random() * width;
  let y = Math.random() * height;
  return {
      x: x,
      y: y,
      startX: x,
      startY: y,
      lastX: x,
      lastY: y,
      age: 0
  };
});

Fast lookups

If we have a look at our dataset, we notice how data is sorted by date and longitude values increase first, and only then latitude values do so. We will take this into account to index our data longitude-wise, creating a multi-index data structure. This is necessary since we will be querying wind vectors values quite often and these operations need to be fast.

The following function returns the corresponding longitude and latitude indexes of our data structure given the longitude and latitude coordinates (in degrees):


function lonLatToGridCoords(lon, lat, bounds, rows, cols, deltaX, deltaY) {
  let x = Math.floor1)lon - bounds.left)/deltaX);
  if (x == cols) x -= 1;
  let y = Math.floor((bounds.top - lat)/deltaY);
  if (y == rows) y -= 1;
  return [x,y];
};


This is used when building the multi-index grid later on:

function createGrid(gridData) { const points = gridData.data.points; const uv = gridData.data.uv; const values = gridData.data.values; if (!(points.length == uv.length && points.length == values.length throw new Error('Points, UV, Values vectors must be of same length'); // Bounding box bounds const bounds = { left: d3.min(_.map(points, d => d[0])), bottom: d3.min(_.map(points, d => d[1])), right: d3.max(_.map(points, d => d[0])), top: d3.max(_.map(points, d => d[1])) }; const extentX = d3.extent(points, d => d[0]); const extentY = d3.extent(points, d => d[1]); const cols = (extentX[1] - extentX[0])/0.5 + 1; const rows = (extentY[1] - extentY[0])/0.5 + 1; // Allocate grid (indexation is row-wise => x=normalized longitude, y=normalized latitude) let grid = _.map(_.range(cols), () => { return _.map(_.range(rows), () => { return {v: [0, 0]}; }); }); // Populate grid with proper values _.forEach(points, (point, idx) => { // Find grid index based on point lon/lat values const pos = lonLatToGridCoords(point[0], point[1], bounds, rows, cols, 0.5, 0.5); // Assign wind vector to this position grid[pos[0]][pos[1]].v[0] = uv[idx][0]; grid[pos[0]][pos[1]].v[1] = -uv[idx][1]; // Check max wind speed for wind color let windSpeed = Math.sqrt(Math.pow(uv[idx][0], 2) + Math.pow(uv[idx][1], 2)); if (windSpeed > MAX_WIND) MAX_WIND = windSpeed; if (windSpeed < MIN_WIND) MIN_WIND = windSpeed; }); return grid; }; // Hourly grids let tGrids = &#91;&#93; // Create grid (`gridData` corresponds to the JSON object created during the data processing step) _.forEach(gridData, value => { tGrids.push(createGrid(value)); }); // Set current grid let currentGrid = tGrids[0];

We are now good to go to implement the draw function, which is what canvas executes at the specified frame rate and is in charge of the wind drawing.

The draw loop

The canvas loop function has different purposes, in general terms:

  • To update the particles positions interpolating the values from the grid multi-index
  • To redraw the particles
  • To reset particles positions if they they are too “old”, meaning they should disappear for aesthetics

At some point, every particle will occupy an specific position in the canvas. This position can be mapped directly to an specific point within the grid of wind vectors and is used to retrieve the four closest wind vectors to such point. To obtain the corresponding wind vector at such position we do interpolate linearly (“lerping”) between these vectors. The resulting vector tells us how this particle should move, in which direction and with what speed.

Since on every draw call we will be drawing the wind particles in the updated positions we need some mechanism to “clean” old wind positions. In this scenario, and in order to acoomplish the wind trail effect, a common practice is to fill the canvas with some low opacity, so that after some time, old particle positions are wiped out.

The following code snippet is in charge of all this work:


// Some variables to control the wind aesthetics and behaviour
const PARTICLE_MAX_AGE = 30;
const LOOP_PER_GRID = 5;
const DELTA_SPEED = 0.001;
let currentFrame = 0;

function draw() {

  const COLS = _.size(tGrids[0]);
  const ROWS = _.size(tGrids[0][0]);
  const X_OFFSET = width/COLS;
  const Y_OFFSET = height/ROWS;

  ctx.save();

  // Apply opacity to accomplish the trail effect
  ctx.fillStyle = 'rgba(0, 0, 0, .05)';
  ctx.fillRect(0, 0, width, height);

  // Draw particles
  _.forEach(particles, (particle, idx) => {

    // Increase particle age
    particles[idx].age += Math.random()*2;

    // Check if particle is in bounds
    if (inBounds(particle, width, height)) {

      // If particle exceeded its age, should be placed somewhere else and start from scratch again
      if (particles[idx].age > PARTICLE_MAX_AGE) {
        let newX = Math.random() * width;
        let newY = Math.random() * height;
        particles[idx].x = newX;
        particles[idx].y = newY;
        particles[idx].startX = newX;
        particles[idx].startY = newY;
        particles[idx].age = 0;
      }

      // Compute new particle position
      let i = Math.floor(particle.x / X_OFFSET);
      let j = Math.floor(particle.y / Y_OFFSET);
      if (i < 0 || i > (COLS - 1)) throw new Error('Invalid i-index for grid' + i);
      if (j < 0 || j > (ROWS - 1)) throw new Error('Invalid j-index for grid' + j);
      let a, b, c, d, u, v, colA, colB;

      if (i < (COLS - 1)) {
        colA = currentGrid&#91;i&#93;;
        colB = currentGrid&#91;i + 1&#93;;
      } else {
        colA = currentGrid&#91;i - 1&#93;;
        colB = currentGrid&#91;i&#93;;
      }
      if (j < (ROWS - 1)) {
        a = colA&#91;j&#93;;
        b = colB&#91;j&#93;;
        c = colB&#91;j + 1&#93;;
        d = colA&#91;j + 1&#93;;
      } else {
        a = colA&#91;j - 1&#93;;
        b = colB&#91;j - 1&#93;;
        c = colB&#91;j&#93;;
        d = colA&#91;j&#93;;
      }
      if (!a.hasOwnProperty('v') ||
          !b.hasOwnProperty('v') ||
          !c.hasOwnProperty('v') ||
          !d.hasOwnProperty('v')) throw new Error('Could not locate related grid vectors');
      u = (particle.x % X_OFFSET) / X_OFFSET;
      v = (particle.y % Y_OFFSET) / Y_OFFSET;
      if (u < 0.0 || u > 1.0) throw new Error('Invalid value for u component');
      if (v < 0.0 || v > 1.0) throw new Error('Invalid value for v component');
      let interpolatedVector = quadLerp(a.v, b.v, c.v, d.v, u, v);

      // Draw particle in current position
      ctx.beginPath();
      let windSpeed = Math.sqrt(Math.pow(interpolatedVector[0], 2) + Math.pow(interpolatedVector[1], 2));
      let windRatio = windSpeed/MAX_WIND;
      let color = chroma(colorScale(windRatio)).alpha(1).css();
      ctx.strokeStyle = color;
      ctx.moveTo(particle.x, particle.y);
      ctx.lineTo(particle.x + interpolatedVector[0]*windRatio, particle.y + interpolatedVector[1]*windRatio);
      ctx.stroke();

      // Update positions
      particles[idx].lastX = particles[idx].x;
      particles[idx].lastY = particles[idx].y;
      particles[idx].x += (interpolatedVector[0]*windRatio);
      particles[idx].y += (interpolatedVector[1]*windRatio);

      // Check if particle got stuck (diff in x/y velocities is too small)
      if (Math.sqrt2)Math.pow(interpolatedVector[0], 2) + Math.pow(interpolatedVector[1], 2) < DELTA_SPEED){
        let newX = Math.random() * width;
        let newY = Math.random() * height;
        particles&#91;idx&#93;.x = newX;
        particles&#91;idx&#93;.y = newY;
        particles&#91;idx&#93;.startX = newX;
        particles&#91;idx&#93;.startY = newY;
        particles&#91;idx&#93;.age = 0;
      }
    } else {
      let newX = Math.random() * width;
      let newY = Math.random() * height;
      particles&#91;idx&#93;.x = newX;
      particles&#91;idx&#93;.y = newY;
      particles&#91;idx&#93;.startX = newX;
      particles&#91;idx&#93;.startY = newY;
      particles&#91;idx&#93;.age = 0;
    }
  });

  if (++currentFrame%LOOP_PER_GRID == 0) {
    let idx = Math.floor(currentFrame/LOOP_PER_GRID);
    if (idx == tGrids.length) {
      currentFrame = 0;
      idx = 0;
    }
    currentGrid = tGrids&#91;idx&#93;;
  }

  ctx.restore();
}

&#91;/code&#93;

We can finally call the <code>draw</code> function at the desired frame rate:



let intervalId = null;
let requestId = null;
let fps = 30;

function loop() {
  intervalId = setTimeout(() => {
    requestId = requestAnimationFrame(loop);
    draw();
  }, 1000/fps)
}

loop();

Final result

This is the final result (wind layer overlaps the Mexico map layer).

Read the whole story on how we measured people’s economic resilience to Hurrican Odile.

 

 

References

References
1 lon - bounds.left)/deltaX); if (x == cols) x -= 1; let y = Math.floor((bounds.top - lat)/deltaY); if (y == rows) y -= 1; return [x,y]; };

This is used when building the multi-index grid later on:


function createGrid(gridData) {
  const points = gridData.data.points;
  const uv = gridData.data.uv;
  const values = gridData.data.values;
  if (!(points.length == uv.length && points.length == values.length
2 Math.pow(interpolatedVector[0], 2) + Math.pow(interpolatedVector[1], 2