Tiling system#
The core class in the Equi7Grid is the Equi7TilingSystem and represents an Equi7Grid continental subgrid as a regular, multi-level (e.g., zoom, tiling level, …) projected tiling system. A projected tiling system is defined as follows:
A projection represented via a CRS definition (EPSG, PROJ4, WKT, …) and a projection zone defining the validity of coordinates.
A regular tiling scheme slicing the projection extent into smaller regular tiles with a given tile and pixel size (= sampling); tiling schemes are stacked in a tiling system to form a hierarchy of tiling levels (generic realisation of OGC’s Tile Matrix Set). Equi7Grid’s standard tiling system consists of three tiling levels: T1 (100km tile size), T3 (300km tile size), and T6 (600km tile size).
On a code-level, Equi7TilingSystem inherits from pytileproj’s RegularProjTilingSystem and extends it with a definition of land masses for each continent to decide if a tile is covering land or not.
Lets dive deeper into how an Equi7TilingSystem instance can be used. As an example, we assume we retrieved the instance for the European Equi7Grid tiling system (“EU”).
e7eu = ...
print(e7eu) # noqa: T201
Equi7TilingSystem
-----------------
Name:
EU
Projection:
+proj=aeqd +lat_0=53 +lon_0=24 +x_0=5837287.82 +y_0=2121415.696 +datum=WGS84 +units=m +no_defs +type=crs
Tilings:
{0: RegularTiling(T6, 500.0), 1: RegularTiling(T3, 500.0)}
Projection system interactions#
Often it is required to transform coordinates back and forth between different projection systems. Equi7TilingSystem provides several methods to do this in a straighforward way. If you want transform geographic coordinates to projected coordinates, you can use lonlat_to_xy:
lon, lat = 16.3926, 48.1674
proj_coord = e7eu.lonlat_to_xy(lon, lat)
proj_coord
ProjCoord(5271748.127753316, 1613228.4676009014)
The same can be done the other way around:
e7eu.xy_to_lonlat(proj_coord.x, proj_coord.y)
GeogCoord(16.392599999999998, 48.1674)
Note that for coordinates outside the (continental) projection zone, an error will be raised:
lon, lat = -100, 30
e7eu.lonlat_to_xy(lon, lat)
---------------------------------------------------------------------------
GeomOutOfZoneError Traceback (most recent call last)
Cell In[6], line 2
1 lon, lat = -100, 30
----> 2 e7eu.lonlat_to_xy(lon, lat)
File ~/work/Equi7Grid/Equi7Grid/.venv/lib/python3.12/site-packages/pytileproj/tiling_system.py:206, in ProjSystem.lonlat_to_xy(self, lon, lat)
204 coord = ProjCoord(x, y, self._crs)
205 else:
--> 206 raise GeomOutOfZoneError(shapely.Point((lon, lat)))
208 return coord
GeomOutOfZoneError: The given Point ('POINT (-100 30)') is not within the zone boundaries of the projection.
You can check in advance if a point or geometry is within the projection by using Python’s in operator:
from pytileproj import GeogCoord
coord = GeogCoord(lon, lat)
coord in e7eu # ty: ignore[unsupported-operator]
False
As you have seen in this example pytileproj offers several classes to store geometries along with projection information: ProjCoord and ProjGeom as the core classes and GeogCoord and GeogGeom to represent geographic geometries (more details can be found here).
If you want to know the units of the projection system, you can retrieve this information via:
e7eu.unit
'metre'
Tiles#
Equi7 tiling systems create tile objects on the fly, since their regular tiling scheme allows an efficient computation of tile properties. A tile can be created by providing a location of interest with the following methods: get_tile_from_lonlat, get_tile_from_xy, and get_tile_from_coord. Here is an example:
lon, lat = 16.3926, 48.1674
e7eu.get_tile_from_lonlat(lon, lat, tiling_id="T6")
EU_E048N012T6
If you are interested in retrieving a list of tiles intersecting with a region of interest in geographic coordinates, you can use get_tiles_in_geog_bbox. It returns a generator to enable lazy tile retrieval.
list(e7eu.get_tiles_in_geog_bbox((16, 48, 18, 50), tiling_id="T6"))
[EU_E048N012T6, EU_E048N018T6, EU_E054N012T6, EU_E054N018T6]
Important
For a geographical bounding box the order of longitude coordinates is important. As an example: (-170, 20, 150, 30) is a large bounding box going from -170° west to 150° east and (150, 20, -170, 30) is a small bounding box crossing the antimeridian.
Similar queries can be done with projected/native coordinates
list(e7eu.get_tiles_in_bbox((5155982, 1505898, 6323863, 2140396), tiling_id="T6"))
[EU_E048N012T6,
EU_E048N018T6,
EU_E054N012T6,
EU_E054N018T6,
EU_E060N012T6,
EU_E060N018T6]
or by using a geometry object.
from pytileproj import GeogGeom
from shapely import Polygon
geom = Polygon([(16, 48), (18, 48), (18, 50), (16, 50)])
ggeom = GeogGeom(geom=geom)
list(e7eu.get_tiles_in_geom(ggeom, tiling_id="T3"))
[EU_E051N015T3, EU_E051N018T3, EU_E054N015T3, EU_E054N018T3]
Since tiling systems represent an hierarchy of different tilings, Equi7TilingSystem also allows to trace parent and children tiles. You can obtain the parent of a tile with
e7eu.get_parent_from_name("EU_E048N012T3")
EU_E048N012T6
and the children with
list(e7eu.get_children_from_name("EU_E054N012T6"))
[EU_E054N012T3, EU_E054N015T3, EU_E057N012T3, EU_E057N015T3]
Export#
Note
The export methods requires additional dependencies, which are not installed by default. You can install them via uv pip install equi7grid[export].
A tiling system offers several export methods, e.g., an export to a GeoDataFrame (to_geodataframe) or a shapefile (to_shapefile). Here is an example how to generate a GeoDataFrame:
e7eu.to_geodataframe()
| crs | n_rows | n_cols | geotrans | px_origin | name | covers_land | tiling_name | tiling_level | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | EPSG:27704 | 1200 | 1200 | (0.0, 500.0, 0.0, 3000000.0, 0.0, -500.0) | ll | EU_E000N024T6 | False | T6 | 0 | POLYGON ((0 2400000, 600000 2400000, 600000 30... |
| 1 | EPSG:27704 | 1200 | 1200 | (0.0, 500.0, 0.0, 2400000.0, 0.0, -500.0) | ll | EU_E000N018T6 | False | T6 | 0 | POLYGON ((0 1800000, 600000 1800000, 600000 24... |
| 2 | EPSG:27704 | 1200 | 1200 | (0.0, 500.0, 0.0, 1800000.0, 0.0, -500.0) | ll | EU_E000N012T6 | False | T6 | 0 | POLYGON ((0 1200000, 600000 1200000, 600000 18... |
| 3 | EPSG:27704 | 1200 | 1200 | (600000.0, 500.0, 0.0, 3600000.0, 0.0, -500.0) | ll | EU_E006N030T6 | False | T6 | 0 | POLYGON ((600000 3000000, 1200000 3000000, 120... |
| 4 | EPSG:27704 | 1200 | 1200 | (600000.0, 500.0, 0.0, 3000000.0, 0.0, -500.0) | ll | EU_E006N024T6 | False | T6 | 0 | POLYGON ((600000 2400000, 1200000 2400000, 120... |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 434 | EPSG:27704 | 600 | 600 | (7800000.0, 500.0, 0.0, 1200000.0, 0.0, -500.0) | ll | EU_E078N009T3 | True | T3 | 1 | POLYGON ((7800000 900000, 8100000 900000, 8100... |
| 435 | EPSG:27704 | 600 | 600 | (7800000.0, 500.0, 0.0, 900000.0, 0.0, -500.0) | ll | EU_E078N006T3 | True | T3 | 1 | POLYGON ((7800000 600000, 8100000 600000, 8100... |
| 436 | EPSG:27704 | 600 | 600 | (8100000.0, 500.0, 0.0, 1500000.0, 0.0, -500.0) | ll | EU_E081N012T3 | False | T3 | 1 | POLYGON ((8100000 1200000, 8400000 1200000, 84... |
| 437 | EPSG:27704 | 600 | 600 | (8100000.0, 500.0, 0.0, 1200000.0, 0.0, -500.0) | ll | EU_E081N009T3 | False | T3 | 1 | POLYGON ((8100000 900000, 8400000 900000, 8400... |
| 438 | EPSG:27704 | 600 | 600 | (8100000.0, 500.0, 0.0, 900000.0, 0.0, -500.0) | ll | EU_E081N006T3 | False | T3 | 1 | POLYGON ((8100000 600000, 8400000 600000, 8400... |
439 rows × 10 columns
If you want to share an OGC compliant definition of the tiling system, you can use to_ogc_standard or to_ogc_json:
import json
import pprint
from pathlib import Path
ogc_ts_path = Path("my_equi7_system.json")
e7eu.to_ogc_json(ogc_ts_path)
with ogc_ts_path.open() as f:
ogc_ts = json.load(f)
pprint.pprint(ogc_ts) # noqa: T203
{'boundingBox': None,
'crs': {'wkt': {'$schema': 'https://proj.org/schemas/v0.7/projjson.schema.json',
'area': 'Europe including Russia west of the Ural Mountains.',
'base_crs': {'coordinate_system': {'axis': [{'abbreviation': 'Lat',
'direction': 'north',
'name': 'Geodetic '
'latitude',
'unit': 'degree'},
{'abbreviation': 'Lon',
'direction': 'east',
'name': 'Geodetic '
'longitude',
'unit': 'degree'}],
'subtype': 'ellipsoidal'},
'datum_ensemble': {'accuracy': '2.0',
'ellipsoid': {'inverse_flattening': 298.257223563,
'name': 'WGS 84',
'semi_major_axis': 6378137},
'id': {'authority': 'EPSG',
'code': 6326},
'members': [{'id': {'authority': 'EPSG',
'code': 1166},
'name': 'World '
'Geodetic '
'System '
'1984 '
'(Transit)'},
{'id': {'authority': 'EPSG',
'code': 1152},
'name': 'World '
'Geodetic '
'System '
'1984 '
'(G730)'},
{'id': {'authority': 'EPSG',
'code': 1153},
'name': 'World '
'Geodetic '
'System '
'1984 '
'(G873)'},
{'id': {'authority': 'EPSG',
'code': 1154},
'name': 'World '
'Geodetic '
'System '
'1984 '
'(G1150)'},
{'id': {'authority': 'EPSG',
'code': 1155},
'name': 'World '
'Geodetic '
'System '
'1984 '
'(G1674)'},
{'id': {'authority': 'EPSG',
'code': 1156},
'name': 'World '
'Geodetic '
'System '
'1984 '
'(G1762)'},
{'id': {'authority': 'EPSG',
'code': 1309},
'name': 'World '
'Geodetic '
'System '
'1984 '
'(G2139)'},
{'id': {'authority': 'EPSG',
'code': 1383},
'name': 'World '
'Geodetic '
'System '
'1984 '
'(G2296)'}],
'name': 'World Geodetic '
'System 1984 '
'ensemble'},
'id': {'authority': 'EPSG', 'code': 4326},
'name': 'WGS 84'},
'bbox': {'east_longitude': 51.73,
'north_latitude': 83.67,
'south_latitude': 29.24,
'west_longitude': -42.52},
'conversion': {'method': {'id': {'authority': 'EPSG',
'code': 1125},
'name': 'Azimuthal Equidistant'},
'name': 'Equi7 projection - Europe',
'parameters': [{'id': {'authority': 'EPSG',
'code': 8801},
'name': 'Latitude of natural '
'origin',
'unit': 'degree',
'value': 53},
{'id': {'authority': 'EPSG',
'code': 8802},
'name': 'Longitude of natural '
'origin',
'unit': 'degree',
'value': 24},
{'id': {'authority': 'EPSG',
'code': 8806},
'name': 'False easting',
'unit': 'metre',
'value': 5837287.82},
{'id': {'authority': 'EPSG',
'code': 8807},
'name': 'False northing',
'unit': 'metre',
'value': 2121415.696}]},
'coordinate_system': {'axis': [{'abbreviation': 'E',
'direction': 'east',
'name': 'Easting',
'unit': 'metre'},
{'abbreviation': 'N',
'direction': 'north',
'name': 'Northing',
'unit': 'metre'}],
'subtype': 'Cartesian'},
'id': {'authority': 'EPSG', 'code': 27704},
'name': 'WGS 84 / Equi7 Europe',
'scope': 'Continental mapping of raster data.',
'type': 'ProjectedCRS'}},
'description': None,
'id': None,
'keywords': None,
'orderedAxes': None,
'tileMatrices': [{'cellSize': 500.0,
'cornerOfOrigin': 'topLeft',
'description': None,
'id': '0',
'keywords': None,
'matrixHeight': 11,
'matrixWidth': 14,
'pointOfOrigin': [0.0, 6000000.0],
'scaleDenominator': 1785714.285714286,
'tileHeight': 1200,
'tileWidth': 1200,
'title': 'T6',
'variableMatrixWidths': None},
{'cellSize': 500.0,
'cornerOfOrigin': 'topLeft',
'description': None,
'id': '1',
'keywords': None,
'matrixHeight': 21,
'matrixWidth': 28,
'pointOfOrigin': [0.0, 5700000.0],
'scaleDenominator': 1785714.285714286,
'tileHeight': 600,
'tileWidth': 600,
'title': 'T3',
'variableMatrixWidths': None}],
'title': None,
'uri': None,
'wellKnownScaleSet': None}
Visualisation#
A Equi7 tiling systems can be also visualised on a map in a similar manner as a Equi7Tile object, if the optional dependencies matplotlib and cartopy are installed.
import matplotlib.pyplot as plt
plt.figure(figsize=(20, 28))
_ = e7eu.plot(
tiling_id="T3",
label_tile=True,
label_size=5,
plot_zone=True,
facecolor="#2cee768f",
alpha=0.6,
extent=[-6e5, -6e5, 9e6, 6e6],
)
