05 isea4t

Working with Open-EAGGR ISEA4T 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

latlon2isea4t¶

In [1]:

from vgrid.conversion.latlon2dggs import latlon2isea4t

lat = 10.775276
lon = 106.706797
res = 12
isea4t_id = latlon2isea4t(lat, lon, res)
isea4t_id

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

Out[1]:

'13102313331320'

ISEA4T to Polygon¶

In [3]:

from vgrid.conversion.dggs2geo.isea4t2geo import isea4t2geo

isea4t_geo = isea4t2geo(isea4t_id)
isea4t_geo

from vgrid.conversion.dggs2geo.isea4t2geo import isea4t2geo isea4t_geo = isea4t2geo(isea4t_id) isea4t_geo

Out[3]:

ISEA4T to GeoJSON¶

In [2]:

from vgrid.conversion.dggs2geo.isea4t2geo import isea4t2geojson

isea4t_geojson = isea4t2geojson(isea4t_id)
# isea4t_geojson

from vgrid.conversion.dggs2geo.isea4t2geo import isea4t2geojson isea4t_geojson = isea4t2geojson(isea4t_id) # isea4t_geojson

Vector to ISEA4T¶

In [2]:

from vgrid.conversion.vector2dggs.vector2isea4t import vector2isea4t

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

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

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

Out[2]:

<Axes: >

ISEA4T Compact¶

In [3]:

from vgrid.conversion.dggscompact.isea4tcompact import isea4tcompact

isea4t_compacted = isea4tcompact(vector_to_isea4t, output_format="gpd")
isea4t_compacted.plot(edgecolor="white")

from vgrid.conversion.dggscompact.isea4tcompact import isea4tcompact isea4t_compacted = isea4tcompact(vector_to_isea4t, output_format="gpd") isea4t_compacted.plot(edgecolor="white")

Out[3]:

<Axes: >

ISEA4T Expand¶

In [4]:

from vgrid.conversion.dggscompact.isea4tcompact import isea4texpand

isea4t_expanded = isea4texpand(vector_to_isea4t, resolution=16, output_format="gpd")
isea4t_expanded.plot(edgecolor="white")

from vgrid.conversion.dggscompact.isea4tcompact import isea4texpand isea4t_expanded = isea4texpand(vector_to_isea4t, resolution=16, output_format="gpd") isea4t_expanded.plot(edgecolor="white")

Out[4]:

<Axes: >

ISEA4T Binning¶

In [8]:

from vgrid.binning.isea4tbin import isea4tbin

file_path = (
    "https://raw.githubusercontent.com/opengeoshub/vopendata/main/csv/dist1_pois.csv"
)
stats = "count"
isea4t_bin = isea4tbin(
    file_path,
    resolution=15,
    stats=stats,
    # numeric_field="confidence",
    # category="category",
    output_format="gpd",
)
isea4t_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.isea4tbin import isea4tbin file_path = ( "https://raw.githubusercontent.com/opengeoshub/vopendata/main/csv/dist1_pois.csv" ) stats = "count" isea4t_bin = isea4tbin( file_path, resolution=15, stats=stats, # numeric_field="confidence", # category="category", output_format="gpd", ) isea4t_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 ISEA4T DGGS: 100%|██████████| 767/767 [00:00<00:00, 7966.79 cells/s]

Out[8]:

<Axes: >

Raster to ISEA4T¶

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 ISEA4T¶

In [10]:

# %pip install folium

# %pip install folium

In [11]:

from vgrid.conversion.raster2dggs.raster2isea4t import raster2isea4t
from vgrid.utils.io import download_file

raster_url = (
    "https://raw.githubusercontent.com/opengeoshub/vopendata/main/raster/rgb.tif"
)
raster_file = download_file(raster_url)

raster_to_isea4t = raster2isea4t(raster_file, output_format="gpd")

# Visualize the output
import folium

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

isea4t_layer = folium.GeoJson(
    raster_to_isea4t,
    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=["isea4t", "band_1", "band_2", "band_3"],
        aliases=["ISEA4T ID", "Band 1", "Band 2", "Band 3"],
        style="""
            background-color: white;
            border: 2px solid black;
            border-radius: 3px;
            box-shadow: 3px;
        """,
    ),
).add_to(m)

m.fit_bounds(isea4t_layer.get_bounds())

# Display the map
m

from vgrid.conversion.raster2dggs.raster2isea4t import raster2isea4t from vgrid.utils.io import download_file raster_url = ( "https://raw.githubusercontent.com/opengeoshub/vopendata/main/raster/rgb.tif" ) raster_file = download_file(raster_url) raster_to_isea4t = raster2isea4t(raster_file, output_format="gpd") # Visualize the output import folium m = folium.Map(tiles="CartoDB positron", max_zoom=28) isea4t_layer = folium.GeoJson( raster_to_isea4t, 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=["isea4t", "band_1", "band_2", "band_3"], aliases=["ISEA4T ID", "Band 1", "Band 2", "Band 3"], style=""" background-color: white; border: 2px solid black; border-radius: 3px; box-shadow: 3px; """, ), ).add_to(m) m.fit_bounds(isea4t_layer.get_bounds()) # Display the map m

rgb.tif already exists. Skip downloading. Set overwrite=True to overwrite.
Cell size: 0.08638527081938627 m2
Nearest ISEA4T resolution determined: 23

Converting raster to ISEA4T: 100%|██████████| 8486/8486 [00:01<00:00, 6673.12 cells/s]

Out[11]:

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

ISEA4T Generator¶

In [12]:

from vgrid.generator.isea4tgrid import isea4tgrid

isea4t_grid = isea4tgrid(
    resolution=2, output_format="gpd", fix_antimeridian="shift_west"
)
# isea4t_grid = isea4tgrid(resolution=16,bbox=[106.699007, 10.762811, 106.717674, 10.778649],output_format="gpd")
isea4t_grid.plot(edgecolor="white")

from vgrid.generator.isea4tgrid import isea4tgrid isea4t_grid = isea4tgrid( resolution=2, output_format="gpd", fix_antimeridian="shift_west" ) # isea4t_grid = isea4tgrid(resolution=16,bbox=[106.699007, 10.762811, 106.717674, 10.778649],output_format="gpd") isea4t_grid.plot(edgecolor="white")

Generating ISEA4T DGGS: 100%|██████████| 320/320 [00:00<00:00, 2721.97 cells/s]

Out[12]:

<Axes: >

ISEA4T Inspect¶

In [13]:

from vgrid.stats.isea4tstats import isea4tinspect

resolution = 5
isea4t_inspect = isea4tinspect(resolution)
isea4t_inspect = isea4t_inspect[isea4t_inspect["crossed"] == False]
isea4t_inspect.head()

from vgrid.stats.isea4tstats import isea4tinspect resolution = 5 isea4t_inspect = isea4tinspect(resolution) isea4t_inspect = isea4t_inspect[isea4t_inspect["crossed"] == False] isea4t_inspect.head()

Generating ISEA4T DGGS: 100%|██████████| 20480/20480 [00:02<00:00, 6998.51 cells/s]

Out[13]:

	isea4t	resolution	center_lat	center_lon	avg_edge_len	cell_area	cell_perimeter	geometry	crossed	norm_area	ipq	zsc	cvh
0	0000000	5	52.798043	-144.000000	249995.286027	2.706746e+10	749985.858081	POLYGON ((-144 51.51071, -145.88123 53.44171, ...	False	1.086805	0.604716	0.777614	1.0
1	0000001	5	53.993894	-144.000000	230922.987391	2.284392e+10	692768.962173	POLYGON ((-144 55.09826, -145.88123 53.44171, ...	False	0.917222	0.598141	0.773379	1.0
2	0000002	5	52.190361	-145.690007	230945.739381	2.284876e+10	692837.218144	POLYGON ((-145.88123 53.44171, -147.18879 51.6...	False	0.917417	0.598150	0.773384	1.0
3	0000003	5	52.190361	-142.309993	230945.739386	2.284876e+10	692837.218159	POLYGON ((-142.11877 53.44171, -144 51.51071, ...	False	0.917417	0.598150	0.773384	1.0
4	0000010	5	50.215849	-144.000000	238726.466768	2.462431e+10	716179.400305	POLYGON ((-144 51.51071, -145.58435 49.56842, ...	False	0.988708	0.603296	0.776702	1.0

ISEA4T Normalized Area Histogram¶

In [14]:

from vgrid.stats.isea4tstats import isea4t_norm_area_hist

isea4t_norm_area_hist(isea4t_inspect)

from vgrid.stats.isea4tstats import isea4t_norm_area_hist isea4t_norm_area_hist(isea4t_inspect)

Distribution of ISEA4T Area Distortions¶

In [15]:

from vgrid.stats.isea4tstats import isea4t_norm_area

isea4t_norm_area(isea4t_inspect)

from vgrid.stats.isea4tstats import isea4t_norm_area isea4t_norm_area(isea4t_inspect)

ISEA4T 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 [16]:

from vgrid.stats.isea4tstats import isea4t_compactness_ipq_hist

isea4t_compactness_ipq_hist(isea4t_inspect)

from vgrid.stats.isea4tstats import isea4t_compactness_ipq_hist isea4t_compactness_ipq_hist(isea4t_inspect)

Distribution of ISEA4T IPQ Compactness¶

In [17]:

from vgrid.stats.isea4tstats import isea4t_compactness_ipq

isea4t_compactness_ipq(isea4t_inspect)

from vgrid.stats.isea4tstats import isea4t_compactness_ipq isea4t_compactness_ipq(isea4t_inspect)

ISEA4T 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 [18]:

from vgrid.stats.isea4tstats import isea4t_compactness_cvh_hist

isea4t_compactness_cvh_hist(isea4t_inspect)

from vgrid.stats.isea4tstats import isea4t_compactness_cvh_hist isea4t_compactness_cvh_hist(isea4t_inspect)

Distribution of ISEA4T Convex hull Compactness¶

In [19]:

from vgrid.stats.isea4tstats import isea4t_compactness_cvh

isea4t_compactness_cvh(isea4t_inspect)

from vgrid.stats.isea4tstats import isea4t_compactness_cvh isea4t_compactness_cvh(isea4t_inspect)

ISEA4T Statistics¶

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

In [1]:

from vgrid.stats.isea4tstats import isea4tstats

isea4tstats()

from vgrid.stats.isea4tstats import isea4tstats isea4tstats()

Out[1]:

	resolution	number_of_cells	avg_edge_len_m	avg_cell_area_m2	cls_m
0	0	20	7.674458e+06	2.550328e+13	5.746990e+06
1	1	80	3.837229e+06	6.375820e+12	2.855171e+06
2	2	320	1.918614e+06	1.593955e+12	1.425344e+06
3	3	1280	9.593072e+05	3.984888e+11	7.123930e+05
4	4	5120	4.796536e+05	9.962219e+10	3.561617e+05
5	5	20480	2.398268e+05	2.490555e+10	1.780765e+05
6	6	81920	1.199134e+05	6.226387e+09	8.903771e+04
7	7	327680	5.995670e+04	1.556597e+09	4.451879e+04
8	8	1310720	2.997835e+04	3.891492e+08	2.225939e+04
9	9	5242880	1.498918e+04	9.728730e+07	1.112969e+04
10	10	20971520	7.494588e+03	2.432182e+07	5.564846e+03
11	11	83886080	3.747294e+03	6.080456e+06	2.782423e+03
12	12	335544320	1.873647e+03	1.520114e+06	1.391211e+03
13	13	1342177280	9.368235e+02	3.800285e+05	6.956057e+02
14	14	5368709120	4.684117e+02	9.500713e+04	3.478029e+02
15	15	21474836480	2.342059e+02	2.375178e+04	1.739014e+02
16	16	85899345920	1.171029e+02	5.937945e+03	8.695071e+01
17	17	343597383680	5.855147e+01	1.484486e+03	4.347536e+01
18	18	1374389534720	2.927573e+01	3.711216e+02	2.173768e+01
19	19	5497558138880	1.463787e+01	9.278040e+01	1.086884e+01
20	20	21990232555520	7.318933e+00	2.319510e+01	5.434420e+00
21	21	87960930222080	3.659467e+00	5.798775e+00	2.717210e+00
22	22	351843720888320	1.829733e+00	1.449694e+00	1.358605e+00
23	23	1407374883553280	9.148667e-01	3.624234e-01	6.793025e-01
24	24	5629499534213120	4.574333e-01	9.060586e-02	3.396512e-01
25	25	22517998136852480	2.287167e-01	2.265146e-02	1.698256e-01
26	26	90071992547409920	1.143583e-01	5.662866e-03	8.491281e-02
27	27	360287970189639680	5.717917e-02	1.415716e-03	4.245640e-02
28	28	1441151880758558720	2.858958e-02	3.539291e-04	2.122820e-02
29	29	5764607523034234880	1.429479e-02	8.848228e-05	1.061410e-02
30	30	23058430092136939520	7.147396e-03	2.212057e-05	5.307050e-03
31	31	92233720368547758080	3.573698e-03	5.530143e-06	2.653525e-03
32	32	368934881474191032320	1.786849e-03	1.382536e-06	1.326763e-03
33	33	1475739525896764129280	8.934245e-04	3.456339e-07	6.633813e-04
34	34	5902958103587056517120	4.467122e-04	8.640848e-08	3.316907e-04
35	35	23611832414348226068480	2.233561e-04	2.160212e-08	1.658453e-04
36	36	94447329657392904273920	1.116781e-04	5.400530e-09	8.292266e-05
37	37	377789318629571617095680	5.583903e-05	1.350132e-09	4.146133e-05
38	38	1511157274518286468382720	2.791952e-05	3.375331e-10	2.073067e-05
39	39	6044629098073145873530880	1.395976e-05	8.438328e-11	1.036533e-05