CODEWITHSUNDEEP
pythonpython-xarray
laba.rababa
laba.rababa

Reputation: 31

Select data by latitude and longitude

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-

Upvotes: 3

Views: 7013

Answers (5)

Ziskin_Ziv
Ziskin_Ziv

Reputation: 116

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.

Upvotes: 0

and-bri
and-bri

Reputation: 1664

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]

Upvotes: 0

ClimateUnboxed
ClimateUnboxed

Reputation: 8057

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

Upvotes: 1

alexamici
alexamici

Reputation: 784

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.

Upvotes: 3

Manmeet Singh
Manmeet Singh

Reputation: 400

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.

Upvotes: 1

Related Questions