04 rhealpix

Working with rHEALpix (Manaaki Whenua – Landcare Research) 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 [1]:

# %pip install vgrid --upgrade

# %pip install vgrid --upgrade

latlon2rhealpix¶

In [2]:

from vgrid.conversion.latlon2dggs import latlon2rhealpix

lat = 10.775276
lon = 106.706797
res = 12
rhealpix_id = latlon2rhealpix(lat, lon, res)
rhealpix_id

from vgrid.conversion.latlon2dggs import latlon2rhealpix lat = 10.775276 lon = 106.706797 res = 12 rhealpix_id = latlon2rhealpix(lat, lon, res) rhealpix_id

Out[2]:

'R312603625535'

rHEALPix to Polygon¶

In [3]:

from vgrid.conversion.dggs2geo.rhealpix2geo import rhealpix2geo

rhealpix_geo = rhealpix2geo(rhealpix_id)
rhealpix_geo

from vgrid.conversion.dggs2geo.rhealpix2geo import rhealpix2geo rhealpix_geo = rhealpix2geo(rhealpix_id) rhealpix_geo

Out[3]:

rHEALPix to GeoJSON¶

In [2]:

from vgrid.conversion.dggs2geo.rhealpix2geo import rhealpix2geojson

rhealpix_geojson = rhealpix2geojson(rhealpix_id)
# rhealpix_geojson

from vgrid.conversion.dggs2geo.rhealpix2geo import rhealpix2geojson rhealpix_geojson = rhealpix2geojson(rhealpix_id) # rhealpix_geojson

Vector to rHEALPix¶

In [5]:

from vgrid.conversion.vector2dggs.vector2rhealpix import vector2rhealpix

file_path = "https://raw.githubusercontent.com/opengeoshub/vopendata/main/shape/polygon2.geojson"
vector_to_rhealpix = vector2rhealpix(
    file_path,
    resolution=10,
    compact=True,
    topology=False,
    predicate="intersects",
    output_format="gpd",
)
# Visualize the output
vector_to_rhealpix.plot(edgecolor="white")

from vgrid.conversion.vector2dggs.vector2rhealpix import vector2rhealpix file_path = "https://raw.githubusercontent.com/opengeoshub/vopendata/main/shape/polygon2.geojson" vector_to_rhealpix = vector2rhealpix( file_path, resolution=10, compact=True, topology=False, predicate="intersects", output_format="gpd", ) # Visualize the output vector_to_rhealpix.plot(edgecolor="white")

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

Out[5]:

<Axes: >

rHEALPix Compact¶

In [6]:

from vgrid.conversion.dggscompact.rhealpixcompact import rhealpixcompact

rhealpix_compacted = rhealpixcompact(
    vector_to_rhealpix, rhealpix_id="rhealpix", output_format="gpd"
)
rhealpix_compacted.plot(edgecolor="white")

from vgrid.conversion.dggscompact.rhealpixcompact import rhealpixcompact rhealpix_compacted = rhealpixcompact( vector_to_rhealpix, rhealpix_id="rhealpix", output_format="gpd" ) rhealpix_compacted.plot(edgecolor="white")

Out[6]:

<Axes: >

rHEALPix Expand¶

In [7]:

from vgrid.conversion.dggscompact.rhealpixcompact import rhealpixexpand

rhealpix_expanded = rhealpixexpand(
    vector_to_rhealpix, resolution=11, output_format="gpd"
)
rhealpix_expanded.plot(edgecolor="white")

from vgrid.conversion.dggscompact.rhealpixcompact import rhealpixexpand rhealpix_expanded = rhealpixexpand( vector_to_rhealpix, resolution=11, output_format="gpd" ) rhealpix_expanded.plot(edgecolor="white")

Out[7]:

<Axes: >

rHEALPix Binning¶

In [8]:

from vgrid.binning.rhealpixbin import rhealpixbin

file_path = (
    "https://raw.githubusercontent.com/opengeoshub/vopendata/main/csv/dist1_pois.csv"
)
stats = "count"
rhealpix_bin = rhealpixbin(
    file_path,
    resolution=10,
    stats=stats,
    # numeric_field="confidence",
    # category="category",
    output_format="gpd",
)
rhealpix_bin.plot(
    column=stats,  # numeric column to base the colors on
    cmap="Spectral_r",  # color scheme (matplotlib colormap)
    legend=True,
    linewidth=0.2,  # boundary width (optional)
)

from vgrid.binning.rhealpixbin import rhealpixbin file_path = ( "https://raw.githubusercontent.com/opengeoshub/vopendata/main/csv/dist1_pois.csv" ) stats = "count" rhealpix_bin = rhealpixbin( file_path, resolution=10, stats=stats, # numeric_field="confidence", # category="category", output_format="gpd", ) rhealpix_bin.plot( column=stats, # numeric column to base the colors on cmap="Spectral_r", # color scheme (matplotlib colormap) legend=True, linewidth=0.2, # boundary width (optional) )

Generating rHEALPix DGGS: 100%|██████████| 726/726 [00:00<00:00, 7611.88 cells/s]

Out[8]:

<Axes: >

Raster to rHEALPix¶

Download and open raster¶

In [9]:

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)

WARNING [rasterio._env:368  open()] CPLE_AppDefined in PROJ: proj_create_from_database: Cannot find proj.db

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["World Geodetic System 1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST]]'), 'transform': Affine(2.6640125000199077e-06, 0.0, 106.708118755,
       0.0, -2.6640136054383103e-06, 10.812568272)}

Out[9]:

<Axes: >

Convert raster to rHEALPix¶

In [10]:

# %pip install folium

# %pip install folium

In [11]:

from vgrid.conversion.raster2dggs.raster2rhealpix import raster2rhealpix


raster_to_rhealpix = raster2rhealpix(raster_file, output_format="gpd")
# Visualize the output
import folium

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

rhealpix_layer = folium.GeoJson(
    raster_to_rhealpix,
    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=["rhealpix", "resolution", "band_1", "band_2", "band_3", "cell_area"],
        aliases=[
            "rHEALPix 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(rhealpix_layer.get_bounds())

# Display the map
m

from vgrid.conversion.raster2dggs.raster2rhealpix import raster2rhealpix raster_to_rhealpix = raster2rhealpix(raster_file, output_format="gpd") # Visualize the output import folium m = folium.Map(tiles="CartoDB positron", max_zoom=28) rhealpix_layer = folium.GeoJson( raster_to_rhealpix, 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=["rhealpix", "resolution", "band_1", "band_2", "band_3", "cell_area"], aliases=[ "rHEALPix 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(rhealpix_layer.get_bounds()) # Display the map m

Cell size: 0.08638527081938627 m2
Nearest rHEALPix resolution determined: 15

Converting raster to rHEALPix: 100%|██████████| 7416/7416 [00:03<00:00, 2365.58 cells/s]

Out[11]:

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

rHEALPix Generator¶

In [24]:

from vgrid.generator.rhealpixgrid import rhealpixgrid

rhealpix_grid = rhealpixgrid(
    resolution=1, output_format="gpd", fix_antimeridian="shift_west"
)
# rhealpix_grid = rhealpixgrid(resolution=10,bbox=[106.699007, 10.762811, 106.717674, 10.778649],output_format="gpd")
# rhealpix_grid.to_file("rhealpix_1_origin.geojson")
rhealpix_grid.plot(edgecolor="white")

from vgrid.generator.rhealpixgrid import rhealpixgrid rhealpix_grid = rhealpixgrid( resolution=1, output_format="gpd", fix_antimeridian="shift_west" ) # rhealpix_grid = rhealpixgrid(resolution=10,bbox=[106.699007, 10.762811, 106.717674, 10.778649],output_format="gpd") # rhealpix_grid.to_file("rhealpix_1_origin.geojson") rhealpix_grid.plot(edgecolor="white")

Generating rHEALPix DGGS: 100%|██████████| 54/54 [00:00<00:00, 1724.09 cells/s]

Out[24]:

<Axes: >

rHEALPix Inspect¶

In [14]:

from vgrid.stats.rhealpixstats import rhealpixinspect

resolution = 4
rhealpix_inspect = rhealpixinspect(resolution)
rhealpix_inspect = rhealpix_inspect[rhealpix_inspect["crossed"] == False]
rhealpix_inspect.head()

from vgrid.stats.rhealpixstats import rhealpixinspect resolution = 4 rhealpix_inspect = rhealpixinspect(resolution) rhealpix_inspect = rhealpix_inspect[rhealpix_inspect["crossed"] == False] rhealpix_inspect.head()

Generating rHEALPix DGGS: 100%|██████████| 39366/39366 [00:20<00:00, 1946.44 cells/s]

Out[14]:

	rhealpix	resolution	center_lat	center_lon	avg_edge_len	cell_area	cell_perimeter	geometry	crossed	norm_area	ipq	zsc	cvh
2	N0002	4	42.571866	177.758908	128917.630146	1.295731e+10	515670.520583	POLYGON ((177.72152 43.20062, 178.86076 43.200...	False	1.000024	0.612323	0.782501	1.0
3	N0003	4	42.571866	-178.884142	129535.796654	1.295728e+10	518143.186616	POLYGON ((-180 43.20062, -178.86076 43.20062, ...	False	1.000021	0.606491	0.778765	1.0
6	N0006	4	42.571866	-177.758908	128917.630146	1.295731e+10	515670.520583	POLYGON ((-178.86076 43.20062, -177.72152 43.2...	False	1.000024	0.612323	0.782501	1.0
7	N0007	4	43.831553	-178.855764	129576.559309	1.295724e+10	518306.237234	POLYGON ((-180 44.45712, -178.83117 44.45712, ...	False	1.000018	0.606107	0.778519	1.0
9	N0010	4	42.571866	176.633673	128311.046607	1.295735e+10	513244.186428	POLYGON ((176.58228 43.20062, 177.72152 43.200...	False	1.000026	0.618127	0.786201	1.0

rHEALPix Normalized Area Histogram¶

In [15]:

from vgrid.stats.rhealpixstats import rhealpix_norm_area_hist

rhealpix_norm_area_hist(rhealpix_inspect)

from vgrid.stats.rhealpixstats import rhealpix_norm_area_hist rhealpix_norm_area_hist(rhealpix_inspect)

Distribution of rHEALPix Area Distortions¶

In [16]:

from vgrid.stats.rhealpixstats import rhealpix_norm_area

rhealpix_norm_area(rhealpix_inspect)

from vgrid.stats.rhealpixstats import rhealpix_norm_area rhealpix_norm_area(rhealpix_inspect)

rHEALPix 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.rhealpixstats import rhealpix_compactness_ipq_hist

rhealpix_compactness_ipq_hist(rhealpix_inspect)

from vgrid.stats.rhealpixstats import rhealpix_compactness_ipq_hist rhealpix_compactness_ipq_hist(rhealpix_inspect)

Distribution of rHEALPix IPQ Compactness¶

In [18]:

from vgrid.stats.rhealpixstats import rhealpix_compactness_ipq

rhealpix_compactness_ipq(rhealpix_inspect)

from vgrid.stats.rhealpixstats import rhealpix_compactness_ipq rhealpix_compactness_ipq(rhealpix_inspect)

rHEALPix 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.rhealpixstats import rhealpix_compactness_cvh_hist

rhealpix_compactness_cvh_hist(rhealpix_inspect)

from vgrid.stats.rhealpixstats import rhealpix_compactness_cvh_hist rhealpix_compactness_cvh_hist(rhealpix_inspect)

Distribution of rHEALPix Convex hull Compactness¶

In [20]:

from vgrid.stats.rhealpixstats import rhealpix_compactness_cvh

rhealpix_compactness_cvh(rhealpix_inspect)

from vgrid.stats.rhealpixstats import rhealpix_compactness_cvh rhealpix_compactness_cvh(rhealpix_inspect)

rHEALPix Statistics¶

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

In [3]:

from vgrid.stats.rhealpixstats import rhealpixstats

rhealpixstats()

from vgrid.stats.rhealpixstats import rhealpixstats rhealpixstats()

Out[3]:

	resolution	number_of_cells	avg_edge_len_m	avg_cell_area_m2	cls_m
0	0	6	9.220138e+06	8.501094e+13	1.071691e+07
1	1	54	3.073379e+06	9.445660e+12	3.478731e+06
2	2	486	1.024460e+06	1.049518e+12	1.156376e+06
3	3	4374	3.414866e+05	1.166131e+11	3.853410e+05
4	4	39366	1.138289e+05	1.295701e+10	1.284427e+05
5	5	354294	3.794295e+04	1.439668e+09	4.281406e+04
6	6	3188646	1.264765e+04	1.599631e+08	1.427135e+04
7	7	28697814	4.215884e+03	1.777368e+07	4.757115e+03
8	8	258280326	1.405295e+03	1.974853e+06	1.585705e+03
9	9	2324522934	4.684315e+02	2.194281e+05	5.285684e+02
10	10	20920706406	1.561438e+02	2.438090e+04	1.761895e+02
11	11	188286357654	5.204795e+01	2.708989e+03	5.872982e+01
12	12	1694577218886	1.734932e+01	3.009987e+02	1.957661e+01
13	13	15251194969974	5.783105e+00	3.344431e+01	6.525535e+00
14	14	137260754729766	1.927702e+00	3.716034e+00	2.175178e+00
15	15	1235346792567894	6.425672e-01	4.128927e-01	7.250595e-01