Boosting Forecast Accuracy: Bilinear Interpolation in Action

Tobias SteinhagenTobias Steinhagen
13 min read

Where we left off

In the last part, we installed our primary tool, grib2json, and learned how to obtain data from it. We also examined the data and explored the most important properties of the header. Additionally, we discussed the need to convert our location into an array index to obtain the forecast values from the data. If you missed any of this, please refer back to Part 3.

Stuff we do now

In this section, we aim to develop the functions for converting the location into the array index, as well as interpolating the forecast values if the location falls between grid points. We'll also create a few helper functions. In the final step, we will save the forecast data in MongoDB.

Unwind with Forecast Object: Simplify

Let's start by streamlining the process of working with all the essential values from the .grib header and the filename. We will consolidate them into a single object. To achieve this, we will create a function that takes the forecast header object as the first argument and destructures it. The second argument will be the filename. We will employ regular expressions to extract the forecast type and forecast name. All other values from the header will be returned with the same names as in the header to prevent confusion over time.

const getforecastHeader = (
  { lo1, lo2, la1, la2, dx, dy, refTime, forecastTime },
  filename,
) => {
  const regex = /(?<=_)[0-9]+_[0-9]+_[0-9]+_[A-Za-z]+(?=.grib)/;
  const regex2 = /(?<=\/)[A-Za-z]+(?=_)/;
  const forecastType = filename.match(regex)[0].split('_')[3];
  const forecastName = filename.match(regex2)[0];
  return {
    forecastName,
    forecastType,
    refTime,
    forecastTime,
    lo1,
    lo2,
    la1,
    la2,
    dx,
    dy,
  };
};

Forecast Fusion: Spot-On Savings

Now, after consolidating all our forecast information into one object, we should proceed to save it with the corresponding forecast values in the spot and forecast models.

The easiest approach is to iterate over the spots array and save the values one by one. First, we need to check if the spot falls within the boundaries of our forecast model. If not, we can easily skip the spot.

The next step is to calculate the forecast value for the spot. Since the spots are most likely not exactly aligned with the grid, we have to interpolate the values for greater accuracy.

To improve the performance of this step, we can map all these calculations into a single array and compute all spots in parallel. By using Promise.all() with that array, we can wait for all promises to resolve and store the resulting array of forecast values in a variable.

The final step is to loop over the spots and update the forecasts. After completing this step, we are ready to save the forecastInfo to MongoDB.

// src/convert_grib/index.js
...
const populateSpots = async (filename, spots) => {
  const forecastJson = await getJson(filename, {
    scriptPath: './src/convert_grib/grib2json/src/bin/grib2json',
    names: true, // (default false): Return descriptive names too
    data: true, // (default false): Return data, not just headers
  });

  // get grib info from header
  const forecastHeader = getforecastHeader(forecastJson[0].header, filename);
  // get forecastInfo document from db or create new one
  const forecastInfo = await getForecastInfo(forecastHeader);

  // Calculate data values for all spots in parallel
  const dataValuePromises = spots.map((spot) => {
    // check if spot is in model borders
    if (inGrid(spot, forecastHeader)) {
      return calculateDataValue(spot, forecastHeader, forecastJson[0].data);
    }
    return null;
  });

  // Wait for all dataValue promises to resolve
  const dataValues = await Promise.all(dataValuePromises);
  // Update spot forecasts with calculated data values
  spots.forEach((spot, index) => {
    if (dataValues[index] !== null) {
      updateSpotForecast(spot, forecastInfo, forecastHeader, dataValues[index]);
    }
  });
  await forecastInfo.save();
};
...

Getting the ForecastInfo Document

When updating the forecasts, it is also necessary to update the forecastInfo document, as preferences may change over time. For this purpose, I created the getForecastInfo() function. This function takes the forecastHeader object as its only argument. The first step is to destructure the forecastHeader, after which we make a findOneAndUpdate query to search for an existing ForecastInfo Document in our database. If we don't find anything, we want to create one. By adding the option object with upsert: 'true', we specify that we want to create a new document if no document is found. By adding new: 'true', we receive the updated document instead of the found one. The last step is to return the forecastInfo.

// src/convert_grib/index.js
...
const getForecastInfo = async (forecastHeader) => {
  const { forecastName, refTime, lo1, lo2, la1, la2, dy, dx } = forecastHeader;

  // check if forecast already exists
  // if not create new forecast
  const forecastInfo = await ForecastInfo.findOneAndUpdate(
    { name: forecastName },
    {
      name: forecastName,
      time: refTime,
      lo1,
      lo2,
      la1,
      la2,
      dy,
      dx,
    },
    { upsert: true, new: true },
  );

  return forecastInfo;
};
...

Validating Spot Coordinates

First, let's take a look at the inGrid() function. We require this function to check whether our spot location is covered by the current forecast model we are using. To perform this check, we simply need to verify that our spot's latitude and longitude values fall between the start and end latitude and longitude values of the forecast model.

One important point to keep in mind is that our longitude value ranges from 0 to 360, and right after 360, it starts again at 0. This behavior is because our Earth is a sphere and not flat. To simplify dealing with this behavior, we just add 360 if our start point is greater than our endpoint.

After making this adjustment, the remainder of the process involves simple comparisons. To improve code readability and minimize repetition, I have created two helper functions. Ultimately, the inGrid() function returns a boolean value.

// src/convert_grib/index.js
...
const getAbsoluteLon = (lonStart, lonEnd) => {
  return lonStart > lonEnd ? lonEnd + 360 : lonEnd;
};

const isBetween = (x, min, max) => {
  return x >= min && x <= max;
};

const inGrid = (spot, forecastHeader) => {
  const lo2 = getAbsoluteLon(forecastHeader.lo1, forecastHeader.lo2);
  const spotLon = getAbsoluteLon(forecastHeader.lo1, spot.lon);

  return (
    isBetween(spot.lat, forecastHeader.la1, forecastHeader.la2) &&
    isBetween(spotLon, forecastHeader.lo1, lo2)
  );
};
...

Interpolation Showdown: Picking the Right Method!

When it comes to interpolation, we have different methods available: 2D-nearest-neighbor, bilinear, and bicubic interpolation. 2D-nearest-neighbor interpolation requires just two points, while bilinear interpolation requires 4 points, and bicubic interpolation requires 16 points to work properly.

I will be writing an article comparing these three interpolation methods for weather forecasts in the future. However, for now, considering the effort and benefit, I have chosen bilinear interpolation because it provides more accurate values than 2D-nearest-neighbor interpolation and requires only 4 points to work. On the other hand, bicubic interpolation may yield slightly better results, but it requires 4 times the points.

In the end, it all comes down to the specific requirements of your application.

Calculate the Data value with bilinear interpolation

Next, we call the calculateDataValue() function, which performs bilinear interpolation. As arguments, we pass our spot object containing the coordinates of the spot, forecastHeader with values like the grid's mesh size, and forecastData, which contains the raw forecast values. The calculation closely follows the Wikipedia article about bilinear interpolation. For easy reference, I retained the same variable names.

To facilitate the process, I've created three helper functions. The first two functions, getMinPoint() and getMaxPoint(), fetches the four points surrounding our desired spot location, while the third one, the getGribIndex() function, converts these locations into a 1D-array index.

In the getGribIndex() function, we calculate the row and position (column) of our value within the forecast grid. By multiplying the number of rows with the row's width and adding the position of the value, we obtain the 1D-array index.

// src/convert_grib/index.js
...
const getMinPoint = (point, delta) => {
  return point % delta === 0 ? point : point - (point % delta);
};

const getMaxPoint = (point, delta) => {
  return point % delta === 0 ? point : point - (point % delta) + delta;
};

const getGribIndex = (forecastHeader, spot) => {
  // check if end value for longitute is lower than start value
  const lo2 = getAbsoluteLon(forecastHeader.lo1, forecastHeader.lo2);
  const spotLon = getAbsoluteLon(forecastHeader.lo1, spot.lon);

  const latRow = (spot.lat - forecastHeader.la1) / forecastHeader.dy;
  const latWidth = (lo2 - forecastHeader.lo1) / forecastHeader.dx + 1;
  const lonPos = (spotLon - forecastHeader.lo1) / forecastHeader.dx;
  return Math.round(latRow * latWidth + lonPos);
};

const calculateDataValue = (spot, forecastHeader, forecastData) => {
  // bilinear interpolation for 4 points around spot position
  // https://en.wikipedia.org/wiki/Bilinear_interpolation
  const x = getAbsoluteLon(forecastHeader.lo1, spot.lon);
  const y = spot.lat;
  const x1 = getMinPoint(x, forecastHeader.dx);
  const x2 = getMaxPoint(x, forecastHeader.dx);
  const y1 = getMinPoint(y, forecastHeader.dy);
  const y2 = getMaxPoint(y, forecastHeader.dy);

  const Q11 = forecastData[getGribIndex(forecastHeader, { lon: x1, lat: y1 })];
  const Q21 = forecastData[getGribIndex(forecastHeader, { lon: x2, lat: y1 })];
  const Q22 = forecastData[getGribIndex(forecastHeader, { lon: x2, lat: y2 })];
  const Q12 = forecastData[getGribIndex(forecastHeader, { lon: x1, lat: y2 })];
  const R1 = ((x2 - x) / (x2 - x1)) * Q11 + ((x - x1) / (x2 - x1)) * Q21;
  const R2 = ((x2 - x) / (x2 - x1)) * Q12 + ((x - x1) / (x2 - x1)) * Q22;
  const P = ((y2 - y) / (y2 - y1)) * R1 + ((y - y1) / (y2 - y1)) * R2;

  return P;
};
...

Update the spots forecasts

Updating the forecast document is relatively simple, but there are a few things to consider. We invoke the updateSpotForecast() function with the following arguments: the spot document, the forecastInfo document, the forecastHeader object, and the dataValue obtained from the bilinear interpolation.

First, we need to determine if a forecast document for this forecast model already exists in the spot.forecasts array. If no forecast document is present in the array, we must create one by calling new Forecast() and push the newly created forecast into the spot.forecasts array. Conversely, if we find a forecast document, we can simply update the existing forecast document with the timestamp and the data value for the corresponding forecast time.

Finally, we save the forecast document and repopulate the spot.forecasts array to include the newly created forecast.

// src/convert_grib/index.js
...
const updateSpotForecast = async (
  spot,
  forecastInfo,
  forecastHeader,
  dataValue,
) => {
  // check if forecast already exists
  const forecastFound = spot.forecasts.find(
    (spotForecast) =>
      spotForecast.forecastInfo.toString() === forecastInfo._id.toString(),
  );

  // if not create new forecast
  if (!forecastFound) {
    const forecastData = new Forecast({
      _id: new mongoose.Types.ObjectId(),
      forecastInfo,
      time: forecastHeader.refTime,
      [forecastHeader.forecastType]: {
        [forecastHeader.forecastTime]: dataValue,
      },
    });
    spot.forecasts.push(forecastData);

    await forecastData.save();
  } else {
    // if forecast exists update data
    forecastFound[forecastHeader.forecastType] = {
      ...forecastFound[forecastHeader.forecastType],
      [forecastHeader.forecastTime]: dataValue,
    };
    await forecastFound.save();
  }
  await spot.populate({
    path: 'forecasts',
    match: { forecastInfo: forecastInfo._id },
  });
};
...

back to convertGrib()

That was the last part of our convertGrib() function. Now we have everything together to run the function and update the spot, forecast and forecastInfo model in MongoDB.

Conclusion

In conclusion, our weather forecast system has made significant strides in accuracy and efficiency. We have streamlined the process, consolidated forecast values, and successfully extracted data using grib2json.

The implementation of bilinear interpolation has provided precise forecast values even for spots between grid points. This method strikes the right balance between accuracy and performance, making it a valuable choice for our system.

We have updated spot forecasts in MongoDB, ensuring that preferences and changes are well-managed over time.

In the next part, we will stitch together everything from downloading the files to saving the forecast values in one single updateForecast function. We will also implement the npm package cron. This will allow us to schedule the execution of our updateForecast function at specific times without manually calling the node app over and over again.

Thank you so much for sticking with me! As always, I genuinely appreciate your comments and feedback. I'll see you all in Part 5.

-Tobias

CODE:

// src/convert_grib/index.js
const mongoose = require('mongoose');
const util = require('util');
const grib2json = require('grib2json').default;
const { Spot, Forecast, ForecastInfo } = require('../models');

const getJson = util.promisify(grib2json);

const getforecastHeader = (
  { lo1, lo2, la1, la2, dx, dy, refTime, forecastTime },
  filename,
) => {
  const regex = /(?<=_)[0-9]+_[0-9]+_[0-9]+_[A-Za-z]+(?=.grib)/;
  const regex2 = /(?<=\/)[A-Za-z]+(?=_)/;
  const forecastType = filename.match(regex)[0].split('_')[3];
  const forecastName = filename.match(regex2)[0];
  return {
    forecastName,
    forecastType,
    refTime,
    forecastTime,
    lo1,
    lo2,
    la1,
    la2,
    dx,
    dy,
  };
};

const getForecastInfo = async (forecastHeader) => {
  const { forecastName, refTime, lo1, lo2, la1, la2, dy, dx } = forecastHeader;

  // check if forecast already exists
  // if not create new forecast
  const forecastInfo = await ForecastInfo.findOneAndUpdate(
    { name: forecastName },
    {
      name: forecastName,
      time: refTime,
      lo1,
      lo2,
      la1,
      la2,
      dy,
      dx,
    },
    { upsert: true, new: true },
  );

  return forecastInfo;
};

const getAbsoluteLon = (lonStart, lonEnd) => {
  return lonStart > lonEnd ? lonEnd + 360 : lonEnd;
};

const isBetween = (x, min, max) => {
  return x >= min && x <= max;
};

const inGrid = (spot, forecastHeader) => {
  const lo2 = getAbsoluteLon(forecastHeader.lo1, forecastHeader.lo2);
  const spotLon = getAbsoluteLon(forecastHeader.lo1, spot.lon);

  return (
    isBetween(spot.lat, forecastHeader.la1, forecastHeader.la2) &&
    isBetween(spotLon, forecastHeader.lo1, lo2)
  );
};

const getMinPoint = (point, delta) => {
  return point % delta === 0 ? point : point - (point % delta);
};

const getMaxPoint = (point, delta) => {
  return point % delta === 0 ? point : point - (point % delta) + delta;
};

const getGribIndex = (forecastHeader, spot) => {
  // check if end value for longitute is lower than start value
  const lo2 = getAbsoluteLon(forecastHeader.lo1, forecastHeader.lo2);
  const spotLon = getAbsoluteLon(forecastHeader.lo1, spot.lon);

  const latRow = (spot.lat - forecastHeader.la1) / forecastHeader.dy;
  const latWidth = (lo2 - forecastHeader.lo1) / forecastHeader.dx + 1;
  const lonPos = (spotLon - forecastHeader.lo1) / forecastHeader.dx;
  return Math.round(latRow * latWidth + lonPos);
};

const calculateDataValue = (spot, forecastHeader, forecastData) => {
  // bilinear interpolation for 4 points around spot position
  // https://en.wikipedia.org/wiki/Bilinear_interpolation
  const x = getAbsoluteLon(forecastHeader.lo1, spot.lon);
  const y = spot.lat;
  const x1 = getMinPoint(x, forecastHeader.dx);
  const x2 = getMaxPoint(x, forecastHeader.dx);
  const y1 = getMinPoint(y, forecastHeader.dy);
  const y2 = getMaxPoint(y, forecastHeader.dy);

  const Q11 = forecastData[getGribIndex(forecastHeader, { lon: x1, lat: y1 })];
  const Q21 = forecastData[getGribIndex(forecastHeader, { lon: x2, lat: y1 })];
  const Q22 = forecastData[getGribIndex(forecastHeader, { lon: x2, lat: y2 })];
  const Q12 = forecastData[getGribIndex(forecastHeader, { lon: x1, lat: y2 })];
  const R1 = ((x2 - x) / (x2 - x1)) * Q11 + ((x - x1) / (x2 - x1)) * Q21;
  const R2 = ((x2 - x) / (x2 - x1)) * Q12 + ((x - x1) / (x2 - x1)) * Q22;
  const P = ((y2 - y) / (y2 - y1)) * R1 + ((y - y1) / (y2 - y1)) * R2;

  return P;
};

const updateSpotForecast = async (
  spot,
  forecastInfo,
  forecastHeader,
  dataValue,
) => {
  // check if forecast already exists
  const forecastFound = spot.forecasts.find(
    (spotForecast) =>
      spotForecast.forecastInfo.toString() === forecastInfo._id.toString(),
  );

  // if not create new forecast
  if (!forecastFound) {
    const forecastData = new Forecast({
      _id: new mongoose.Types.ObjectId(),
      forecastInfo,
      time: forecastHeader.refTime,
      [forecastHeader.forecastType]: {
        [forecastHeader.forecastTime]: dataValue,
      },
    });
    spot.forecasts.push(forecastData);

    await forecastData.save();
  } else {
    // if forecast exists update data
    forecastFound[forecastHeader.forecastType] = {
      ...forecastFound[forecastHeader.forecastType],
      [forecastHeader.forecastTime]: dataValue,
    };
    await forecastFound.save();
  }
  await spot.populate({
    path: 'forecasts',
    match: { forecastInfo: forecastInfo._id },
  });
};

const populateSpots = async (filename, spots) => {
  const forecastJson = await getJson(filename, {
    scriptPath: './src/convert_grib/grib2json/src/bin/grib2json',
    names: true, // (default false): Return descriptive names too
    data: true, // (default false): Return data, not just headers
  });

  // get grib info from header
  const forecastHeader = getforecastHeader(forecastJson[0].header, filename);
  // get forecastInfo document from db or create new one
  const forecastInfo = await getForecastInfo(forecastHeader);

  // Calculate data values for all spots in parallel
  const dataValuePromises = spots.map((spot) => {
    // check if spot is in model borders
    if (inGrid(spot, forecastHeader)) {
      return calculateDataValue(spot, forecastHeader, forecastJson[0].data);
    }
    return null;
  });

  // Wait for all dataValue promises to resolve
  const dataValues = await Promise.all(dataValuePromises);

  // Update spot forecasts with calculated data values
  spots.forEach((spot, index) => {
    if (dataValues[index] !== null) {
      updateSpotForecast(spot, forecastInfo, forecastHeader, dataValues[index]);
    }
  });

  await forecastInfo.save();
};

const convertGrib = async (filenames, path) => {
  const spots = await Spot.find({}).populate('forecasts').exec();
  if (!spots) {
    return false;
  }
  try {
    for (const filename of filenames) {
      await populateSpots(`${path}/${filename}`, spots);
    }
    for (const spot of spots) {
      await spot.save();
    }
  } catch (err) {
    console.log(err);
    return false;
  }
  return true;
};

module.exports = {
  convertGrib,
};
10
Subscribe to my newsletter

Read articles from Tobias Steinhagen directly inside your inbox. Subscribe to the newsletter, and don't miss out.

JavaScriptNode.jsMongoDBMERN Stackforecast

Written by

Tobias Steinhagen
Tobias Steinhagen

My name is Tobias Steinhagen. I'm a self tought full-stack web developer based in Luebeck, Germany. I am also an educated mechanical engineer. I am pursue a career as a web developer