{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Extract a value from a curvilinear grid\n", "\n", "A small notebook showing how to lookup the ij index from a given longitude and latitude in a curvilinear grid" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import xarray as xr\n", "import numpy as np" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset>\n",
       "Dimensions:               (axis_nbounds: 2, deptht: 31, tbnds: 2, time_counter: 12, x: 182, y: 149)\n",
       "Coordinates:\n",
       "  * time_counter          (time_counter) object 0000-01-16 00:00:00 ... 0000-...\n",
       "    nav_lat               (y, x) float32 -78.19 -78.19 -78.19 ... 59.91 59.91\n",
       "    nav_lon               (y, x) float32 105.0 107.0 109.0 ... 106.0 106.0 106.0\n",
       "  * deptht                (deptht) float32 5.0 15.0 25.0 ... 4.75e+03 5.25e+03\n",
       "    time_centered         (time_counter) object 4679-07-16 00:00:00 ... 4680-...\n",
       "Dimensions without coordinates: axis_nbounds, tbnds, x, y\n",
       "Data variables: (12/64)\n",
       "    time_counter_bnds     (time_counter, tbnds) object 0000-01-01 00:00:00 .....\n",
       "    deptht_bounds         (deptht, axis_nbounds) float32 0.0 10.0 ... 5.501e+03\n",
       "    time_centered_bounds  (time_counter, axis_nbounds) object 4679-07-01 00:0...\n",
       "    time_counter_bounds   (time_counter, axis_nbounds) float64 8.645e+10 ... ...\n",
       "    e3t                   (time_counter, deptht, y, x) float32 ...\n",
       "    pbo                   (time_counter, y, x) float32 ...\n",
       "    ...                    ...\n",
       "    rain                  (time_counter, y, x) float32 ...\n",
       "    snow_ao_cea           (time_counter, y, x) float32 ...\n",
       "    snow_ai_cea           (time_counter, y, x) float32 ...\n",
       "    evap_ao_cea           (time_counter, y, x) float32 ...\n",
       "    subl_ai_cea           (time_counter, y, x) float32 ...\n",
       "    sosflxdo              (time_counter, y, x) float32 ...\n",
       "Attributes:\n",
       "    name:                      C10M-PNG-560_1m_grid_T\n",
       "    description:               Created by xios\n",
       "    title:                     Created by xios\n",
       "    Conventions:               CF-1.6\n",
       "    timeStamp:                 2019-Dec-10 12:48:46 GMT\n",
       "    uuid:                      0b59c577-6754-47de-a26a-a13d8a514539\n",
       "    LongName:                  IPSL 10 Ma,2X, BC from PlioMIP except for Aust...\n",
       "    history:                   Wed Dec 11 05:53:06 2019: ncrcat -C --buffer_s...\n",
       "    NCO:                       "4.6.0"\n",
       "    nco_openmp_thread_number:  1
" ], "text/plain": [ "\n", "Dimensions: (axis_nbounds: 2, deptht: 31, tbnds: 2, time_counter: 12, x: 182, y: 149)\n", "Coordinates:\n", " * time_counter (time_counter) object 0000-01-16 00:00:00 ... 0000-...\n", " nav_lat (y, x) float32 ...\n", " nav_lon (y, x) float32 ...\n", " * deptht (deptht) float32 5.0 15.0 25.0 ... 4.75e+03 5.25e+03\n", " time_centered (time_counter) object ...\n", "Dimensions without coordinates: axis_nbounds, tbnds, x, y\n", "Data variables: (12/64)\n", " time_counter_bnds (time_counter, tbnds) object ...\n", " deptht_bounds (deptht, axis_nbounds) float32 ...\n", " time_centered_bounds (time_counter, axis_nbounds) object ...\n", " time_counter_bounds (time_counter, axis_nbounds) float64 ...\n", " e3t (time_counter, deptht, y, x) float32 ...\n", " pbo (time_counter, y, x) float32 ...\n", " ... ...\n", " rain (time_counter, y, x) float32 ...\n", " snow_ao_cea (time_counter, y, x) float32 ...\n", " snow_ai_cea (time_counter, y, x) float32 ...\n", " evap_ao_cea (time_counter, y, x) float32 ...\n", " subl_ai_cea (time_counter, y, x) float32 ...\n", " sosflxdo (time_counter, y, x) float32 ...\n", "Attributes:\n", " name: C10M-PNG-560_1m_grid_T\n", " description: Created by xios\n", " title: Created by xios\n", " Conventions: CF-1.6\n", " timeStamp: 2019-Dec-10 12:48:46 GMT\n", " uuid: 0b59c577-6754-47de-a26a-a13d8a514539\n", " LongName: IPSL 10 Ma,2X, BC from PlioMIP except for Aust...\n", " history: Wed Dec 11 05:53:06 2019: ncrcat -C --buffer_s...\n", " NCO: \"4.6.0\"\n", " nco_openmp_thread_number: 1" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ds = xr.open_dataset('./Data/C10M-PNG-560/OCE/Analyse/SE/C10M-PNG-560_SE_4655_4704_1M_grid_T.nc')\n", "ds" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Small function to find the closest ij coordinates in the curvilinear longitudes and latitudes and returns the ij coordinates.\n", "\n", "
\n", " It would be good to be able to do this with multiple longitude and latitude values but the np.substract doesn't seem to be working correctly.\n", "
" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "def find_nearest(lons, lats, lon0,lat0):\n", " idx = ((lons - lon0)**2+(lats - lat0)**2).argmin()\n", " return np.unravel_index(idx, lons.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Get the ij coordinates" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(90, 133)" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lon = 12\n", "lat = 10\n", "find_nearest(ds.nav_lon.values, ds.nav_lat.values, lon, lat)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `find_nearest` returns the indexs of the closest point so we use them with `isel` to get the value of the given variable (here tos)" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'tos' (time_counter: 12)>\n",
       "array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan],\n",
       "      dtype=float32)\n",
       "Coordinates:\n",
       "  * time_counter   (time_counter) object 0000-01-16 00:00:00 ... 0000-12-16 0...\n",
       "    nav_lat        float32 10.51\n",
       "    nav_lon        float32 11.0\n",
       "    time_centered  (time_counter) object 4679-07-16 00:00:00 ... 4680-06-16 0...\n",
       "Attributes:\n",
       "    standard_name:       sea_surface_temperature\n",
       "    long_name:           sea_surface_temperature\n",
       "    units:               degC\n",
       "    online_operation:    average\n",
       "    interval_operation:  4800 s\n",
       "    interval_write:      1 month\n",
       "    cell_methods:        time: mean (interval: 4800 s) time_counter: mean
" ], "text/plain": [ "\n", "array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan],\n", " dtype=float32)\n", "Coordinates:\n", " * time_counter (time_counter) object 0000-01-16 00:00:00 ... 0000-12-16 0...\n", " nav_lat float32 10.51\n", " nav_lon float32 11.0\n", " time_centered (time_counter) object 4679-07-16 00:00:00 ... 4680-06-16 0...\n", "Attributes:\n", " standard_name: sea_surface_temperature\n", " long_name: sea_surface_temperature\n", " units: degC\n", " online_operation: average\n", " interval_operation: 4800 s\n", " interval_write: 1 month\n", " cell_methods: time: mean (interval: 4800 s) time_counter: mean" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y, x = find_nearest(ds.nav_lon.values, ds.nav_lat.values, lon, lat)\n", "ds['tos'].isel(x=x, y=y)" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "lon: -170\tlat: 80\t value: 0.07544342428445816\n", "lon: -100\tlat: -80\t value: nan\n" ] } ], "source": [ "# Form [lon, lat]\n", "coords = [\n", " [-170, 80],\n", " [-100, -80],\n", "]\n", "\n", "values = []\n", "for coord in coords:\n", " lon, lat = coord\n", " y, x = find_nearest(ds.nav_lon.values, ds.nav_lat.values, lon, lat)\n", " value = float(ds[\"tos\"].isel(x=x, y=y).mean(\"time_counter\"))\n", " values.append(value)\n", " print(f'lon: {lon}\\tlat: {lat}\\t value: {value}')" ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[0.07544342428445816, nan]" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "values" ] } ], "metadata": { "kernelspec": { "display_name": "Python (pangeo)", "language": "python", "name": "python-pangeo" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.9" } }, "nbformat": 4, "nbformat_minor": 4 }