Geographic to image coords howto

_quarto.yml

@@ -24,7 +24,25 @@ website:
     collapse-level: 1
     contents:
       - text: "Welcome"
-        href: nsidc-data-cookbook/introduction.qmd
+        href: index.qmd
+      - text: "Our Cookbook"
+        href: nsidc-data-cookbook/our-cookbook.qmd
+      - text: "Tools we use"
+        href: tools-we-use.qmd
+      - text: "Data Structures - types of data"
+        href: reference-guides/data-structures.qmd
+      - text: "Dataset search and discovery"
+        href: reference-guides/dataset-search-and-discovery.qmd
+      - text: "Common NSIDC file formats"
+        href: reference-guides/nsidc-file-formats.qmd
+      - text: "NSIDC map projections and grids"
+        href: nsidc-data-cookbook/reference-guides/projections.qmd
+      - text: "Plotting NSIDC Data"
+        href: nsidc-data-cookbook/reference-guides/plotting-data.qmd
+      - text: "Datasets"
+        href: nsidc-data-cookbook/reference-guides/datasets.qmd
+      - text: "What is an API?"
+      - text: "NASA collections and granules"
       - section: "How do I..."
         href: nsidc-data-cookbook/how-to-guides/overview.qmd
         contents:
@@ -35,17 +53,6 @@ website:
               href: nsidc-data-cookbook/how-to-guides/netcdf_cf.qmd
             - text: "get the latitude and longitudes for grid cells"
               href: nsidc-data-cookbook/how-to-guides/get_latitude_and_longitude.qmd
-      - section: "Tutorials"
-        contents:
-      - section: "Reference Guides"
-        href: nsidc-data-cookbook/reference-guides/overview.qmd
-        contents:
-          - text: "What is an API?"
-          - text: "NASA collections and granules"
-          - text: "Map projections and grids"
-            href: nsidc-data-cookbook/reference-guides/projections.qmd
-      - section: "Examples"
-        contents:
 
 format:
   html:

nsidc-data-cookbook/how-to-guides/get_column_and_row_indices.ipynb

@@ -0,0 +1,151 @@
+{
+ "cells": [
+  {
+   "cell_type": "raw",
+   "id": "f783c0ec-fe8a-41c3-8997-9d4a70a9a373",
+   "metadata": {},
+   "source": [
+    "---\n",
+    "title: \"How to get grid column and row indices for latitude and longitude\"\n",
+    "author: Andrew P. Barrett\n",
+    "date: last-modified\n",
+    "---"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "676acb3d-0bbe-456e-b9ee-c69a3c8ed5be",
+   "metadata": {},
+   "source": [
+    "## Problem\n",
+    "\n",
+    "I want to get the column and row indices of a grid that corresponds to latitude and longitude. "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "12a5c03b-05ea-4f1d-935a-4c302384fce5",
+   "metadata": {},
+   "source": [
+    "## Solution\n",
+    "\n",
+    "In Python you can use the `pyproj` package to transform latitude and longitude (geographic coordinates) to x and y (projected coordinates) in the EASE Grid CRS, and then use the `affine` package to get the column and row indices (`(i,j)` image coordinates)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "99e45ce1-3e02-4525-8722-7056888d4839",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "414.4236938226952 205.11385157648664\n"
+     ]
+    }
+   ],
+   "source": [
+    "from pyproj import Transformer\n",
+    "from affine import Affine\n",
+    "import numpy as np\n",
+    "\n",
+    "# Define the coordinate transform from WGS84 (lat, lon) \n",
+    "# to CRS of your grid using EPSG codes.\n",
+    "# WGS84 is EPSG:4326\n",
+    "# NSIDC Northern Hemisphere EASE Grid is EPSG:3408\n",
+    "transformer = Transformer.from_crs(4326, 3408)\n",
+    "\n",
+    "# Define the affine transformation matrix\n",
+    "grid_cell_width = 50057.094722222224\n",
+    "grid_cell_height = -50057.094722222224\n",
+    "upper_left_x = -9036842.76\n",
+    "upper_left_y = 9036842.76\n",
+    "geotransform = Affine(grid_cell_width, 0.0, upper_left_x,\n",
+    "                      0.0, grid_cell_height, upper_left_y)\n",
+    "\n",
+    "x, y = transformer.transform(-45., 84.)  # latitude and longitude\n",
+    "i, j = ~geotransform * (x, y)  # ~geotransform is the inverse operation\n",
+    "\n",
+    "print(i,j)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ca295ab6-44e3-48de-a4fa-62accf6694dd",
+   "metadata": {},
+   "source": [
+    "`i` and `j` are in image coordinates with the origin as the upper-left corner of the upper-left grid-cell.  To get the column and row indices of the cell you can use `floor`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "id": "b3fa03fc-c38e-48e5-bff3-2d30b9f8e5d2",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "414 205\n"
+     ]
+    }
+   ],
+   "source": [
+    "col, row = np.floor(i).astype(int), np.floor(j).astype(int)\n",
+    "print(col, row)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "98617747-c928-4a1f-8b63-1635d41546da",
+   "metadata": {},
+   "source": [
+    "## Discussion\n",
+    "\n",
+    "Many grids, whether they are images or gridded data are in _projected coordinate systems_ so that they can be displayed on a 2D screen.  Converting from geographic coordinates (_latitude_, _longitude_) to image coordinates (_column_, _row_) is usually a two step process.  In the first step, the geographic coordinates are transformed from a Geographic Coordinate System with coordinates _latitude_ and _longitude_ to the _projected coordinate system_ of the image or grid.  This is a _Cartesian_ coordinate system with coordinates _x_ and _y_, usually in meters.  In the second step, the _projected coordinates_ are converted into _image coordinates_.  This is another _Cartesian_ coordinate system with coordinates _column_ and _row_ or _i_ and _j_ with units grid-cells.  Image coordinate systems often have the origin in the upper-left corner of the upper-left grid cell.\n",
+    "\n",
+    "In Python the first step is easily acheived using the `pyproj` package.  If the Coordinate Reference System (CRS) of the grid is registered in the EPSG database, defining the transformation is easily done with the EPSG codes for the CRS for the _latitude_ and _longitude_ coordinates and the CRS of the destination grid.  Usually, _latitude_ and _longitude_ are in the Ensemble WGS84 CRS, which is EPSG:4326.  For NSIDC gridded data, the CRS is some flavour of EASE Grid or Polar Stereographic.  See the help pages for [EASE Grid](https://nsidc.org/data/user-resources/help-center/guide-ease-grids) and [Polar Stereographic](https://nsidc.org/data/user-resources/help-center/guide-nsidcs-polar-stereographic-projection) for more information. \n",
+    "\n",
+    ":::{note}\n",
+    "Both EASE Grid and Polar Stereographic CRS come in two main flavours _Original_ and WGS84.  For EASE Grids, these two flavours are EASE-Grid and EASE-Grid v2.0.\n",
+    ":::\n",
+    "\n",
+    "TODO: Describe projected to image coordinate conversion\n",
+    "\n",
+    "Add image showing the relationship between the three coordinate systems"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9d688cd2-1b28-4dba-851c-fae51a476c46",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "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.12.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

our-cookbook.qmd

@@ -0,0 +1,11 @@
+---
+title: "Our Cookbook"
+---
+
+## How to use this cookbook
+
+
+## Get Started
+
+
+## Tools we use