Select data by latitude and longitude

Question

I am using a dataset from DWD (Deutscher Wetterdienst) and want to select data by latitude and longitude. The import works so far. So no problem there. Now I want to select data by latitude and longitude. It works when I try to select data with sel when I use x and y.

But not with lat and long. I tried all the answer which I could find, like:

ds.sel(latitude=50, longitude=14, method='nearest')

but I am getting the error

ValueError: dimensions or multi-index levels ['latitude', 'longitude'] do not exist

That's my code:

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import xarray as xr


ds = xr.open_dataset(
            'cosmo-d2_germany_rotated-lat-lon_single-level_2019061721_012_ASWDIFD_S.grib2',
            engine='cfgrib',
            backend_kwargs={'filter_by_keys': {'stepUnits': 1}}
        )

print(ds)

Output:

<xarray.Dataset>
Dimensions:     (x: 651, y: 716)
Coordinates:
    time        datetime64[ns] ...
    step        timedelta64[ns] ...
    surface     int32 ...
    latitude    (y, x) float64 ...
    longitude   (y, x) float64 ...
    valid_time  datetime64[ns] ...
Dimensions without coordinates: x, y
Data variables:
    ASWDIFD_S   (y, x) float32 ...
Attributes:
    GRIB_edition:            2
    GRIB_centre:             edzw
    GRIB_centreDescription:  Offenbach 
    GRIB_subCentre:          255
    Conventions:             CF-1.7
    institution:             Offenbach 
    history:                 2019-07-22T13:35:33 GRIB to CDM+CF via cfgrib-

I'm unable to repeat this problem with well formatted grib or grib2 files I have on hand. Can you post a link to the grib2 file you are working with? — , Jul 22 '19 at 19:22
You can download recent files under https://opendata.dwd.de/weather/nwp/cosmo-d2/grib/21/aswdifd_s/ — laba.rababa, Jul 23 '19 at 07:53

score 3 · Answer 1 · edited Oct 03 '19 at 22:44

3

In your file latitude and longitude are not dimensions but rather helper 2D variables containing coordinate data. In xarray parlance they are called non-dimension coordinates and you cannot slice on them. See also Working with Multidimensional Coordinates.

edited Oct 03 '19 at 22:44

Asclepius

57,944
17
167
143

answered Jul 23 '19 at 06:42

alexamici

754
5
10

Thanks for reference. But i cant find any hints there how solve to problem. When I select one dataset by x, y it returns me the corresponding lat and long. Do have to convert the non-dimension coordnates? I think the problem is that i have 2-dimensional lats and longs. – laba.rababa Jul 23 '19 at 07:13

score 1 · Answer 2 · answered Aug 08 '19 at 11:41

It would be better if you regrid the data to a regular grid inside python so that you have latitudes and longitudes as 1D vectors, you would have to make a grid and then interpolate the data over that grid.

Also you need to check https://www.ecmwf.int/sites/default/files/elibrary/2018/18727-cfgrib-easy-and-efficient-grib-file-access-xarray.pdf to see the way to access grib files in xarray. If you dont want to use xarray for this purpose pygrib is another option.

score 1 · Answer 3 · answered Nov 08 '19 at 09:10

I can't test the solution as I don't have the cfgrib engine installed, but could you try to use

numpy.find_nearest(lonarray, lonvalue)

to find the lon and lat indexes near your point as per this soln: Find nearest value in numpy array

And then select the point using the index directly on the x,y coordinates?

http://xarray.pydata.org/en/stable/indexing.html

and-bri · Answer 4 · 2019-11-24T11:18:37.247

i wrote a function for the files from the DWD:

import pygrib   # https://jswhit.github.io/pygrib/docs/
import numpy as np

def get_grib_data_nearest_point(grib_file, inp_lat, inp_lon):
    """
    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
    :return:            scalar
    """

    # open the grib file, get the coordinates and values
    grbs = pygrib.open(grib_file)
    grb = grbs[1]
    lats, lons = grb.latlons()
    values = grb.values
    grbs.close()

    # check if user coords are valide
    if inp_lat > max(grb.distinctLatitudes): return np.nan
    if inp_lat < min(grb.distinctLatitudes): return np.nan
    if inp_lon > max(grb.distinctLongitudes): return np.nan
    if inp_lon < min(grb.distinctLongitudes): return np.nan

    # find index for  closest lat (x)
    diff_save = 999
    for x in range(0, len(lats)):
        diff = abs(lats[x][0] - inp_lat)
        if diff < diff_save:
            diff_save = diff
        else:
            break

    # find index for closest lon (y)
    diff_save = 999
    for y in range(0, len(lons[x])):
        diff = abs(lons[x][y] - inp_lon)
        if diff < diff_save:
            diff_save = diff
        else:
            break

    # index the array to return the correspondent value
    return values[x][y]

I believe this code is incorrect. The answer is off by 1. the return should be `return values[ x - 1 ][ y - 1 ]` — kdbaseball8, Mar 28 '22 at 21:00

score 0 · Answer 5 · answered Dec 16 '19 at 15:18

As noted above, you can re-grid your data (probably given in curvilinear grid i.e., lat and lon in 2D arrays) to your desired resolution of 1-D array (lat/lon) , after which you can use .sel directly on the lat/lon coords to slice the data.

Check out xESMF(https://xesmf.readthedocs.io/en/latest/notebooks/Curvilinear_grid.html).

Easy, fast interpolation and regridding of Xarray fields with good examples and documentation.

Select data by latitude and longitude

5 Answers5