euroargodev / terrain-following-example

0 stars 0 forks source link

Retrieve bathymetry with argopy #2

Open gmaze opened 2 years ago

gmaze commented 2 years ago

There is no need for users to retrieve the multi-Gb GEBCO topography files, argopy has a utility function to retrieve it:

from argopy.utilities import TopoFetcher

The topo fetcher require a box with [lon_min, lon_max, lat_min, lat_max]. One could use the Pandas dataframe with the float trajectory coordinates to retrieve topo over a slightly larger than needed area:

box = [df.where(df['lons']<999)['lons'].min(), 
       df.where(df['lons']<999)['lons'].max(), 
       df.where(df['lats']<999)['lats'].min(), 
       df.where(df['lats']<999)['lats'].max()]
box = np.array(box) + np.array([-5, 5, -5, 5])  # Enlarge the box
print(box)

We can then retrieve topography:

# Retrieve full resolution topo (1/240 deg):
topo = TopoFetcher(box, ds='gebco', cache=True).to_xarray()  
# Or retrieve lower resolution (1/60 deg):
topo = TopoFetcher(box, ds='gebco', stride=[4, 4], cache=True).to_xarray()

The topo xarray dataset could then be used in the example notebook to get topographic estimates at specific lat/lon points

imab4bsh commented 2 years ago

To make this work I need to first understand what the function getDepth does (defined in the script https://github.com/euroargodev/terrain-following/blob/main/interpolation_v2.py). Being lazy and trying to retrieve the entire grid, by defining the box as box = [-180, 180, -90, 90] gives a time out error :-P

gmaze commented 2 years ago

"Being lazy and trying to retrieve the entire grid" won't be a good solution in here 😨

gmaze commented 2 years ago

and since the topo array will not be the global array but a sub-domain around the float, it's the entire getDepth function that should be modified

gmaze commented 2 years ago

but using xarray, this could go easily along something like:

def getDepth(lon,lat,topo,message=False): # this function returns distance-weighted average of depths
    return topo.interp(latitude=lat, longitude=lon, method='cubic')['elevation'].data

with topo from argopy TopoFetcher

gmaze commented 2 years ago

@imab4bsh when you'll have a working notebook of the procedure based on argopy and not relying on some files, it would be great if you give some thoughts into adding this as a new feature in argopy ! a great "new first issues" ! (and obviously I will help you along the way if necessary)