# Extracting data from ACCESS-OM2-01 using coordinate pairs

**URL:** <https://forum.access-hive.org.au/t/extracting-data-from-access-om2-01-using-coordinate-pairs/830>\
**Category:** Technical\
**Created:** [2 June 2023 07:08 UTC](https://forum.access-hive.org.au/t/extracting-data-from-access-om2-01-using-coordinate-pairs/830 "2023-06-02T07:08:08Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![lidefi87](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/lidefi87/32/91_2.png) [@lidefi87](https://forum.access-hive.org.au/u/lidefi87)\
**Post date:** [2 June 2023 07:08 UTC](https://forum.access-hive.org.au/t/extracting-data-from-access-om2-01-using-coordinate-pairs/830/1 "2023-06-02T07:08:08Z")

</div>

I am planning to model the habitat range of a seal in the Southern Ocean, but to do this I need to extract data for some model outputs for each location where a seal has been recorded.

Initially, I decided to assign a reference system ([EPSG 4326](https://epsg.io/4326)) to model outputs and save the results as a georeference image (tif file). The problem is that given the model grid is not even (size/area of grid cells change with latitude), I get some large differences between model outputs and observations.

The image below shows this issue. The data is the bathymetry (from GEBCO, but already regridded to match ACCESS-OM2-01). The areas in grey are the land masks in the bathymetry data. The dots are all the observations for which I need to extract data, and the purple polygons are the location of continents (i.e., areas where land mask should overlap).

 ![image](https://us1.discourse-cdn.com/flex020/uploads/access1/original/1X/bc8c3b44bbb48f38338fd34473cec464c51dd8f3.jpeg)

I have a couple of questions:

1. Has anyone been able to extract data for grid cell closest to a coordinate pairs (i.e., seal location in my case) from the netcdf files directly?, **OR**
2. Does anyone know how I can save my uneven grid in the netcdf files into a georeferenced raster? Once I have this, I can extract the data I need really easily.

Any ideas are appreciated. I feel a little stuck with this.

Thanks!

Denisse

---

<div class="post-metadata">

**Author:** ![dougiesquire](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/dougiesquire/32/16_2.png) [@dougiesquire](https://forum.access-hive.org.au/u/dougiesquire)\
**Post date:** [2 June 2023 09:16 UTC](https://forum.access-hive.org.au/t/extracting-data-from-access-om2-01-using-coordinate-pairs/830/2 "2023-06-02T09:16:27Z")

</div>

Hi @lidefi87. I possibly haven’t properly grokked what you’re wanting to do here, but [`xoak`](https://xoak.readthedocs.io/en/latest/) might be helpful here. E.g. see the [example in the introduction](https://xoak.readthedocs.io/en/latest/examples/introduction.html).

Note that the sort of functionality that `xoak` provides is/will be available natively in `xarray` with the ongoing [flexible indexes refactor](https://github.com/pydata/xarray/blob/main/design_notes/flexible_indexes_notes.md)

---

<div class="post-metadata">

**Author:** ![dougiesquire](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/dougiesquire/32/16_2.png) [@dougiesquire](https://forum.access-hive.org.au/u/dougiesquire)\
**Post date:** [2 June 2023 09:32 UTC](https://forum.access-hive.org.au/t/extracting-data-from-access-om2-01-using-coordinate-pairs/830/3 "2023-06-02T09:32:34Z")

</div>

Alternatively, you could interpolate the data at the seal locations using something like [xesmf](https://xesmf.readthedocs.io/en/latest/)?

---

<div class="post-metadata">

**Author:** ![lidefi87](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/lidefi87/32/91_2.png) [@lidefi87](https://forum.access-hive.org.au/u/lidefi87)\
**Post date:** [5 June 2023 00:57 UTC](https://forum.access-hive.org.au/t/extracting-data-from-access-om2-01-using-coordinate-pairs/830/4 "2023-06-05T00:57:05Z")

</div>

Thanks @dougiesquire. I have not managed to test this option because `xoak` is not currently available in the conda environments in gadi, and I am having trouble installing it. I think I will have to stick to the using `da.sel()` and providing the coordinates for each point.

---

<div class="post-metadata">

**Author:** ![CloLanglais](https://avatars.discourse-cdn.com/v4/letter/c/f04885/32.png) [@CloLanglais](https://forum.access-hive.org.au/u/CloLanglais)\
**Post date:** [5 June 2023 01:46 UTC](https://forum.access-hive.org.au/t/extracting-data-from-access-om2-01-using-coordinate-pairs/830/5 "2023-06-05T01:46:20Z")

</div>

Hi Denisse,  
for the interpolation solution, you can also use LinearNDInterpolator which performs a Delaunay triangulation interpolation using Qhull. The inputs are given as vectors (Lon/lat pairs + data values at that location), so the origin grid does not matter anymore (it does not have to be regular, and can be even unstructured). You can also get rid of the land points and only keep the ocean points, which solve any mask issue. the outputs are a list of lon/lat pairs as well, so works well for observations locations.  
from scipy.interpolate import LinearNDInterpolator

Note that because it is based on triangulation, extrapolation outside of the convex hull is not possible.  
If you have obs which are on the edge of a land point and end up outside of the convex hull of the ocean point, this will not work well.

When I need to extrapolate, I use xesmf regridder. I use the same orgin and destination grid, to “fill” the land point with ocean values.

cheers  
Clo

---

<div class="post-metadata">

**Author:** ![Aidan](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/aidan/32/42_2.png) [@Aidan](https://forum.access-hive.org.au/u/Aidan)\
**Post date:** [5 June 2023 09:23 UTC](https://forum.access-hive.org.au/t/extracting-data-from-access-om2-01-using-coordinate-pairs/830/6 "2023-06-05T09:23:32Z")

</div>

> [@lidefi87](#):
>
> I think I will have to stick to the using `da.sel()` and providing the coordinates for each point.

I’m pretty sure you can do lookups for entire vectors of positions with [`xarray` advanced interpolation](https://docs.xarray.dev/en/stable/user-guide/interpolation.html#advanced-interpolation).

See this [stackoverflow answer](https://stackoverflow.com/a/53952202/4727812) for an example of advanced interpolation.

I thought this is what is used for the langrangian particle tracking work, but @hrsdawson would know better.

---

<div class="post-metadata">

**Author:** ![dougiesquire](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/dougiesquire/32/16_2.png) [@dougiesquire](https://forum.access-hive.org.au/u/dougiesquire)\
**Post date:** [5 June 2023 10:03 UTC](https://forum.access-hive.org.au/t/extracting-data-from-access-om2-01-using-coordinate-pairs/830/7 "2023-06-05T10:03:22Z")

</div>

Good point @Aidan - no need to worry about the curvilinear grid down there

---

<div class="post-metadata">

**Author:** ![Scott](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/scott/32/37_2.png) [@Scott](https://forum.access-hive.org.au/u/Scott)\
**Post date:** [5 June 2023 10:20 UTC](https://forum.access-hive.org.au/t/extracting-data-from-access-om2-01-using-coordinate-pairs/830/8 "2023-06-05T10:20:23Z")

</div>

You should be able to interpolate to points with xesmf using the methods that it supports

- [regrid\_to\_points.ipynb · GitHub](https://gist.github.com/ScottWales/7f5eb086cf28ccabcfdbe6ced31a2b57)

---

<div class="post-metadata">

**Author:** ![dougrichardson](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/dougrichardson/32/523_2.png) [@dougrichardson](https://forum.access-hive.org.au/u/dougrichardson)\
**Post date:** [6 July 2023 05:07 UTC](https://forum.access-hive.org.au/t/extracting-data-from-access-om2-01-using-coordinate-pairs/830/9 "2023-07-06T05:07:15Z")

</div>

Hi @Scott, will your notebook work with a little adapting to go the other way i.e. from a `.tif` to an `xarray` object in lat/lon coords?

---

<div class="post-metadata">

**Author:** ![Scott](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/scott/32/37_2.png) [@Scott](https://forum.access-hive.org.au/u/Scott)\
**Post date:** [6 July 2023 05:32 UTC](https://forum.access-hive.org.au/t/extracting-data-from-access-om2-01-using-coordinate-pairs/830/10 "2023-07-06T05:32:31Z")

</div>

If .tif is a list of lat, lon, value data then yeah I guess you could do it in a similar way, though I’ve not tried it. You may want to turn on [extrapolation](https://xesmf.readthedocs.io/en/latest/notebooks/Masking.html#Extrapolation) to fill in the space between the input data points
