diff --git a/WGS84_to_HEALPix_Conversion.ipynb b/WGS84_to_HEALPix_Conversion.ipynb new file mode 100644 index 0000000..dbeb8e9 --- /dev/null +++ b/WGS84_to_HEALPix_Conversion.ipynb @@ -0,0 +1,287 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "6fb64788", + "metadata": {}, + "source": [ + "\n", + "# WGS84 to HEALPix Conversion Using Python\n", + "\n", + "In this notebook, we will demonstrate how to convert geographic coordinates (latitude, longitude) based on the WGS84 reference system into HEALPix pixel indices. This involves transforming coordinates into a geocentric Cartesian system and subsequently into spherical coordinates that are compatible with HEALPix.\n", + "\n", + "---\n", + "\n", + "## What is WGS84?\n", + "\n", + "WGS84 (World Geodetic System 1984) is a global reference system for geographic coordinates (latitude, longitude, and altitude). It is widely used for mapping, GPS, and geospatial data applications.\n", + "\n", + "For more details, you can refer to the [WGS84 definition](https://en.wikipedia.org/wiki/World_Geodetic_System).\n", + "\n", + "---\n", + "\n", + "## What is HEALPix?\n", + "For more information on HEALPix, visit the [HEALPix Documentation](https://healpix.sourceforge.io/).\n", + "\n", + "\n", + "---\n", + "\n", + "## Setup\n", + "\n", + "To run this notebook, you need to install the following libraries:\n", + "1. **Pyproj**: For coordinate transformations.\n", + "2. **Healpy**: For working with HEALPix, a hierarchical spatial indexing system.\n", + "3. **XDGGS**: For working with HEALPix together with Xarray to handle simpler data manipulation.\n", + "\n", + "Additionally, if you are using XDGGS for spatial data applications, you will need to install the XDGGS library along with Lonboard. You can find the repository here: [XDGGS GitHub Repository](https://github.com/xdggs/xdggs).\n", + "\n", + "To install the necessary dependencies, use the following command:\n", + "\n", + "```bash\n", + "pip install pyproj healpy xdggs\n", + "\n" + ] + }, + { + "cell_type": "raw", + "id": "a197793c-552b-4973-a124-8482ca50e2a5", + "metadata": { + "scrolled": true + }, + "source": [ + "pip install pyproj healpy xdggs" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "a54d14d4", + "metadata": {}, + "outputs": [], + "source": [ + "# Import necessary libraries\n", + "import pyproj\n", + "import numpy as np\n", + "import healpy as hp # HEALPix library\n" + ] + }, + { + "cell_type": "markdown", + "id": "0f66516b", + "metadata": {}, + "source": [ + "## Initialize Coordinate Reference Systems (CRS)\n", + "\n", + "Here, we initialize two CRS:\n", + "1. **WGS84**: A geographic CRS (latitude, longitude).\n", + "2. **Geocentric Cartesian**: A 3D Cartesian system for transformations involving Earth-centered coordinates." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "a7e56800", + "metadata": {}, + "outputs": [], + "source": [ + "# Initialize CRS (Coordinate Reference Systems)\n", + "wgs84 = pyproj.CRS(\"EPSG:4326\") # WGS84 (geographic latitude, longitude)\n", + "geocentric = pyproj.CRS(\"EPSG:4978\") # Geocentric Cartesian coordinate system (based on WGS84)" + ] + }, + { + "cell_type": "markdown", + "id": "648efe43", + "metadata": {}, + "source": [ + "## Create a Transformer for WGS84 to Geocentric Conversion\n", + "\n", + "Using `pyproj`, we create a transformer that converts geographic coordinates (latitude, longitude) into a geocentric Cartesian system. This is necessary for further spherical transformations." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "f8f80d0e", + "metadata": {}, + "outputs": [], + "source": [ + "# Create a transformer to convert WGS84 -> geocentric\n", + "transformer = pyproj.Transformer.from_crs(wgs84, geocentric, always_xy=True)" + ] + }, + { + "cell_type": "markdown", + "id": "28d4a81d", + "metadata": {}, + "source": [ + "## Define Function: WGS84 to Spherical Coordinates\n", + "\n", + "This function converts WGS84 geographic coordinates into HEALPix-compatible spherical coordinates (theta, phi). These coordinates are essential for identifying the correct HEALPix pixel index." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "3fec33f0", + "metadata": {}, + "outputs": [], + "source": [ + "# Define a function to convert WGS84 (lat, lon) to spherical coordinates (HEALPix compatible)\n", + "def wgs84_to_spherical(lat, lon):\n", + " \"\"\"\n", + " Converts geographic coordinates (latitude, longitude) to HEALPix-compatible spherical coordinates.\n", + " \n", + " Parameters:\n", + " lat (float): Latitude in degrees.\n", + " lon (float): Longitude in degrees.\n", + " \n", + " Returns:\n", + " tuple: Spherical coordinates (theta, phi) in radians.\n", + " \"\"\"\n", + " # Convert to geocentric Cartesian coordinates\n", + " x, y, z = transformer.transform(lon, lat, 0) # Altitude = 0 (on the ellipsoid)\n", + "\n", + " # Convert to spherical angles\n", + " r = np.sqrt(x**2 + y**2 + z**2) # Radius\n", + " theta = np.arccos(z / r) # Co-latitude in radians\n", + " phi = np.arctan2(y, x) # Longitude in radians\n", + " return theta, phi" + ] + }, + { + "cell_type": "markdown", + "id": "c18f3dfe", + "metadata": {}, + "source": [ + "## Define Function: HEALPix Pixel Index\n", + "\n", + "This function calculates the HEALPix pixel index for a given set of geographic coordinates using the resolution parameter `nside`." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "3a540902", + "metadata": {}, + "outputs": [], + "source": [ + "# Define a function to calculate the HEALPix pixel index\n", + "def ang2pix_wgs84(lat, lon, nside):\n", + " \"\"\"\n", + " Converts WGS84 coordinates to a HEALPix pixel index.\n", + " \n", + " Parameters:\n", + " lat (float): Latitude in degrees.\n", + " lon (float): Longitude in degrees.\n", + " nside (int): HEALPix resolution parameter.\n", + " \n", + " Returns:\n", + " int: Pixel index corresponding to the given coordinates.\n", + " \"\"\"\n", + " # Convert WGS84 to spherical coordinates\n", + " theta, phi = wgs84_to_spherical(lat, lon)\n", + " \n", + " # Return the HEALPix pixel index\n", + " return hp.ang2pix(nside, theta, phi)" + ] + }, + { + "cell_type": "markdown", + "id": "7284e4e2", + "metadata": {}, + "source": [ + "## Example Usage\n", + "\n", + "Finally, we use the functions defined above to calculate the HEALPix pixel index for a specific location, in this case, Paris." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "51aa92d4", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Pixel index: 5941\n" + ] + } + ], + "source": [ + "# Example usage\n", + "nside = 64 # HEALPix resolution\n", + "lat, lon = 48.8566, 2.3522 # Latitude and longitude of Paris\n", + "\n", + "# Get the HEALPix pixel index for the given coordinates\n", + "pix = ang2pix_wgs84(lat, lon, nside)\n", + "\n", + "# Print the result\n", + "print(f\"Pixel index: {pix}\")" + ] + }, + { + "cell_type": "markdown", + "id": "b0e09199", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "## Conclusion\n", + "\n", + "This notebook demonstrates the process of converting WGS84 coordinates to HEALPix indices. The workflow involves:\n", + "1. Transforming WGS84 coordinates into a geocentric Cartesian system.\n", + "2. Converting the Cartesian coordinates into spherical coordinates.\n", + "3. Calculating the HEALPix pixel index using the spherical coordinates.\n", + "\n", + "\n", + "---" + ] + }, + { + "cell_type": "markdown", + "id": "fc28efe3-98d0-445f-925c-2b18c06a9ccb", + "metadata": {}, + "source": [ + "## Next Step\n", + "\n", + "We will show:\n", + "1. The transformation does not loose the benefit specific of HealPIx, such as spherical harmonics properties. \n", + "2. Comparision with rHealPix for err geolocalisation. \n", + "\n", + "---" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "70c863c0-ab05-455a-89a8-6f6959778dcf", + "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.5" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}