Compute the diameter size (km) of a healpix cell

healpix_diameter_km.py

# install dependencies : pip install geopandas pandas numpy matplotlib descartes healpy

import pandas as pd
import numpy as np

import geopandas as gpd
from shapely.geometry import Polygon,MultiPolygon,Point,MultiLineString,LineString
import healpy

import matplotlib.pyplot as plt
import matplotlib
import pylab

from math import cos, asin, sqrt, pi


def distance(lat1, lon1, lat2, lon2):
    p = pi/180
    a = 0.5 - cos((lat2-lat1)*p)/2 + cos(lat1*p) * cos(lat2*p) * (1-cos((lon2-lon1)*p))/2
    return 12742 * asin(sqrt(a)) #2*R*asin...


def latlon2healpix( lat , lon , nside ):
    lat = np.radians(lat)
    lon = np.radians(lon)
    xs = ( np.cos(lat) * np.cos(lon) ) #
    ys = ( np.cos(lat) * np.sin(lon) ) # -> Sphere coordinates: https://vvvv.org/blog/polar-spherical-and-geographic-coordinates
    zs = ( np.sin(lat) ) # 
    return healpy.vec2pix( int(nside) , xs , ys , zs )


coord_div = 20 # precision (1 = 1°, 4 = 0.25°)
nside= 128 # check number of cell formula : Nb_cell = 12 * (nside^2)

LAT, LON = 180*coord_div,360*coord_div
mat = np.zeros((LAT,LON))


for i in range(LAT):
    for j in range(LON):
        mat[i,j] = latlon2healpix((i/coord_div)-90,(j/coord_div)-180,nside)

# A random colormap for matplotlib
cmap = matplotlib.colors.ListedColormap ( np.random.rand ( 12*(nside**2),3))


def foo(x,y):
    x,y = np.asarray(x), np.asarray(y)
    x+=180
    x*=coord_div
    y+=90
    y*=coord_div
    X = np.concatenate((x.reshape(-1,1),y.reshape(-1,1)),axis=1)
    return X

def transform_polygon(boundary):
    shapes_array=[]
    if isinstance(boundary,MultiLineString):
        for b in boundary:
            shapes_array.append(LineString(foo(*b.xy)))
        return Polygon(b)
    else:
        return Polygon(foo(*boundary.xy))

def transform_multi_polygon(g):
    new_poly = [transform_polygon(boundary) for boundary in g.boundary]
    return MultiPolygon(new_poly)
        

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
world["geometry"] = world.geometry.apply(lambda x: transform_polygon(x.boundary) if isinstance(x,Polygon) else transform_multi_polygon(x))


fig, ax = plt.subplots(1,figsize=(20,10))
ax.imshow(mat,cmap=cmap)
world.boundary.plot(ax=ax,color="black")
ax.set_ylim(0,180*coord_div)
plt.axis("off")
plt.show()


for i in range(50):
    hp_id = np.random.choice(mat.flatten(),1)[0]
    x_min,x_max,y_min,y_max=180,-180,90,-90
    for i in range(LAT):
        for j in range(LON):
            if mat[i,j] == hp_id:
                x=(j/coord_div)-180
                y=(i/coord_div) - 90
                x_min=min(x_min,x)
                x_max=max(x_max,x)
                y_min=min(y_min,y)
                y_max=max(y_max,y)
    print(x_min,y_min,x_max,y_max)
    print(distance(y_min,x_min,y_max,x_max))