# Tag Info

8

Xarray is (intentionally) ignorant of coordinate systems, so it has no special handling for cyclic coordinates such as longitude. Because your longitude array has only increasing values, xarray interprets selections like slice(40, -80) in the same way that x[i:j] works if x is a NumPy array and i > j >= 0, and thus returns an empty selection. The ...

4

The problem turned out to be a bug in GDAL, which is now fixed: https://trac.osgeo.org/gdal/ticket/6870.

4

This documentation from xarray outlines quite simply the solution to the problem. xarray allows you to interpolate in multiple dimensions and specify another Dataset's x and y dimensions as the output dimensions. So in this case it is done with # interpolation based on http://xarray.pydata.org/en/stable/interpolation.html # interpolation can't be done ...

4

Upon inspecting the dataset, I realized that the units of the data are in radians. import xarray import rioxarray from pyproj import CRS xds = xarray.open_dataset("OR_ABI-L2-LSTF-M6_G17_s20200341900321_e20200341909388_c20200341910038.nc") Inside the x variable the attributes say the data is in radians: xds.x.attrs {'units': 'rad', 'axis': 'X', ...

4

You can use xarray.where. See: http://xarray.pydata.org/en/stable/indexing.html xds = xarray.open_dataset("TT_199501_daymean.nc") xds.where((xds.lon==5.8252) & (xds.lat==46.9359), drop=True) UPDATE: You can find the nearest neighbor using a KDTree. See: https://stackoverflow.com/questions/10818546/finding-index-of-nearest-point-in-numpy-arrays-...

4

3

xarray.Dataset has a to_netcdf method...so it should be clipped.to_netcdf(path_on_disk) http://xarray.pydata.org/en/stable/generated/xarray.Dataset.to_netcdf.html

3

To open the data with the projection information you need to open the sub-datasets individually. I will use a MODIS dataset I have to hand as an example, MOD11A1, but it will be the same for yours. You can get the filename of the subdatasets using rasterio for example: import rasterio filename = '/data/MOD11A1.A2019225.h17v03.006.2019226085002.hdf' with ...

2

I believe what you are looking for is rioxarray. An example of what you want to do is at: https://corteva.github.io/rioxarray/stable/examples/clip_geom.html import rioxarray geometries = [ { 'type': 'Polygon', 'coordinates': [[ [425499.18381405267, 4615331.540546387], [425499.18381405267, 4615478.540546387], ...

2

It's impossible to do this kind of reshape [(100, 256, 256) -> (100, 256, 256,3)]. It's only possible a compatible reshape. This is an example: >>>import numpy as np >>>list = range(100*256*256) >>>array = np.reshape(list, (100, 256, 256)) >>>array = np.reshape(list, (100, 256, 256, 3)) #your reshape: I got an error! ...

2

You could use rioxarray: https://corteva.github.io/rioxarray/stable/examples/clip_box.html import rioxarray min_lon = -24.995 min_lat = 25.05 max_lon = 45.50 max_lat = 71.55 subset = band.rio.clip_box(minx=min_lon, miny=min_lat, maxx=max_lon, maxy=max_lat)

2

One way is to create a boolean mask for the dataset coordinates using the extent you specified and then using the .where() method on the dataset. Here is one example using a tutorial dataset that comes with xarray. First, load the dataset (passing the decode_times=False argument because, at least in my case, it raises an error otherwise) and inspect it. ...

2

I would use xarray.Dataset.resample. xarray is very handy for time series and comes with several temporal dimension handling methods, some of which taken from the comprehensive pandas. Check the docs here. If this does not solve the issue, please add some more context as Rowan asked in the comment. import numpy as np import pandas as pd import xarray as xr ...

2

You can try using Digital Earth Australia's xr_rasterize function to convert your geopandas geodataframe into an xarray object, and then use xarray's .where() method to mask you're array. e.g. mask = xr_rasterize(gdf, da) masked_da = da.where(mask) If you would prefer to use rasterio.features.geometry_mask, then the following code should work. You may ...

2

Since a Dataset can contain more than one variable, it would not be trivial for software to know what you want it to display. If you select one of its variables as DataArray however, that can easily be displayed (if it is in an appropriate format). The variables of a Dataset are available via its data_vars attribute. This exposes the variables as DataArrays ...

2

The error hints at xarray trying to find a variable called "zlib" in your data. The correct structure for the encoding dict would be something like: encode = {"precipitation": {'zlib': True, ...}} But due to the way the data was loaded this is tricky. What you have after loading a file using open_rasterio is a DataArray. A DataArray does ...

1

Using this code: import geopandas import rioxarray import xarray from shapely.geometry import mapping ERA5_daily = xarray.open_dataarray('ERA5_VIMD_2018_08.nc') ERA5_daily.rio.write_crs("EPSG:4326", inplace=True) shapefile = geopandas.read_file('Cluster_34586.shp', crs="EPSG:4326") clipped = ERA5_daily.rio.clip(shapefile.geometry.apply(...

1

If you add the backend_kwargs line, it will work properly, for example:- ds = xr.open_dataset('gfs_4_20191010_0000_006.grb2', engine='cfgrib', backend_kwargs={'filter_by_keys': {'typeOfLevel': 'surface'}}) You can specify type of level you want to do the analysis specifically, the developers says cfgrib and its ...

1

I found a reasonably fast solution combining the answers in https://stackoverflow.com/questions/44681828/efficient-computation-of-minimum-of-haversine-distances and Finding closest point to shapefile coastline Python The code that works now looks like this: import geopandas as gpd from shapely.geometry import Point, box from random import uniform from ...

1

You could use the haversine package, its quite easy to use. From their documentation: from haversine import haversine, Unit lyon = (45.7597, 4.8422) # (lat, lon) paris = (48.8567, 2.3508) haversine(lyon, paris) # in kilometers so for what you want you would need: haversine(lyon, paris, unit=Unit.METERS) # in meters

1

The problem with your dataset is that the lat/lon data is masked: <xarray.DataArray 'lon' (north_south: 320, east_west: 300)> array([[ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, ...

1

It ended up being a bug that was fixed in version 0.0.27

1

How to mask NetCDF time series data from a shapefile in Python? You can use rioxarray. Here is an example: https://corteva.github.io/rioxarray/stable/examples/clip_geom.html import rioxarray import geopandas geodf = geopandas.read_file(...) xds = rioxarray.open_rasterio(...) clipped = xds.rio.clip(geodf.geometry.apply(mapping), geodf.crs)

1

I think this should fix it: ds = ds.assign_coords(x=(((ds.x + 180) % 360) - 180)).sortby('x') The reason is that rioxarray uses the 1-d coordinate and lon is 2-d. <xarray.Dataset> Dimensions: (bnds: 2, time: 2407, x: 360, y: 180) Coordinates: ... * x (x) float64 0.5 1.5 2.5 3.5 4.5 ... 356.5 357.5 358.5 359.5 ... lon (...

1

Since there is none answer since question was posted, it might be helpful to tell how I dealt with this issue. First of all, I'm still not familiar with .netcdf format and therefore it took me some time to suspect that the problem wasn't exclusively this lack of familiarity nor some intrinsic complexity of the format, but the apparent absence of an explicit ...

1

I think one of the most evident advantages is indexing. Consider the following example where you have data and both longitude and latitude stored in 2D arrays: data = np.random.randint(100, 1000, size=(4, 4)) data [[176, 479, 713, 973], [992, 259, 969, 355], [182, 139, 633, 938], [761, 911, 124, 855]] x = np.linspace(-76, -74.5, 4) y = np.linspace(-5.0,...

1

I ended up writing a function to parse the structured metadata attribute and add coordinates to a dataset. Here is my code to read in the dataset, create the coordinates, and apply it to the dataset: import xarray as xr import pyproj import numpy as np from collections import OrderedDict def parse_hdfeos_metadata(string): out = OrderedDict() lines = [i....

Only top voted, non community-wiki answers of a minimum length are eligible