2

I have a grib file. I would like to be able to do something like this:

grib_data = grib_data_nearest_point(grib_file='grib_file.grib2', lat=15.5, lon=23.3)

Is there a library that can accomplish this in python, or some library that has a workaround that can accomplish this? I am currently reading in every point in the grib file, and performing the distance calculation on every single point, which is incredibly slow.

Here is my comparison of my old method and the new method suggested:

import timeit
import pygrib
import math
import gdal
import numpy as np
from geopy.distance import geodesic

class SpeedGribOld:

    def get_closest_datapoint(self, location_to_match, grib_file):

        grib_file_handler = pygrib.open(grib_file)
        data = grib_file_handler[1].data()[0].flatten()
        grib_lats = grib_file_handler[1].data()[1].flatten()
        grib_lons = grib_file_handler[1].data()[2].flatten()

        match_lat = location_to_match[0]
        match_lon = location_to_match[1]

        distances = []
        grib_indices = []
        for index, grib_lat in enumerate(grib_lats):
            grib_lon = grib_lons[index]

            # Only check grib points within 2 degrees lat/lon of the powerplant
            if abs(grib_lat - match_lat) < 2:
                if abs(grib_lon - match_lon) < 2:
                    distances.append(geodesic((grib_lat, grib_lon), (match_lat, match_lon)))
                    grib_indices.append(index)
        # Get min distance index
        index_min_dist = np.argmin(distances)
        closest_data = data[grib_indices[index_min_dist]]
        print("OLD: Closest data: {}".format(closest_data))

class SpeedGribNew:
    def get_closest_datapoint(self, location_to_match, grib_file, band_i=1):
        """
        Gets the correspondent value to a latitude-longitude pair of coordinates in
        a grib file.

        :param grib_file:   path to the grib file in disk
        :param match_lat:         latitude
        :param match_lon:         longitude
        :param band_i:      band index
        :return:            scalar
        """

        match_lat = location_to_match[0]
        match_lon = location_to_match[1]

        # open the grib file, get the specified band and read it as an array
        ds = gdal.Open(grib_file, 0)
        band = ds.GetRasterBand(band_i)
        arr = band.ReadAsArray()

        # get origin's coordinates, pixel width and pixel height
        # the GetGeoTransform method returns the skew in the x and y axis but you
        # can ignore these values
        ox, pw, xskew, oy, yskew, ph = ds.GetGeoTransform()

        # calculate the indices (row and column)
        i = math.floor((oy - match_lat) / ph)
        j = math.floor((match_lon - ox) / pw)

        # close the file
        del ds, band

        # index the array to return the correspondent value
        closest_data = arr[i, j]
        print("NEW: Closest data: {}".format(closest_data))


grib_files = ['/path/to/GFS_UREL.grib']



start = timeit.default_timer()
SpeedGribOld().get_closest_datapoint(
    location_to_match = np.array([10, 30]),
    grib_file = grib_files[0]
)
stop = timeit.default_timer()
print('Old Time: ', stop - start)


start = timeit.default_timer()
SpeedGribNew().get_closest_datapoint(
    location_to_match = np.array([10, 30]),
    grib_file = grib_files[0]
)
stop = timeit.default_timer()
print('New Time: ', stop - start)

Output:

OLD: Closest data: 2.2722020149230957
Old Time:  2.053238352003973
NEW: Closest data: -3.877797842025757
New Time:  0.03821680700639263
1
  • Can you edit your question to include the code you are using to read the data and perform the distance calculation? – Marcelo Villa-Piñeros Jun 16 '19 at 19:42
1

gdal supports the GRIB format. You can use it to open your file and read a specific band as a numpy array. Then, assuming both your grib file and the x y coordinates share the same spatial reference and that your raster is not skewed, you can calculate the row and column (indices) of the array the coordinates correspond to. To do this you need to know the origin coordinates in x and y and the raster's pixel width and height, which can be obtained from the raster's geotransform. Following this answer, this calculation can be done with the following function:

def get_indices(x, y, ox, oy, pw, ph):
    """
    Gets the row (i) and column (j) indices in an array for a specific point.
    Based on https://gis.stackexchange.com/a/92015/86131

    :param x:   point x coordinate (longitude)
    :param y:   point y coordinate (latitude)
    :param ox:  raster x origin
    :param oy:  raster y origin
    :param pw:  raster pixel width
    :param ph:  raster pixel height
    :return:    row (i) and column (j) indices
    """

    i = math.floor((oy-y) / ph)
    j = math.floor((x-ox) / pw)

    return i, j

Knowing this, you can easily incoporate this functionality into a function like the one you asked for. Here is an example:

import math

import gdal


def grib_data_nearest_point(grib_file, lat, lon, band_i=1):
    """
    Gets the correspondent value to a latitude-longitude pair of coordinates in
    a grib file.

    :param grib_file:   path to the grib file in disk
    :param lat:         latitude
    :param lon:         longitude
    :param band_i:      band index
    :return:            scalar
    """
    # open the grib file, get the specified band and read it as an array
    ds = gdal.Open(grib_file, 0)
    band = ds.GetRasterBand(band_i)
    arr = band.ReadAsArray()

    # get origin's coordinates, pixel width and pixel height
    # the GetGeoTransform method returns the skew in the x and y axis but you
    # can ignore these values
    ox, pw, xskew, oy, yskew, ph = ds.GetGeoTransform()

    # calculate the indices (row and column)
    i = math.floor((oy - lat) / ph)
    j = math.floor((lon - ox) / pw)

    # close the file
    del ds, band

    # index the array to return the correspondent value
    return arr[i, j]
9
  • This is extremely helpful! What is the 'band' of the grib file? – P.V. Jun 16 '19 at 18:26
  • Well, that depends on your data. Rasters can have multiple bands. Think of RGB images: they usually have three bands (one for the red values, one for the green values and one for the blue values). It is possible that your grib file has just one band too. By default, the function I wrote above opens the first band but you can change it by specifying an integer in the band_n parameter. To see how many bands ome of your file has simply call print(ds.RasterCount). – Marcelo Villa-Piñeros Jun 16 '19 at 19:15
  • That makes sense. When I run my code through your function, I get the following error: "Warning: Inside GRIB2Inventory, Message # 2 ERROR: Couldn't find 'GRIB' or 'TDLP'", and arr returns as an array in which every value is 9999.0.. The grib file I am looking at is '000' on this page: awsopendata.s3-website-us-west-2.amazonaws.com/noaa-hrrr. What am I doing wrong? – P.V. Jun 16 '19 at 20:53
  • Can you provide a full link to the download or the grib file? I do not see any 000 file listed on the link you provided. – Marcelo Villa-Piñeros Jun 16 '19 at 21:28
  • Ah sorry, to get to the grib file, I followed this file path on that link: VGRD/80 m HGHT/20190613/1500/000 – P.V. Jun 16 '19 at 23:22
0

To my opinion gives the round function a more accurate result:

i = round((oy - lat) / ph, 0)

j = round((lon - ox) / pw, 0)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.