<p>First of all: <strong>NEVER, EVER, EVER</strong> (I repeat, <strong>EVER</strong>) loop in <code>Python</code> over an array to subset some values :)</p>
<p><a href="https://docs.xarray.dev/en/latest/index.html" rel="nofollow noreferrer">xarray</a> was created exactly to avoid having to write something like this by exploiting the native coordinates/dimensions indexing. For a brief introduction on how to subset an array</p>
<p>List item</p>
<p>have a look <a href="https://docs.xarray.dev/en/latest/user-guide/indexing.html" rel="nofollow noreferrer">here</a></p>
<p>Regarding your specific issue, it is not clear what you are trying to achieve because you don't state what you want to &quot;extract&quot; in the description. Judging from the code, it seems you're trying to downsample the grid to a specific resolution (<code>res</code>). In pure <code>xarray</code> going from e.g. 0.25 to 0.5 degrees resolution would be as simple as doing</p>
<pre><code>import xarray as xr

ds = xr.open_dataset(&quot;test.nc&quot;)

ds.sel(lon=ds.lon[::2], lat=ds.lat[::2])
</code></pre>
<p>There are many other ways of subsetting data to extract a single point or set of points using the <code>sel</code> operator, for example</p>
<pre><code>ds.sel(lon=43, lat=0, method='nearest') # single point
</code></pre>
<p>There are many functions that are written to do operations on geospatial multi-dimensional arrays (interpolation, subsetting,...). As a general rule of thumb, if you're writing your own low-level function to access element one by one in Python...you're doing it wrong :D There's almost surely another function that does that in a vectorized way!</p>
<p>Please take some time to learn how to use <code>xarray</code> to work with Netcdf files and subset data, you'll save a lot of time in the future!</p>
<p>P.S. Even when using <code>xarray</code>, <code>Python</code> remains quite slow with large dataset. If you need to do operations on large files I'd suggest you to pre-process the file BEFORE loading it into <code>Python</code> using something like <a href="https://code.mpimet.mpg.de/projects/cdo" rel="nofollow noreferrer">CDO</a></p>

CC BY-SA 4.0
2024-02-05T14:13:36.887Z
As mentioned above, I want to extract ocean current zonal and meridional vectors which are variables vozocrte and vomecrtn in this dataset. Your selection statement gives me a subset of lons/lats. This will not yield a regular grid. I understand this is not clear from seeing the code alone. That's why I mentioned the dataset and type of dataset itself.
mRcSchwering
CC BY-SA 4.0
2024-02-05T17:14:12.763Z
So you basically want to subset in lat/lon, right? Why not use `xr.where`? This will return a regular grid in both cases, unless you specify `drop=True`. You just have to write your condition where you want to keep the data (it is still not clear to me what this condition is :)). But honestly, for your own sake, try to avoid loops as much as possible. Also, try to make a self-contained and self-explanatory MWE next time, it would be easier for people to respond.
Droid
CC BY-SA 4.0
2024-02-05T13:26:04.680Z
2024-02-05T13:26:04.680Z
Droid
<p>There is a very simple way to do this.
The point here is that these arrays are all ordered consistently.
The dataset needs only to be subset by depth if all <code>lon,lat</code> pairs should be extracted.
Lat and lon values can still be corrected later on to create a specific resolution.</p>
<p>What I ended up doing is in essence this:</p>
<pre><code>import numpy as np
import pandas as pd
from netCDF4 import Dataset


def _rm_mask(masked: np.ma.MaskedArray, fill=np.nan) -&gt; np.ndarray:
    return np.where(masked.mask, fill, masked.data).astype(masked.dtype)


def _digitize(df: pd.DataFrame, var: str, interval: tuple[float, float], res: float):
    for target in np.arange(interval[0], interval[1] + res, res):
        mask = (df[var] &gt;= target - res / 2) &amp; (df[var] &lt; target + res / 2)
        df.loc[mask, var] = target


ds = Dataset(file, &quot;r&quot;, format=&quot;NETCDF4&quot;)
depths = _rm_mask(ds.variables[&quot;depth&quot;][:])  # (75,) depth in m
lats = _rm_mask(ds.variables[&quot;nav_lat&quot;][:])  # (1021, 1442) lat in degrees
lons = _rm_mask(ds.variables[&quot;nav_lon&quot;][:])  # (1021, 1442) lon in degrees
vals = _rm_mask(ds.variables[&quot;vozocrte&quot;][:])  # (76, 1021, 1442) in m/s
vals = vals[0, depths &lt; 5]

dfdict = {&quot;lat&quot;: lats.flatten().tolist(), &quot;lon&quot;: lons.flatten().tolist()}
for layer_i in range(vals.shape[0]):
    dfdict[f&quot;l{layer_i}&quot;] = vals[layer_i].flatten().tolist()

df = pd.DataFrame(dfdict)
_digitize(df=df, var=&quot;lon&quot;, interval=(-180, 180), res=0.25)
_digitize(df=df, var=&quot;lat&quot;, interval=(-70, 70), res=0.25)
df = df.groupby([&quot;lon&quot;, &quot;lat&quot;]).mean()
df.sort_index(inplace=True)
</code></pre>
<p>Here, most heavy lifting is done by numpy and is multithreaded.
<code>_digitize</code> can also be written using <a href="https://numpy.org/doc/stable/reference/generated/numpy.digitize.html" rel="nofollow noreferrer">np.digitize</a> but it looks a bit confusing (mapping to bin indices) and doesn't make such a huge difference.
Overall this computes ~36 times faster (now just a few seconds).</p>
<p>However, as Droid mentioned <a href="https://docs.xarray.dev/en/latest/index.html" rel="nofollow noreferrer">xarray</a> might be even faster.</p>

CC BY-SA 4.0
2024-02-05T14:38:38.480Z
2024-02-05T14:38:38.480Z
mRcSchwering
<p>I am working on the <a href="https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-oras5?tab=overview" rel="nofollow noreferrer">ORAS5</a> dataset which comes as an archive of netCDF files. Currently, I am using python and <a href="https://unidata.github.io/netcdf4-python" rel="nofollow noreferrer">netcdf4</a> to process it.</p>
<p>My goal is to get surface level ocean currents across the globe (almost, 70°N to 70°S).
The dataset itself has several variables.
These are relevant for my problem:</p>
<ul>
<li><strong>depth</strong> 1D (75,) array describing increasing depths</li>
<li><strong>lat</strong> 2D (1021, 1442) array describing latitude</li>
<li><strong>lon</strong> 2D (1021, 1442) array describing longitude</li>
<li><strong>vozocrte/vomecrtn</strong> 3D (75, 1021, 1442) array describing zonal or meridional velocity</li>
</ul>
<p>Currently, I am creating a lat-lon grid with desired resolution and then loop over this grid to extract <code>vozocrte</code> or <code>vomecrtn</code>.
In principle it looks like this:</p>
<pre><code>import numpy as np
from netCDF4 import Dataset

def mask_map(vals: np.ndarray, range: tuple[int, int], res=0.25) -&gt; dict[float, np.ndarray]:
    return {
        d: (vals &gt;= d - res / 2) &amp; (vals &lt;= d + res / 2)
        for d in np.arange(range[0], range[1] + res, res)
    }

ds = Dataset(file, &quot;r&quot;, format=&quot;NETCDF4&quot;)
depths = ds.variables[&quot;depth&quot;][:]  # (75,) depth in m
lats = ds.variables[&quot;nav_lat&quot;][:]  # (1021, 1442) lat in degrees
lons = ds.variables[&quot;nav_lon&quot;][:]  # (1021, 1442) lon in degrees
vals = ds.variables[&quot;vozocrte&quot;][:]  # (76, 1021, 1442) in m/s
vals = vals[0, depths &lt; 5]
lon_2_mask = mask_map(vals=lons, range=lon_range)
lat_2_mask = mask_map(vals=lats, range=lat_range)

records = []
for lon, lons_mask in lon_2_mask.items():
    for lat, lats_mask in lat_2_mask.items():
        masked = vals[..., lats_mask &amp; lons_mask]
        if not masked.mask.all():
            records.append({&quot;component&quot;: masked.mean(), &quot;lat&quot;: lat, &quot;lon&quot;: lon})
</code></pre>
<p>It works but it is very slow.</p>
<p>I am unfamiliar with this type of dataset (mapping depth-lon-lat to z-y-x, then using z-y-x to get some values).
Is there a more efficient way of doing this?</p>
<p><strong>Edit</strong></p>
<p>I realized I should explain the data better.
The variables <code>depth</code>, <code>nav_lat</code>, <code>nav_lon</code>, <code>vozocrte</code>/<code>vomecrtn</code> mentioned above all contain masked numpy arrays of different shape.</p>
<ul>
<li><code>vozocrte</code>/<code>vomecrtn</code> describe a velocity in m/s as <code>velocity(z, y, x)</code>.</li>
<li><code>depth</code> describes depth in m as <code>depth(z)</code>.</li>
<li><code>nav_lon</code> describes longitude in degree as <code>longitude(y, x)</code>.</li>
<li><code>nav_lat</code> describes latitude in degree as <code>latitude(y, x)</code>.</li>
</ul>
<p><em>E.g.</em> <code>nav_lat</code> can contain a masked array of <code>[[87.021, 86.721, 86.101, ...], ...]</code>.
And <code>nav_lon</code> can contain a masked array of <code>[[178.131, 177.809, 176.980, ...], ...]</code>.
<code>depth</code> a masked array of <em>e.g.</em> <code>[0.023, 1.201, ...]</code>.</p>
<p>My naive approach above follows the logic of first getting a mask for <code>z, y, x</code> from the variables depth, nav_lat, nav_lon, then using them to get velocities.
<em>E.g.</em> if I wanted 10°-11°N, 15°-16°E, less than 5m deep:</p>
<pre><code>depth_mask = depths &lt; 5
lon_mask = (lons &gt;= 15) &amp; (lons &lt; 16)
lat_mask = (lats &gt;= 10) &amp; (lats &lt; 11)
vals[depth_mask, ..., lon_mask &amp; lat_mask]
</code></pre>

<ocean-currents><python><netcdf>
Efficient way to parse lat-lon grid at certain depth in tripolar model from nc file?
mRcSchwering
# Efficient way to parse lat-lon grid at certain depth in tripolar model from nc file?

(Asked by: username_0 on 2024-02-03)

I am working on the [ORAS5](https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-oras5?tab=overview) dataset which comes as an archive of netCDF files. Currently, I am using python and [netcdf4](https://unidata.github.io/netcdf4-python) to process it.

My goal is to get surface level ocean currents across the globe (almost, 70°N to 70°S).
The dataset itself has several variables.
These are relevant for my problem:

* **depth** 1D (75,) array describing increasing depths
* **lat** 2D (1021, 1442) array describing latitude
* **lon** 2D (1021, 1442) array describing longitude
* **vozocrte/vomecrtn** 3D (75, 1021, 1442) array describing zonal or meridional velocity

Currently, I am creating a lat-lon grid with desired resolution and then loop over this grid to extract `vozocrte` or `vomecrtn`.
In principle it looks like this:

```
import numpy as np
from netCDF4 import Dataset

def mask_map(vals: np.ndarray, range: tuple[int, int], res=0.25) -> dict[float, np.ndarray]:
    return {
        d: (vals >= d - res / 2) & (vals <= d + res / 2)
        for d in np.arange(range[0], range[1] + res, res)
    }

ds = Dataset(file, "r", format="NETCDF4")
depths = ds.variables["depth"][:]  # (75,) depth in m
lats = ds.variables["nav_lat"][:]  # (1021, 1442) lat in degrees
lons = ds.variables["nav_lon"][:]  # (1021, 1442) lon in degrees
vals = ds.variables["vozocrte"][:]  # (76, 1021, 1442) in m/s
vals = vals[0, depths < 5]
lon_2_mask = mask_map(vals=lons, range=lon_range)
lat_2_mask = mask_map(vals=lats, range=lat_range)

records = []
for lon, lons_mask in lon_2_mask.items():
    for lat, lats_mask in lat_2_mask.items():
        masked = vals[..., lats_mask & lons_mask]
        if not masked.mask.all():
            records.append({"component": masked.mean(), "lat": lat, "lon": lon})

```

It works but it is very slow.

I am unfamiliar with this type of dataset (mapping depth-lon-lat to z-y-x, then using z-y-x to get some values).
Is there a more efficient way of doing this?

**Edit**

I realized I should explain the data better.
The variables `depth`, `nav_lat`, `nav_lon`, `vozocrte`/`vomecrtn` mentioned above all contain masked numpy arrays of different shape.

* `vozocrte`/`vomecrtn` describe a velocity in m/s as `velocity(z, y, x)`.
* `depth` describes depth in m as `depth(z)`.
* `nav_lon` describes longitude in degree as `longitude(y, x)`.
* `nav_lat` describes latitude in degree as `latitude(y, x)`.

*E.g.* `nav_lat` can contain a masked array of `[[87.021, 86.721, 86.101, ...], ...]`.
And `nav_lon` can contain a masked array of `[[178.131, 177.809, 176.980, ...], ...]`.
`depth` a masked array of *e.g.* `[0.023, 1.201, ...]`.

My naive approach above follows the logic of first getting a mask for `z, y, x` from the variables depth, nav\_lat, nav\_lon, then using them to get velocities.
*E.g.* if I wanted 10°-11°N, 15°-16°E, less than 5m deep:

```
depth_mask = depths < 5
lon_mask = (lons >= 15) & (lons < 16)
lat_mask = (lats >= 10) & (lats < 11)
vals[depth_mask, ..., lon_mask & lat_mask]

```

**Tags:** ocean-currents, python, netcdf


---

## 2 answers

### Answer by username_1 on 2024-02-05. Score: 3

First of all: **NEVER, EVER, EVER** (I repeat, **EVER**) loop in `Python` over an array to subset some values :)

[xarray](https://docs.xarray.dev/en/latest/index.html) was created exactly to avoid having to write something like this by exploiting the native coordinates/dimensions indexing. For a brief introduction on how to subset an array

List item

have a look [here](https://docs.xarray.dev/en/latest/user-guide/indexing.html)

Regarding your specific issue, it is not clear what you are trying to achieve because you don't state what you want to "extract" in the description. Judging from the code, it seems you're trying to downsample the grid to a specific resolution (`res`). In pure `xarray` going from e.g. 0.25 to 0.5 degrees resolution would be as simple as doing

```
import xarray as xr

ds = xr.open_dataset("test.nc")

ds.sel(lon=ds.lon[::2], lat=ds.lat[::2])

```

There are many other ways of subsetting data to extract a single point or set of points using the `sel` operator, for example

```
ds.sel(lon=43, lat=0, method='nearest') # single point

```

There are many functions that are written to do operations on geospatial multi-dimensional arrays (interpolation, subsetting,...). As a general rule of thumb, if you're writing your own low-level function to access element one by one in Python...you're doing it wrong :D There's almost surely another function that does that in a vectorized way!

Please take some time to learn how to use `xarray` to work with Netcdf files and subset data, you'll save a lot of time in the future!

P.S. Even when using `xarray`, `Python` remains quite slow with large dataset. If you need to do operations on large files I'd suggest you to pre-process the file BEFORE loading it into `Python` using something like [CDO](https://code.mpimet.mpg.de/projects/cdo)

> **Comments:**
> **username_0**: As mentioned above, I want to extract ocean current zonal and meridional vectors which are variables vozocrte and vomecrtn in this dataset. Your selection statement gives me a subset of lons/lats. This will not yield a regular grid. I understand this is not clear from seeing the code alone. That's why I mentioned the dataset and type of dataset itself.  
> **username_1**: So you basically want to subset in lat/lon, right? Why not use `xr.where`? This will return a regular grid in both cases, unless you specify `drop=True`. You just have to write your condition where you want to keep the data (it is still not clear to me what this condition is :)). But honestly, for your own sake, try to avoid loops as much as possible. Also, try to make a self-contained and self-explanatory MWE next time, it would be easier for people to respond.  

---

### Answer by username_0 on 2024-02-05. Score: 0

There is a very simple way to do this.
The point here is that these arrays are all ordered consistently.
The dataset needs only to be subset by depth if all `lon,lat` pairs should be extracted.
Lat and lon values can still be corrected later on to create a specific resolution.

What I ended up doing is in essence this:

```
import numpy as np
import pandas as pd
from netCDF4 import Dataset

def _rm_mask(masked: np.ma.MaskedArray, fill=np.nan) -> np.ndarray:
    return np.where(masked.mask, fill, masked.data).astype(masked.dtype)

def _digitize(df: pd.DataFrame, var: str, interval: tuple[float, float], res: float):
    for target in np.arange(interval[0], interval[1] + res, res):
        mask = (df[var] >= target - res / 2) & (df[var] < target + res / 2)
        df.loc[mask, var] = target

ds = Dataset(file, "r", format="NETCDF4")
depths = _rm_mask(ds.variables["depth"][:])  # (75,) depth in m
lats = _rm_mask(ds.variables["nav_lat"][:])  # (1021, 1442) lat in degrees
lons = _rm_mask(ds.variables["nav_lon"][:])  # (1021, 1442) lon in degrees
vals = _rm_mask(ds.variables["vozocrte"][:])  # (76, 1021, 1442) in m/s
vals = vals[0, depths < 5]

dfdict = {"lat": lats.flatten().tolist(), "lon": lons.flatten().tolist()}
for layer_i in range(vals.shape[0]):
    dfdict[f"l{layer_i}"] = vals[layer_i].flatten().tolist()

df = pd.DataFrame(dfdict)
_digitize(df=df, var="lon", interval=(-180, 180), res=0.25)
_digitize(df=df, var="lat", interval=(-70, 70), res=0.25)
df = df.groupby(["lon", "lat"]).mean()
df.sort_index(inplace=True)

```

Here, most heavy lifting is done by numpy and is multithreaded.
`_digitize` can also be written using [np.digitize](https://numpy.org/doc/stable/reference/generated/numpy.digitize.html) but it looks a bit confusing (mapping to bin indices) and doesn't make such a huge difference.
Overall this computes ~36 times faster (now just a few seconds).

However, as username_1 mentioned [xarray](https://docs.xarray.dev/en/latest/index.html) might be even faster.


---

