docs/examples/creating_and_working_with_grid_objects.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pygeogrids.grids as grids\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's create a simple regular 10x10 degree grid with grid points at the center of each 10x10 degree cell.\n",
    "\n",
    "First by hand to understand what is going on underneath"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-175 -165 -155 -145 -135 -125 -115 -105  -95  -85  -75  -65  -55  -45\n",
      "  -35  -25  -15   -5    5   15   25   35   45   55   65   75   85   95\n",
      "  105  115  125  135  145  155  165  175]\n",
      "[ 85  75  65  55  45  35  25  15   5  -5 -15 -25 -35 -45 -55 -65 -75 -85]\n"
     ]
    }
   ],
   "source": [
    "# create the longitudes\n",
    "lons = np.arange(-180 + 5, 180, 10)\n",
    "print(lons)\n",
    "lats = np.arange(90 - 5, -90, -10)\n",
    "print(lats)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are just the dimensions or we can also call them the \"sides\" of the array that defines all the gridpoints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create all the grid points by using the numpy.meshgrid function\n",
    "longrid, latgrid = np.meshgrid(lons, lats)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "now we can create a BasicGrid. We can also define the shape of the grid. The first part of the shape must be in longitude direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 1 2 3 4 5 6 7 8 9] [-175 -165 -155 -145 -135 -125 -115 -105  -95  -85] [85 85 85 85 85 85 85 85 85 85]\n"
     ]
    }
   ],
   "source": [
    "manualgrid = grids.BasicGrid(longrid.flatten(), latgrid.flatten(), shape=(18, 36))\n",
    "\n",
    "# Each point of the grid automatically got a grid point number\n",
    "gpis, gridlons, gridlats = manualgrid.get_grid_points()\n",
    "print(gpis[:10], gridlons[:10], gridlats[:10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The grid point indices or numbers are useful when creating lookup tables between grids."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now use the manualgrid instance to find the nearest gpi to any longitude and latitude"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "235 424808.5131778209\n",
      "(15, 25)\n"
     ]
    }
   ],
   "source": [
    "ngpi, distance = manualgrid.find_nearest_gpi(15.84, 28.76)\n",
    "print(ngpi, distance)\n",
    "# convert the gpi to longitude and latitude\n",
    "print(manualgrid.gpi2lonlat(ngpi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The same grid can also be created by a method for creating regular grids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "autogrid = grids.genreg_grid(10, 10)\n",
    "autogrid == manualgrid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If your grid has a 2D shape like the ones we just created then you can also get the row and the column of a grid point.\n",
    "This can be useful if you know that you have data stored on a specific grid and you want to read the data from a grid point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6 19\n"
     ]
    }
   ],
   "source": [
    "row, col = autogrid.gpi2rowcol(ngpi)\n",
    "print(row, col)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Iteration over gridpoints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 -175.0 85.0\n",
      "1 -165.0 85.0\n",
      "2 -155.0 85.0\n",
      "3 -145.0 85.0\n",
      "4 -135.0 85.0\n",
      "5 -125.0 85.0\n",
      "6 -115.0 85.0\n",
      "7 -105.0 85.0\n",
      "8 -95.0 85.0\n",
      "9 -85.0 85.0\n",
      "10 -75.0 85.0\n"
     ]
    }
   ],
   "source": [
    "for i, (gpi, lon, lat) in enumerate(autogrid.grid_points()):\n",
    "    print(gpi, lon, lat)\n",
    "    if i==10: # this is just to keep the example output short\n",
    "        break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculation of lookup tables\n",
    "\n",
    "If you have a two grids and you know that you want to get the nearest neighbors for all of its grid points in the second grid you can calculate a lookup table once and reuse it later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-30.89674861  59.08649779 -49.27545779 -69.52856569  -1.39290501\n",
      " -10.88163146  53.67514799  60.14998259  55.15218786 -29.79836856]\n",
      "[ 125.68580631 -169.379169    105.7270734    46.42576672  -28.28325787\n",
      "   47.49761201  173.67561423   86.59161319 -178.39301312  115.68242529]\n"
     ]
    }
   ],
   "source": [
    "# lets generate a second grid with 10 random points on the Earth surface.\n",
    "\n",
    "randlat = np.random.random(10) * 180 - 90\n",
    "randlon = np.random.random(10) * 360 - 180\n",
    "print(randlat)\n",
    "print(randlon)\n",
    "# This grid has no meaningful 2D shape so none is given\n",
    "randgrid = grids.BasicGrid(randlon, randlat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets calculate a lookup table to the regular 10x10° grid we created earlier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[462 109 496 562 339 382 143  98 108 425]\n"
     ]
    }
   ],
   "source": [
    "lut = randgrid.calc_lut(autogrid)\n",
    "print(lut)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The lookup table contains the grid point indices of the other grid, autogrid in this case. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-35.  55. -45. -65.  -5. -15.  55.  65.  55. -25.]\n",
      "[ 125. -165.  105.   45.  -25.   45.  175.   85. -175.  115.]\n"
     ]
    }
   ],
   "source": [
    "lut_lons, lut_lats = autogrid.gpi2lonlat(lut)\n",
    "print(lut_lats)\n",
    "print(lut_lons)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Storing and loading grids\n",
    "Grids can be stored to disk as CF compliant netCDF files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pygeogrids.netcdf as nc\n",
    "nc.save_grid('example.nc', randgrid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "loadedgrid = nc.load_grid('example.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<pygeogrids.grids.BasicGrid at 0x7f4110fe7700>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "loadedgrid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<pygeogrids.grids.BasicGrid at 0x7f41128d92b0>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "randgrid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define geodetic datum for grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid_WGS84 = grids.BasicGrid(randlon, randlat, geodatum='WGS84')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid_GRS80 = grids.BasicGrid(randlon, randlat, geodatum='GRS80')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6378137.0"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid_WGS84.geodatum.geod.a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6378137.0"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid_GRS80.geodatum.geod.a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid_WGS84.kdTree.geodatum.geod.sphere"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}