02 s2

Working with S2 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¶

In [1]:

# %pip install vgrid --upgrade

# %pip install vgrid --upgrade

latlon2s2¶

In [2]:

from vgrid.conversion.latlon2dggs import latlon2s2

lat = 10.775276
lon = 106.706797
res = 12
s2_token = latlon2s2(lat, lon, res)
s2_token

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

Out[2]:

'31752f5'

S2 to Polygon¶

In [3]:

from vgrid.conversion.dggs2geo.s22geo import s22geo

s2_geo = s22geo(s2_token)
s2_geo

from vgrid.conversion.dggs2geo.s22geo import s22geo s2_geo = s22geo(s2_token) s2_geo

Out[3]:

S2 to GeoJSON¶

In [3]:

from vgrid.conversion.dggs2geo.s22geo import s22geojson

s2_geojson = s22geojson(s2_token)
# s2_geojson

from vgrid.conversion.dggs2geo.s22geo import s22geojson s2_geojson = s22geojson(s2_token) # s2_geojson

Vector to S2¶

In [5]:

from vgrid.conversion.vector2dggs.vector2s2 import vector2s2

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

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

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

Out[5]:

<Axes: >

S2 Compact¶

In [6]:

from vgrid.conversion.dggscompact.s2compact import s2compact

s2_compacted = s2compact(vector_to_s2, s2_token="s2", output_format="gpd")
s2_compacted.plot(edgecolor="white")

from vgrid.conversion.dggscompact.s2compact import s2compact s2_compacted = s2compact(vector_to_s2, s2_token="s2", output_format="gpd") s2_compacted.plot(edgecolor="white")

Out[6]:

<Axes: >

S2 Expand¶

In [7]:

from vgrid.conversion.dggscompact.s2compact import s2expand

s2_expanded = s2expand(vector_to_s2, resolution=18, output_format="gpd")
s2_expanded.plot(edgecolor="white")

from vgrid.conversion.dggscompact.s2compact import s2expand s2_expanded = s2expand(vector_to_s2, resolution=18, output_format="gpd") s2_expanded.plot(edgecolor="white")

Out[7]:

<Axes: >

S2 Binning¶

In [8]:

from vgrid.binning.s2bin import s2bin

file_path = (
    "https://raw.githubusercontent.com/opengeoshub/vopendata/main/csv/dist1_pois.csv"
)
stats = "count"
s2_bin = s2bin(
    file_path,
    resolution=16,
    stats=stats,
    # numeric_field="confidence",
    # category="category",
    output_format="gpd",
)
s2_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.s2bin import s2bin file_path = ( "https://raw.githubusercontent.com/opengeoshub/vopendata/main/csv/dist1_pois.csv" ) stats = "count" s2_bin = s2bin( file_path, resolution=16, stats=stats, # numeric_field="confidence", # category="category", output_format="gpd", ) s2_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 DGGS: 100%|██████████| 737/737 [00:00<00:00, 2173.65 cells/s]

Out[8]:

<Axes: >

Raster to S2¶

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

In [10]:

# %pip install folium

# %pip install folium

In [11]:

from vgrid.conversion.raster2dggs.raster2s2 import raster2s2


raster_to_s2 = raster2s2(raster_file, output_format="gpd")
# Visualize raster_to_s2
import folium

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

s2_layer = folium.GeoJson(
    raster_to_s2,
    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=["s2", "resolution", "band_1", "band_2", "band_3", "cell_area"],
        aliases=["S2 Token", "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(s2_layer.get_bounds())

# Display the map
m

from vgrid.conversion.raster2dggs.raster2s2 import raster2s2 raster_to_s2 = raster2s2(raster_file, output_format="gpd") # Visualize raster_to_s2 import folium m = folium.Map(tiles="CartoDB positron", max_zoom=28) s2_layer = folium.GeoJson( raster_to_s2, 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=["s2", "resolution", "band_1", "band_2", "band_3", "cell_area"], aliases=["S2 Token", "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(s2_layer.get_bounds()) # Display the map m

Cell size: 0.08638527081938627 m2
Nearest S2 resolution determined: 24

Converting raster to S2: 100%|██████████| 8231/8231 [00:06<00:00, 1250.98 cells/s]

Out[11]:

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

S2 Generator¶

In [12]:

from vgrid.generator.s2grid import s2grid

s2_grid = s2grid(resolution=2, output_format="gpd", fix_antimeridian="shift_east")
# s2_grid = s2grid(resolution=16,bbox=[106.699007, 10.762811, 106.717674, 10.778649],output_format="gpd")
# s2_grid.to_crs('proj=moll').plot(edgecolor="white")
# s2_grid.to_file("s2_2_shift_east.geojson")
s2_grid.plot(edgecolor="white")

from vgrid.generator.s2grid import s2grid s2_grid = s2grid(resolution=2, output_format="gpd", fix_antimeridian="shift_east") # s2_grid = s2grid(resolution=16,bbox=[106.699007, 10.762811, 106.717674, 10.778649],output_format="gpd") # s2_grid.to_crs('proj=moll').plot(edgecolor="white") # s2_grid.to_file("s2_2_shift_east.geojson") s2_grid.plot(edgecolor="white")

Generating DGGS: 100%|██████████| 96/96 [00:00<00:00, 355.71 cells/s]

Out[12]:

<Axes: >

S2 Inspect¶

In [13]:

from vgrid.stats.s2stats import s2inspect

resolution = 6
s2_inspect = s2inspect(resolution)
s2_inspect = s2_inspect[
    ~s2_inspect["crossed"]
]  # remove cells that cross the Antimeridian

from vgrid.stats.s2stats import s2inspect resolution = 6 s2_inspect = s2inspect(resolution) s2_inspect = s2_inspect[ ~s2_inspect["crossed"] ] # remove cells that cross the Antimeridian

Generating DGGS: 100%|██████████| 24576/24576 [00:21<00:00, 1143.54 cells/s]

In [14]:

s2_inspect_sorted = s2_inspect.sort_values("cvh", ascending=False, na_position="last")
s2_inspect_sorted.head()

s2_inspect_sorted = s2_inspect.sort_values("cvh", ascending=False, na_position="last") s2_inspect_sorted.head()

Out[14]:

	s2	resolution	center_lat	center_lon	avg_edge_len	cell_area	cell_perimeter	geometry	crossed	norm_area	ipq	zsc	cvh
0	0001	6	-34.976691	-44.394709	126254.074677	1.391628e+10	505016.298708	POLYGON ((-45 -35.26439, -43.79085 -35.82442, ...	False	0.670514	0.685681	0.828047	1.0
1	0003	6	-35.529669	-43.170222	127426.101397	1.427186e+10	509704.405587	POLYGON ((-43.79085 -35.82442, -42.55096 -36.3...	False	0.687647	0.690325	0.830846	1.0
2	0005	6	-34.378795	-43.170143	128679.032199	1.464789e+10	514716.128798	POLYGON ((-43.79085 -34.68431, -42.55096 -35.2...	False	0.705765	0.694783	0.833525	1.0
3	0007	6	-33.834099	-44.394632	127398.731047	1.426656e+10	509594.924188	POLYGON ((-45 -34.13233, -43.79085 -34.68431, ...	False	0.687392	0.690365	0.830871	1.0
4	0009	6	-32.679004	-44.394555	128450.519201	1.460280e+10	513802.076805	POLYGON ((-45 -32.98768, -43.79085 -33.53064, ...	False	0.703593	0.695111	0.833721	1.0

S2 Normalized Area Histogram¶

In [15]:

from vgrid.stats.s2stats import s2_norm_area_hist

s2_norm_area_hist(s2_inspect)

from vgrid.stats.s2stats import s2_norm_area_hist s2_norm_area_hist(s2_inspect)

S2 IPQ Compactness Histogram¶

Distribution of S2 Area Distortions¶

In [16]:

from vgrid.stats.s2stats import s2_norm_area

s2_norm_area(s2_inspect)

from vgrid.stats.s2stats import s2_norm_area s2_norm_area(s2_inspect)

S2 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.s2stats import s2_compactness_ipq_hist

s2_compactness_ipq_hist(s2_inspect)

from vgrid.stats.s2stats import s2_compactness_ipq_hist s2_compactness_ipq_hist(s2_inspect)

Distribution of S2 IPQ Compactness¶

In [18]:

from vgrid.stats.s2stats import s2_compactness_ipq

s2_compactness_ipq(s2_inspect)

from vgrid.stats.s2stats import s2_compactness_ipq s2_compactness_ipq(s2_inspect)

S2 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.s2stats import s2_compactness_cvh_hist

s2_compactness_cvh_hist(s2_inspect)

from vgrid.stats.s2stats import s2_compactness_cvh_hist s2_compactness_cvh_hist(s2_inspect)

Distribution of S2 Convex hull Compactness¶

In [20]:

from vgrid.stats.s2stats import s2_compactness_cvh

s2_compactness_cvh(s2_inspect)

from vgrid.stats.s2stats import s2_compactness_cvh s2_compactness_cvh(s2_inspect)

S2 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.s2stats import s2stats

s2stats()

from vgrid.stats.s2stats import s2stats s2stats()

Out[21]:

	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	24	4.610069e+06	2.125273e+13	5.238725e+06
2	2	96	2.305034e+06	5.313184e+12	2.605490e+06
3	3	384	1.152517e+06	1.328296e+12	1.301042e+06
4	4	1536	5.762586e+05	3.320740e+11	6.503088e+05
5	5	6144	2.881293e+05	8.301849e+10	3.251279e+05
6	6	24576	1.440646e+05	2.075462e+10	1.625607e+05
7	7	98304	7.203232e+04	5.188656e+09	8.127991e+04
8	8	393216	3.601616e+04	1.297164e+09	4.063990e+04
9	9	1572864	1.800808e+04	3.242910e+08	2.031995e+04
10	10	6291456	9.004041e+03	8.107275e+07	1.015997e+04
11	11	25165824	4.502020e+03	2.026819e+07	5.079986e+03
12	12	100663296	2.251010e+03	5.067047e+06	2.539993e+03
13	13	402653184	1.125505e+03	1.266762e+06	1.269996e+03
14	14	1610612736	5.627525e+02	3.166904e+05	6.349982e+02
15	15	6442450944	2.813763e+02	7.917260e+04	3.174991e+02
16	16	25769803776	1.406881e+02	1.979315e+04	1.587496e+02
17	17	103079215104	7.034407e+01	4.948288e+03	7.937478e+01
18	18	412316860416	3.517203e+01	1.237072e+03	3.968739e+01
19	19	1649267441664	1.758602e+01	3.092680e+02	1.984369e+01
20	20	6597069766656	8.793008e+00	7.731700e+01	9.921847e+00
21	21	26388279066624	4.396504e+00	1.932925e+01	4.960924e+00
22	22	105553116266496	2.198252e+00	4.832312e+00	2.480462e+00
23	23	422212465065984	1.099126e+00	1.208078e+00	1.240231e+00
24	24	1688849860263936	5.495630e-01	3.020195e-01	6.201155e-01
25	25	6755399441055744	2.747815e-01	7.550488e-02	3.100577e-01
26	26	27021597764222976	1.373908e-01	1.887622e-02	1.550289e-01
27	27	108086391056891904	6.869538e-02	4.719055e-03	7.751443e-02
28	28	432345564227567616	3.434769e-02	1.179764e-03	3.875722e-02
29	29	1729382256910270464	1.717384e-02	2.949409e-04	1.937861e-02
30	30	6917529027641081856	8.586922e-03	7.373523e-05	9.689304e-03