|
| 1 | +{ |
| 2 | + "cells": [ |
| 3 | + { |
| 4 | + "cell_type": "raw", |
| 5 | + "id": "f783c0ec-fe8a-41c3-8997-9d4a70a9a373", |
| 6 | + "metadata": {}, |
| 7 | + "source": [ |
| 8 | + "---\n", |
| 9 | + "title: \"How to get grid column and row indices for latitude and longitude\"\n", |
| 10 | + "author: Andrew P. Barrett\n", |
| 11 | + "date: last-modified\n", |
| 12 | + "---" |
| 13 | + ] |
| 14 | + }, |
| 15 | + { |
| 16 | + "cell_type": "markdown", |
| 17 | + "id": "676acb3d-0bbe-456e-b9ee-c69a3c8ed5be", |
| 18 | + "metadata": {}, |
| 19 | + "source": [ |
| 20 | + "## Problem\n", |
| 21 | + "\n", |
| 22 | + "I want to get the column and row indices of a grid that corresponds to latitude and longitude. " |
| 23 | + ] |
| 24 | + }, |
| 25 | + { |
| 26 | + "cell_type": "markdown", |
| 27 | + "id": "12a5c03b-05ea-4f1d-935a-4c302384fce5", |
| 28 | + "metadata": {}, |
| 29 | + "source": [ |
| 30 | + "## Solution\n", |
| 31 | + "\n", |
| 32 | + "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)." |
| 33 | + ] |
| 34 | + }, |
| 35 | + { |
| 36 | + "cell_type": "code", |
| 37 | + "execution_count": 5, |
| 38 | + "id": "99e45ce1-3e02-4525-8722-7056888d4839", |
| 39 | + "metadata": {}, |
| 40 | + "outputs": [ |
| 41 | + { |
| 42 | + "name": "stdout", |
| 43 | + "output_type": "stream", |
| 44 | + "text": [ |
| 45 | + "414.4236938226952 205.11385157648664\n" |
| 46 | + ] |
| 47 | + } |
| 48 | + ], |
| 49 | + "source": [ |
| 50 | + "from pyproj import Transformer\n", |
| 51 | + "from affine import Affine\n", |
| 52 | + "import numpy as np\n", |
| 53 | + "\n", |
| 54 | + "# Define the coordinate transform from WGS84 (lat, lon) \n", |
| 55 | + "# to CRS of your grid using EPSG codes.\n", |
| 56 | + "# WGS84 is EPSG:4326\n", |
| 57 | + "# NSIDC Northern Hemisphere EASE Grid is EPSG:3408\n", |
| 58 | + "transformer = Transformer.from_crs(4326, 3408)\n", |
| 59 | + "\n", |
| 60 | + "# Define the affine transformation matrix\n", |
| 61 | + "grid_cell_width = 50057.094722222224\n", |
| 62 | + "grid_cell_height = -50057.094722222224\n", |
| 63 | + "upper_left_x = -9036842.76\n", |
| 64 | + "upper_left_y = 9036842.76\n", |
| 65 | + "geotransform = Affine(grid_cell_width, 0.0, upper_left_x,\n", |
| 66 | + " 0.0, grid_cell_height, upper_left_y)\n", |
| 67 | + "\n", |
| 68 | + "x, y = transformer.transform(-45., 84.) # latitude and longitude\n", |
| 69 | + "i, j = ~geotransform * (x, y) # ~geotransform is the inverse operation\n", |
| 70 | + "\n", |
| 71 | + "print(i,j)" |
| 72 | + ] |
| 73 | + }, |
| 74 | + { |
| 75 | + "cell_type": "markdown", |
| 76 | + "id": "ca295ab6-44e3-48de-a4fa-62accf6694dd", |
| 77 | + "metadata": {}, |
| 78 | + "source": [ |
| 79 | + "`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`." |
| 80 | + ] |
| 81 | + }, |
| 82 | + { |
| 83 | + "cell_type": "code", |
| 84 | + "execution_count": 9, |
| 85 | + "id": "b3fa03fc-c38e-48e5-bff3-2d30b9f8e5d2", |
| 86 | + "metadata": {}, |
| 87 | + "outputs": [ |
| 88 | + { |
| 89 | + "name": "stdout", |
| 90 | + "output_type": "stream", |
| 91 | + "text": [ |
| 92 | + "414 205\n" |
| 93 | + ] |
| 94 | + } |
| 95 | + ], |
| 96 | + "source": [ |
| 97 | + "col, row = np.floor(i).astype(int), np.floor(j).astype(int)\n", |
| 98 | + "print(col, row)" |
| 99 | + ] |
| 100 | + }, |
| 101 | + { |
| 102 | + "cell_type": "markdown", |
| 103 | + "id": "98617747-c928-4a1f-8b63-1635d41546da", |
| 104 | + "metadata": {}, |
| 105 | + "source": [ |
| 106 | + "## Discussion\n", |
| 107 | + "\n", |
| 108 | + "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", |
| 109 | + "\n", |
| 110 | + "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", |
| 111 | + "\n", |
| 112 | + ":::{note}\n", |
| 113 | + "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", |
| 114 | + ":::\n", |
| 115 | + "\n", |
| 116 | + "TODO: Describe projected to image coordinate conversion\n", |
| 117 | + "\n", |
| 118 | + "Add image showing the relationship between the three coordinate systems" |
| 119 | + ] |
| 120 | + }, |
| 121 | + { |
| 122 | + "cell_type": "code", |
| 123 | + "execution_count": null, |
| 124 | + "id": "9d688cd2-1b28-4dba-851c-fae51a476c46", |
| 125 | + "metadata": {}, |
| 126 | + "outputs": [], |
| 127 | + "source": [] |
| 128 | + } |
| 129 | + ], |
| 130 | + "metadata": { |
| 131 | + "kernelspec": { |
| 132 | + "display_name": "Python 3 (ipykernel)", |
| 133 | + "language": "python", |
| 134 | + "name": "python3" |
| 135 | + }, |
| 136 | + "language_info": { |
| 137 | + "codemirror_mode": { |
| 138 | + "name": "ipython", |
| 139 | + "version": 3 |
| 140 | + }, |
| 141 | + "file_extension": ".py", |
| 142 | + "mimetype": "text/x-python", |
| 143 | + "name": "python", |
| 144 | + "nbconvert_exporter": "python", |
| 145 | + "pygments_lexer": "ipython3", |
| 146 | + "version": "3.12.3" |
| 147 | + } |
| 148 | + }, |
| 149 | + "nbformat": 4, |
| 150 | + "nbformat_minor": 5 |
| 151 | +} |
0 commit comments