Tiling system#
pytileproj’s tiling_system module defines several classes representing a multi-level (e.g., zoom, tiling level, …) projected tiling scheme. A tiling system is a collection of tilings defined in pytileproj.tiling, i.e., regular and irregular tilings. These tilings can then be ordered by their tiling/zoom level and forced to only allow certain pixel spacings at a specific level. In the following section, we will look at the most common use case of a tiling system, i.e., a regular projected tiling system.
Regular projected tiling system#
pytileproj’s regular tiling system class RegularProjTilingSystem exptects the following arguments:
name: Name of the tiling system.crs: A spatial reference system represented by anything pyproj supports, e.g., EPSG code, PROJ4 string, WKT string, …tilings: Dictionary linking tiling/zoom levels to regular tilings.proi_zone_geog(optional): Projection zone in geographic coordinates. If not given, the zone is fetched from the EPSG registry.
Initialisation#
Lets dive deeper into how a RegularProjTilingSystem can be used. As an example, we want to create tiling system for the European Equi7Grid projection. First, we need to define the tilings we want to use within the system:
from pytileproj.tiling import RegularTiling
extent = (0, 0, 8_700_000, 6_000_000)
axis_orientation = ("E", "S")
coarse_tiling = RegularTiling(
name="my_coarse_tiling",
extent=extent,
sampling=1_000,
tile_shape=(300_000, 300_000),
tiling_level=1,
axis_orientation=axis_orientation,
)
fine_tiling = RegularTiling(
name="my_fine_tiling",
extent=extent,
sampling=10,
tile_shape=(100_000, 100_000),
tiling_level=2,
axis_orientation=axis_orientation,
)
Then, we can create a RegularProjTilingSystem by attaching the respective EPSG code:
from pathlib import Path
from pytileproj import RegularProjTilingSystem
name = "e7eu"
epsg = 27704
# this argument is only optional, but in case the EPSG API is down,
# it needs to be provided
proj_zone_geog = Path("../../tests/data/eu_zone.parquet")
rpts = RegularProjTilingSystem(
name=name,
tilings={tiling.tiling_level: tiling for tiling in [coarse_tiling, fine_tiling]},
crs=epsg,
proj_zone_geog=proj_zone_geog,
)
print(rpts) # noqa: T201
RegularProjTilingSystem
-----------------------
Name:
e7eu
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:
{1: RegularTiling(my_coarse_tiling, 1000.0), 2: RegularTiling(my_fine_tiling, 10.0)}
RegularProjTilingSystem also provides another method to create an object from tiling definitions, which directly assume that a regular projected tiling system has the same extent and axis orientation. In addition, the user has control over the tile size in projection units rather than the tile shape in pixels:
from pytileproj.tiling_system import ProjSystemDefinition, RegularTilingDefinition
rpts_def = ProjSystemDefinition(
name=name,
crs=epsg,
min_xy=extent[:2],
max_xy=extent[2:],
axis_orientation=axis_orientation,
proj_zone_geog=proj_zone_geog,
)
tiling_defs = {
1: RegularTilingDefinition(name="my_coarse_tiling", tile_shape=300_000),
2: RegularTilingDefinition(name="my_fine_tiling", tile_shape=100_000),
}
rpts = RegularProjTilingSystem.from_sampling(
{1: 1_000, 2: 10},
proj_def=rpts_def,
tiling_defs=tiling_defs,
)
Projection system interactions#
Often it is required to transform coordinates back and forth between different projection systems. RegularProjTilingSystem 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 = rpts.lonlat_to_xy(lon, lat)
proj_coord
ProjCoord(5271748.127753316, 1613228.4676009014)
The same can be done the other way around:
rpts.xy_to_lonlat(proj_coord.x, proj_coord.y)
GeogCoord(16.392599999999998, 48.1674)
Note that for coordinates outside the projection zone, an error will be raised:
lon, lat = -100, 30
rpts.lonlat_to_xy(lon, lat)
---------------------------------------------------------------------------
GeomOutOfZoneError Traceback (most recent call last)
Cell In[7], line 2
1 lon, lat = -100, 30
----> 2 rpts.lonlat_to_xy(lon, lat)
File ~/work/pytileproj/pytileproj/src/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:
import pyproj
from pytileproj import ProjCoord
coord = ProjCoord(lon, lat, pyproj.CRS.from_epsg(4326))
coord in rpts
False
If you want to know the units of the projection system, you can retrieve this information via:
rpts.unit
'metre'
Tiles#
Regular projected tiling systems create tiles on the fly, since their regular tiling scheme allows an efficient computation of tile properties. If you want to know if the tilings within the tiling systems are congruent, i.e., tiles at a higher tiling level (fine tiling) are multiples of tiles at a lower tiling level (coarse tiling), you can use the is_congruent property.
rpts.is_congruent
True
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
rpts.get_tile_from_lonlat(lon, lat, tiling_id=1)
RasterTile(e7eu_X17Y14T01)
In addition, you can create a tile by using an OGC tile index (x, y, z):
rpts.get_tile_from_index(17, 5, 1)
RasterTile(e7eu_X17Y05T01)
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(rpts.get_tiles_in_geog_bbox((16, 48, 18, 50), tiling_id="my_fine_tiling"))
[RasterTile(e7eu_X52Y41T02),
RasterTile(e7eu_X52Y42T02),
RasterTile(e7eu_X52Y43T02),
RasterTile(e7eu_X52Y44T02),
RasterTile(e7eu_X53Y41T02),
RasterTile(e7eu_X53Y42T02),
RasterTile(e7eu_X53Y43T02),
RasterTile(e7eu_X53Y44T02),
RasterTile(e7eu_X54Y41T02),
RasterTile(e7eu_X54Y42T02)]
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 native coordinates
list(
rpts.get_tiles_in_bbox(
(5155982, 1505898, 6323863, 2140396), tiling_id="my_coarse_tiling"
)
)
[RasterTile(e7eu_X17Y12T01),
RasterTile(e7eu_X17Y13T01),
RasterTile(e7eu_X17Y14T01),
RasterTile(e7eu_X18Y12T01),
RasterTile(e7eu_X18Y13T01),
RasterTile(e7eu_X18Y14T01),
RasterTile(e7eu_X19Y12T01),
RasterTile(e7eu_X19Y13T01),
RasterTile(e7eu_X19Y14T01),
RasterTile(e7eu_X20Y12T01),
RasterTile(e7eu_X20Y13T01),
RasterTile(e7eu_X20Y14T01),
RasterTile(e7eu_X21Y12T01),
RasterTile(e7eu_X21Y13T01),
RasterTile(e7eu_X21Y14T01)]
or by using a geometry.
from shapely import Polygon
from pytileproj import GeogGeom
geom = Polygon([(16, 48), (18, 48), (18, 50), (16, 50)])
ggeom = GeogGeom(geom=geom)
list(rpts.get_tiles_in_geom(ggeom, tiling_id="my_fine_tiling"))
[RasterTile(e7eu_X52Y41T02),
RasterTile(e7eu_X52Y42T02),
RasterTile(e7eu_X52Y43T02),
RasterTile(e7eu_X52Y44T02),
RasterTile(e7eu_X53Y41T02),
RasterTile(e7eu_X53Y42T02),
RasterTile(e7eu_X53Y43T02),
RasterTile(e7eu_X53Y44T02),
RasterTile(e7eu_X54Y41T02),
RasterTile(e7eu_X54Y42T02)]
Note
All geometries passed to pytileproj functions and class methods are segmentised to ensure an adequate representation of the geometry in a different projection system. A maximum segment length is defined as
0.1 degree for geographic projections
10 kilometres for metric projections
Since tiling systems represent an hierarchy of different tilings, pytileproj also allows to trace parent and children tiles. You can obtain the parent of a tile with
rpts.get_parent_from_name("X53Y17T02")
RasterTile(e7eu_X17Y05T01)
and the children with
list(rpts.get_children_from_name("X17Y05T01"))
[RasterTile(e7eu_X51Y15T02),
RasterTile(e7eu_X51Y16T02),
RasterTile(e7eu_X51Y17T02),
RasterTile(e7eu_X52Y15T02),
RasterTile(e7eu_X52Y16T02),
RasterTile(e7eu_X52Y17T02),
RasterTile(e7eu_X53Y15T02),
RasterTile(e7eu_X53Y16T02),
RasterTile(e7eu_X53Y17T02)]
Export#
Note
The export methods requires additional depencies, which are not installed by default. You can install them via uv pip install pytileproj[geo].
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:
rpts.to_geodataframe()
| crs | n_rows | n_cols | geotrans | px_origin | name | tiling_name | tiling_level | geometry | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 27704 | 300 | 300 | (0.0, 1000.0, 0.0, 2700000.0, 0.0, -1000.0) | ul | e7eu_X00Y11T01 | my_coarse_tiling | 1 | POLYGON ((0 2400000, 300000 2400000, 300000 27... |
| 1 | 27704 | 300 | 300 | (0.0, 1000.0, 0.0, 2400000.0, 0.0, -1000.0) | ul | e7eu_X00Y12T01 | my_coarse_tiling | 1 | POLYGON ((0 2100000, 300000 2100000, 300000 24... |
| 2 | 27704 | 300 | 300 | (0.0, 1000.0, 0.0, 2100000.0, 0.0, -1000.0) | ul | e7eu_X00Y13T01 | my_coarse_tiling | 1 | POLYGON ((0 1800000, 300000 1800000, 300000 21... |
| 3 | 27704 | 300 | 300 | (300000.0, 1000.0, 0.0, 3000000.0, 0.0, -1000.0) | ul | e7eu_X01Y10T01 | my_coarse_tiling | 1 | POLYGON ((300000 2700000, 600000 2700000, 6000... |
| 4 | 27704 | 300 | 300 | (300000.0, 1000.0, 0.0, 2700000.0, 0.0, -1000.0) | ul | e7eu_X01Y11T01 | my_coarse_tiling | 1 | POLYGON ((300000 2400000, 600000 2400000, 6000... |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 3098 | 27704 | 10000 | 10000 | (8100000.0, 10.0, 0.0, 1100000.0, 0.0, -10.0) | ul | e7eu_X81Y49T02 | my_fine_tiling | 2 | POLYGON ((8100000 1000000, 8200000 1000000, 82... |
| 3099 | 27704 | 10000 | 10000 | (8100000.0, 10.0, 0.0, 1000000.0, 0.0, -10.0) | ul | e7eu_X81Y50T02 | my_fine_tiling | 2 | POLYGON ((8100000 900000, 8200000 900000, 8200... |
| 3100 | 27704 | 10000 | 10000 | (8100000.0, 10.0, 0.0, 900000.0, 0.0, -10.0) | ul | e7eu_X81Y51T02 | my_fine_tiling | 2 | POLYGON ((8100000 800000, 8200000 800000, 8200... |
| 3101 | 27704 | 10000 | 10000 | (8200000.0, 10.0, 0.0, 1000000.0, 0.0, -10.0) | ul | e7eu_X82Y50T02 | my_fine_tiling | 2 | POLYGON ((8200000 900000, 8300000 900000, 8300... |
| 3102 | 27704 | 10000 | 10000 | (8200000.0, 10.0, 0.0, 900000.0, 0.0, -10.0) | ul | e7eu_X82Y51T02 | my_fine_tiling | 2 | POLYGON ((8200000 800000, 8300000 800000, 8300... |
3103 rows × 9 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_tiling_system.json")
rpts.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': 1000.0,
'cornerOfOrigin': 'topLeft',
'description': None,
'id': '1',
'keywords': None,
'matrixHeight': 20,
'matrixWidth': 29,
'pointOfOrigin': [0.0, 6000000.0],
'scaleDenominator': 3571428.571428572,
'tileHeight': 300,
'tileWidth': 300,
'title': 'my_coarse_tiling',
'variableMatrixWidths': None},
{'cellSize': 10.0,
'cornerOfOrigin': 'topLeft',
'description': None,
'id': '2',
'keywords': None,
'matrixHeight': 60,
'matrixWidth': 87,
'pointOfOrigin': [0.0, 6000000.0],
'scaleDenominator': 35714.28571428572,
'tileHeight': 10000,
'tileWidth': 10000,
'title': 'my_fine_tiling',
'variableMatrixWidths': None}],
'title': None,
'uri': None,
'wellKnownScaleSet': None}
Visualisation#
A RegularProjTilingSystem can be also visualised on a map in a similar manner as a RasterTile object, if the optional dependencies matplotlib and cartopy are installed.
import matplotlib.pyplot as plt
plt.figure(figsize=(20, 28))
_ = rpts.plot(
tiling_id=1,
label_tile=True,
label_size=5,
plot_zone=True,
facecolor="#2cee768f",
alpha=0.6,
extent=extent,
)
