04 healpix

Working with HEALPix in Vgrid DGGS¶

Full Vgrid DGGS documentation is available at vgrid document.

To work with Vgrid DGGS directly in GeoPandas and Pandas, please use vgridpandas. Full Vgridpandas DGGS documentation is available at vgridpandas document.

To work with Vgrid DGGS in QGIS, install the Vgrid Plugin.

To visualize DGGS in Maplibre GL JS, try the vgrid-maplibre library.

For an interactive demo, visit the Vgrid Homepage.

Install vgrid¶

Uncomment the following line to install vgrid.

In [ ]:

# %pip install vgrid --upgrade

# %pip install vgrid --upgrade

latlon2healpix¶

In [ ]:

from vgrid.conversion.latlon2dggs import latlon2healpix

lat = 10.775276
lon = 180.0
res = 0
healpix_id = latlon2healpix(lat, lon, res)
healpix_id

from vgrid.conversion.latlon2dggs import latlon2healpix lat = 10.775276 lon = 180.0 res = 0 healpix_id = latlon2healpix(lat, lon, res) healpix_id

Out[ ]:

5

HEALPix to Polygon¶

In [2]:

from vgrid.conversion.dggs2geo.healpix2geo import healpix2geo

healpix_geo = healpix2geo(healpix_id)
healpix_geo

from vgrid.conversion.dggs2geo.healpix2geo import healpix2geo healpix_geo = healpix2geo(healpix_id) healpix_geo

Out[2]:

HEALPix to GeoJSON¶

In [3]:

from vgrid.conversion.dggs2geo.healpix2geo import healpix2geojson

healpix_geojson = healpix2geojson(healpix_id)
healpix_geojson

from vgrid.conversion.dggs2geo.healpix2geo import healpix2geojson healpix_geojson = healpix2geojson(healpix_id) healpix_geojson

Out[3]:

{'type': 'FeatureCollection',
 'features': [{'type': 'Feature',
   'geometry': {'type': 'Polygon',
    'coordinates': (((90.0, 41.81031489577859),
      (45.0, 0.0),
      (90.0, -41.81031489577861),
      (134.99999999999994, 0.0),
      (90.0, 41.81031489577859)),)},
   'properties': {'healpix': '5',
    'resolution': 0,
    'center_lat': -7.105427357601001e-15,
    'center_lon': 89.99999999999997,
    'avg_edge_len': 6462621.584028445,
    'cell_area': 50810264447065.44,
    'cell_perimeter': 25850486.33611378}}]}

Vector to HEALPix¶

(Multi)Point to HEALPix¶

In [4]:

from vgrid.conversion.vector2dggs.vector2healpix import vector2healpix
import geopandas as gpd

file_path = (
    "https://raw.githubusercontent.com/opengeoshub/vopendata/main/shape/point.geojson"
)

gdf = gpd.read_file(file_path)

point_to_healpix = vector2healpix(
    gdf,
    # resolution=9,
    topology=True,
    output_format="gpd",
)
point_to_healpix.plot()

from vgrid.conversion.vector2dggs.vector2healpix import vector2healpix import geopandas as gpd file_path = ( "https://raw.githubusercontent.com/opengeoshub/vopendata/main/shape/point.geojson" ) gdf = gpd.read_file(file_path) point_to_healpix = vector2healpix( gdf, # resolution=9, topology=True, output_format="gpd", ) point_to_healpix.plot()

Processing features: 100%|██████████| 12/12 [00:00<00:00, 3906.52it/s]

Out[4]:

<Axes: >

(Multi)Polyline to HEALPix¶

In [8]:

from vgrid.conversion.vector2dggs.vector2healpix import vector2healpix

file_path = (
    "https://raw.githubusercontent.com/opengeoshub/vopendata/main/shape/polyline2.geojson"
)
polyline_to_healpix = vector2healpix(
    file_path,
    resolution=18,
    output_format="gpd",
)
polyline_to_healpix.plot(edgecolor="white")

from vgrid.conversion.vector2dggs.vector2healpix import vector2healpix file_path = ( "https://raw.githubusercontent.com/opengeoshub/vopendata/main/shape/polyline2.geojson" ) polyline_to_healpix = vector2healpix( file_path, resolution=18, output_format="gpd", ) polyline_to_healpix.plot(edgecolor="white")

Processing features: 100%|██████████| 1/1 [03:55<00:00, 235.19s/it]

Out[8]:

<Axes: >

(Multi)Polygon to HEALPix¶

In [ ]:

from vgrid.conversion.vector2dggs.vector2healpix import vector2healpix

file_path = (
    "https://raw.githubusercontent.com/opengeoshub/vopendata/main/shape/polygon2.geojson"
)
vector_to_healpix = vector2healpix(
    file_path,
    resolution=18,
    compact=False,
    predicate="intersects",
    output_format="gpd",
)
vector_to_healpix.plot(edgecolor="white")

from vgrid.conversion.vector2dggs.vector2healpix import vector2healpix file_path = ( "https://raw.githubusercontent.com/opengeoshub/vopendata/main/shape/polygon2.geojson" ) vector_to_healpix = vector2healpix( file_path, resolution=18, compact=False, predicate="intersects", output_format="gpd", ) vector_to_healpix.plot(edgecolor="white")

Processing features:   0%|          | 0/1 [00:00<?, ?it/s]

HEALPix Compact¶

In [ ]:

from vgrid.conversion.dggscompact.healpixcompact import healpixcompact

healpix_compacted = healpixcompact(
    vector_to_healpix, healpix_id="healpix", output_format="gpd"
)
healpix_compacted.plot(edgecolor="white")

from vgrid.conversion.dggscompact.healpixcompact import healpixcompact healpix_compacted = healpixcompact( vector_to_healpix, healpix_id="healpix", output_format="gpd" ) healpix_compacted.plot(edgecolor="white")

HEALPix Expand¶

In [ ]:

from vgrid.conversion.dggscompact.healpixcompact import healpixexpand

healpix_expanded = healpixexpand(
    vector_to_healpix, resolution=11, output_format="gpd"
)
healpix_expanded.plot(edgecolor="white")

from vgrid.conversion.dggscompact.healpixcompact import healpixexpand healpix_expanded = healpixexpand( vector_to_healpix, resolution=11, output_format="gpd" ) healpix_expanded.plot(edgecolor="white")

Raster to HEALPix¶

In [13]:

from vgrid.utils.io import download_file
import rasterio
from rasterio.plot import show

raster_url = (
    "https://raw.githubusercontent.com/opengeoshub/vopendata/main/raster/rgb.tif"
)
raster_file = download_file(raster_url)
src = rasterio.open(raster_file, "r")
print(src.meta)
show(src)

from vgrid.utils.io import download_file import rasterio from rasterio.plot import show raster_url = ( "https://raw.githubusercontent.com/opengeoshub/vopendata/main/raster/rgb.tif" ) raster_file = download_file(raster_url) src = rasterio.open(raster_file, "r") print(src.meta) show(src)

rgb.tif already exists. Skip downloading. Set overwrite=True to overwrite.
{'driver': 'GTiff', 'dtype': 'uint8', 'nodata': None, 'width': 240, 'height': 147, 'count': 3, 'crs': CRS.from_wkt('GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]'), 'transform': Affine(2.6640125000199077e-06, 0.0, 106.708118755,
       0.0, -2.6640136054383103e-06, 10.812568272)}

Out[13]:

<Axes: >

In [23]:

from vgrid.conversion.raster2dggs.raster2healpix import raster2healpix


raster_to_healpix = raster2healpix(
    raster_file,
    method="binning",  # nearest, binning
    stats="mean",
    output_format="gpd",
)
# Visualize the output
import folium

m = folium.Map(tiles="CartoDB positron", max_zoom=28)

healpix_layer = folium.GeoJson(
    raster_to_healpix,
    style_function=lambda x: {
        "fillColor": f"rgb({x['properties']['band_1']}, {x['properties']['band_2']}, {x['properties']['band_3']})",
        "fillOpacity": 1,
        "color": "black",
        "weight": 1,
    },
    popup=folium.GeoJsonPopup(
        fields=["healpix", "resolution", "band_1", "band_2", "band_3", "cell_area"],
        aliases=[
            "HEALPix ID",
            "Resolution",
            "Band 1",
            "Band 2",
            "Band 3",
            "Area (m²)",
        ],
        style="""
            background-color: white;
            border: 2px solid black;
            border-radius: 3px;
            box-shadow: 3px;
        """,
    ),
).add_to(m)

m.fit_bounds(healpix_layer.get_bounds())

# Display the map
m

from vgrid.conversion.raster2dggs.raster2healpix import raster2healpix raster_to_healpix = raster2healpix( raster_file, method="binning", # nearest, binning stats="mean", output_format="gpd", ) # Visualize the output import folium m = folium.Map(tiles="CartoDB positron", max_zoom=28) healpix_layer = folium.GeoJson( raster_to_healpix, style_function=lambda x: { "fillColor": f"rgb({x['properties']['band_1']}, {x['properties']['band_2']}, {x['properties']['band_3']})", "fillOpacity": 1, "color": "black", "weight": 1, }, popup=folium.GeoJsonPopup( fields=["healpix", "resolution", "band_1", "band_2", "band_3", "cell_area"], aliases=[ "HEALPix ID", "Resolution", "Band 1", "Band 2", "Band 3", "Area (m²)", ], style=""" background-color: white; border: 2px solid black; border-radius: 3px; box-shadow: 3px; """, ), ).add_to(m) m.fit_bounds(healpix_layer.get_bounds()) # Display the map m

Cell size: 0.08638527081938627 m2
Nearest HEALPix resolution determined: 24
Method: binning
Stats: mean

Binning raster blocks to HEALPix: 100%|██████████| 14/14 [00:00<00:00, 21.66 block/s]
Converting raster to HEALPix: 100%|██████████| 19948/19948 [00:02<00:00, 9379.17 cells/s]

Out[23]:

Make this Notebook Trusted to load map: File -> Trust Notebook

HEALPix Generator¶

In [26]:

from vgrid.generator.healpixgrid import healpixgrid

healpix_grid = healpixgrid(
    resolution=2,
    # bbox=[106.699007, 10.762811, 106.717674, 10.778649],
    compact=False,
    output_format="gpd",
)
healpix_grid.plot(edgecolor="white")

from vgrid.generator.healpixgrid import healpixgrid healpix_grid = healpixgrid( resolution=2, # bbox=[106.699007, 10.762811, 106.717674, 10.778649], compact=False, output_format="gpd", ) healpix_grid.plot(edgecolor="white")

Generating HEALPix DGGS: 100%|██████████| 192/192 [00:00<00:00, 8403.05 cells/s]

Out[26]:

<Axes: >

HEALPix Inspect¶

In [ ]:

from vgrid.stats.healpixstats import healpixinspect

resolution = 4
healpix_inspect = healpixinspect(resolution)
healpix_inspect.head()

from vgrid.stats.healpixstats import healpixinspect resolution = 4 healpix_inspect = healpixinspect(resolution) healpix_inspect.head()

Generating HEALPix DGGS: 100%|██████████| 3072/3072 [00:00<00:00, 10312.30 cells/s]

Out[ ]:

	healpix	resolution	center_lat	center_lon	avg_edge_len	cell_area	cell_perimeter	geometry	crossed	norm_area	ipq	zsc	cvh
0	1020	4	2.389402	45.0000	409446.133674	1.654049e+11	1.637785e+06	POLYGON ((45 4.78019, 42.1875 2.38802, 45 0, 4...	False	0.996242	0.774898	0.880141	1.0
1	1021	4	4.782988	47.8125	409283.696729	1.654171e+11	1.637135e+06	POLYGON ((47.8125 7.18076, 45 4.78019, 47.8125...	False	0.996316	0.775571	0.880522	1.0
2	1022	4	4.782988	42.1875	409283.696729	1.654171e+11	1.637135e+06	POLYGON ((42.1875 7.18076, 39.375 4.78019, 42....	False	0.996316	0.775571	0.880522	1.0
3	1023	4	7.185005	45.0000	409023.550892	1.654374e+11	1.636094e+06	POLYGON ((45 9.59407, 42.1875 7.18076, 45 4.78...	False	0.996438	0.776653	0.881136	1.0
4	1024	4	7.185005	50.6250	409023.550892	1.654374e+11	1.636094e+06	POLYGON ((50.625 9.59407, 47.8125 7.18076, 50....	False	0.996438	0.776653	0.881136	1.0

HEALPix Normalized Area Histogram¶

In [15]:

from vgrid.stats.healpixstats import healpix_norm_area_hist

healpix_norm_area_hist(healpix_inspect)

from vgrid.stats.healpixstats import healpix_norm_area_hist healpix_norm_area_hist(healpix_inspect)

Distribution of HEALPix Area Distortions¶

In [16]:

from vgrid.stats.healpixstats import healpix_norm_area

healpix_norm_area(healpix_inspect)

from vgrid.stats.healpixstats import healpix_norm_area healpix_norm_area(healpix_inspect)

HEALPix IPQ Compactness Histogram¶

Isoperimetric Inequality (IPQ) Compactness (suggested by Osserman, 1978):

$$C_{IPQ} = \frac{4 \pi A}{p^2}$$ The range of the IPQ compactness metric is (0,1].

A circle represents the maximum compactness with a value of 1.

As shapes become more irregular or elongated, their compactness decreases toward 0.

In [17]:

from vgrid.stats.healpixstats import healpix_compactness_ipq_hist

healpix_compactness_ipq_hist(healpix_inspect)

from vgrid.stats.healpixstats import healpix_compactness_ipq_hist healpix_compactness_ipq_hist(healpix_inspect)

Distribution of HEALPix IPQ Compactness¶

In [18]:

from vgrid.stats.healpixstats import healpix_compactness_ipq

healpix_compactness_ipq(healpix_inspect)

from vgrid.stats.healpixstats import healpix_compactness_ipq healpix_compactness_ipq(healpix_inspect)

HEALPix Convex hull Compactness Histogram:¶

$$C_{CVH} = \frac{A}{A_{CVH}}$$

The range of the convex hull compactness metric is (0,1].

As shapes become more concave, their convex hull compactness decreases toward 0.

In [19]:

from vgrid.stats.healpixstats import healpix_compactness_cvh_hist

healpix_compactness_cvh_hist(healpix_inspect)

from vgrid.stats.healpixstats import healpix_compactness_cvh_hist healpix_compactness_cvh_hist(healpix_inspect)

Distribution of HEALPix Convex hull Compactness¶

In [20]:

from vgrid.stats.healpixstats import healpix_compactness_cvh

healpix_compactness_cvh(healpix_inspect)

from vgrid.stats.healpixstats import healpix_compactness_cvh healpix_compactness_cvh(healpix_inspect)

HEALPix Statistics¶

Characteristic Length Scale (CLS - suggested by Ralph Kahn): the diameter of a spherical cap of the same cell's area

In [21]:

from vgrid.stats.healpixstats import healpixstats

healpixstats()

from vgrid.stats.healpixstats import healpixstats healpixstats()

Out[21]:

	resolution	number_of_cells	avg_edge_len_m	avg_cell_area_m2	cls_m
0	0	12	6.519622e+06	4.250547e+13	7.462814e+06
1	1	48	3.259811e+06	1.062637e+13	3.691196e+06
2	2	192	1.629905e+06	2.656592e+12	1.840752e+06
3	3	768	8.149527e+05	6.641479e+11	9.197754e+05
4	4	3072	4.074764e+05	1.660370e+11	4.598128e+05
5	5	12288	2.037382e+05	4.150925e+10	2.298970e+05
6	6	49152	1.018691e+05	1.037731e+10	1.149473e+05
7	7	196608	5.093455e+04	2.594328e+09	5.747353e+04
8	8	786432	2.546727e+04	6.485820e+08	2.873675e+04
9	9	3145728	1.273364e+04	1.621455e+08	1.436837e+04
10	10	12582912	6.366818e+03	4.053637e+07	7.184185e+03
11	11	50331648	3.183409e+03	1.013409e+07	3.592093e+03
12	12	201326592	1.591705e+03	2.533523e+06	1.796046e+03
13	13	805306368	7.958523e+02	6.333808e+05	8.980231e+02
14	14	3221225472	3.979261e+02	1.583452e+05	4.490116e+02
15	15	12884901888	1.989631e+02	3.958630e+04	2.245058e+02
16	16	51539607552	9.948153e+01	9.896576e+03	1.122529e+02
17	17	206158430208	4.974077e+01	2.474144e+03	5.612645e+01
18	18	824633720832	2.487038e+01	6.185360e+02	2.806322e+01
19	19	3298534883328	1.243519e+01	1.546340e+02	1.403161e+01
20	20	13194139533312	6.217596e+00	3.865850e+01	7.015806e+00
21	21	52776558133248	3.108798e+00	9.664625e+00	3.507903e+00
22	22	211106232532992	1.554399e+00	2.416156e+00	1.753951e+00
23	23	844424930131968	7.771995e-01	6.040390e-01	8.769757e-01
24	24	3377699720527872	3.885997e-01	1.510098e-01	4.384879e-01
25	25	13510798882111488	1.942999e-01	3.775244e-02	2.192439e-01
26	26	54043195528445952	9.714994e-02	9.438110e-03	1.096220e-01
27	27	216172782113783808	4.857497e-02	2.359527e-03	5.481098e-02
28	28	864691128455135232	2.428748e-02	5.898819e-04	2.740549e-02
29	29	3458764513820540928	1.214374e-02	1.474705e-04	1.370275e-02