Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add JM's WGS84_to_HEALPix_Conversion #4

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

add JM's WGS84_to_HEALPix_Conversion #4

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
287 changes: 287 additions & 0 deletions WGS84_to_HEALPix_Conversion.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
}