Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[Backport release/3.10] MVT: allow generating tilesets with more than 1 tile at zoom level 0 #11785

Merged
merged 1 commit into from
Feb 3, 2025
Merged

[Backport release/3.10] MVT: allow generating tilesets with more than 1 tile at zoom level 0 #11785

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
44 changes: 43 additions & 1 deletion autotest/ogr/ogr_mvt.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1754,4 +1754,46 @@ def test_ogr_mvt_write_reuse_temp_db():


###############################################################################
#


@pytest.mark.require_driver("SQLite")
@pytest.mark.require_geos
@pytest.mark.parametrize(
"TILING_SCHEME",
["EPSG:4326,-180,90,180", "EPSG:4326,-180,90,180,2,1", "EPSG:4326,-180,90,180,1,1"],
)
def test_ogr_mvt_write_custom_tiling_scheme_WorldCRS84Quad(tmp_vsimem, TILING_SCHEME):

src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
srs = osr.SpatialReference()
srs.SetFromUserInput("WGS84")
lyr = src_ds.CreateLayer("mylayer", srs=srs)

f = ogr.Feature(lyr.GetLayerDefn())
f.SetGeometry(ogr.CreateGeometryFromWkt("POINT(120 40)"))
lyr.CreateFeature(f)

filename = str(tmp_vsimem / "out")
out_ds = gdal.VectorTranslate(
filename,
src_ds,
format="MVT",
datasetCreationOptions=["TILING_SCHEME=" + TILING_SCHEME],
)
assert out_ds is not None
out_ds = None

if TILING_SCHEME == "EPSG:4326,-180,90,180,1,1":
# If explictly setting tile_matrix_width_zoom_0 == 1,
# we have no tiles beyond longitude 0 degree
with pytest.raises(Exception):
ogr.Open(filename + "/0")
else:
out_ds = ogr.Open(filename + "/0")
assert out_ds is not None
out_lyr = out_ds.GetLayerByName("mylayer")
assert out_lyr.GetSpatialRef().ExportToWkt().find("4326") >= 0
out_f = out_lyr.GetNextFeature()
ogrtest.check_feature_geometry(
out_f, "MULTIPOINT ((120.0146484375 39.990234375))"
)
23 changes: 21 additions & 2 deletions doc/source/drivers/vector/mvt.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -144,7 +144,8 @@ scanning the features of the tile(s).

As an extension, OGR handles in reading and writing custom tiling
schemes by using the *crs*, *tile_origin_upper_left_x*,
*tile_origin_upper_left_y* and *tile_dimension_zoom_0* metadata items.
*tile_origin_upper_left_y*, *tile_dimension_zoom_0*, *tile_matrix_width_zoom_0*
and *tile_matrix_height_zoom_0* metadata items.
For example, for the Finnish ETRS-TM35FIN (EPSG:3067) tiling scheme:

.. code-block:: json
Expand All @@ -157,6 +158,21 @@ For example, for the Finnish ETRS-TM35FIN (EPSG:3067) tiling scheme:
"tile_dimension_zoom_0":2097152.0,
}
Or for a ``WorldCRS84Quad`` tiling scheme with 2 tiles in the horizontal
direction at zoom level 0:

.. code-block:: json
{
"...": "...",
"crs":"EPSG:4326",
"tile_origin_upper_left_x":-180.0,
"tile_origin_upper_left_y":90.0,
"tile_dimension_zoom_0":180.0,
"tile_matrix_width_zoom_0":2,
"tile_matrix_height_zoom_0":1
}
Opening options
---------------

Expand Down Expand Up @@ -351,7 +367,7 @@ The following dataset creation options are supported:
metadata item, which is the center of :co:`BOUNDS` at minimum zoom level.

- .. co:: TILING_SCHEME
:choices: <crs\,tile_origin_upper_left_x\,tile_origin_upper_left_y\,tile_dimension_zoom_0>
:choices: <crs\,tile_origin_upper_left_x\,tile_origin_upper_left_y\,tile_dimension_zoom_0[\,tile_matrix_width_zoom_0\,tile_matrix_height_zoom_0]>

Define a custom tiling scheme with a CRS
(typically given as EPSG:XXXX), the coordinates of the upper-left
Expand All @@ -365,6 +381,9 @@ The following dataset creation options are supported:
scheme, the 'crs', 'tile_origin_upper_left_x',
'tile_origin_upper_left_y' and 'tile_dimension_zoom_0' entries are
added to the metadata.json, and are honoured by the OGR MVT reader.
Starting with GDAL 3.10.2, 'tile_matrix_width_zoom_0' (resp.
'tile_matrix_height_zoom_0') can be specified to indicate the number of
tiles along the X (resp. Y) axis at zoom level 0.

Layer configuration
-------------------
Expand Down
132 changes: 106 additions & 26 deletions ogr/ogrsf_frmts/mvt/ogrmvtdataset.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -292,9 +292,14 @@ class OGRMVTDataset final : public GDALDataset
bool m_bClip = true;
CPLString m_osTileExtension{"pbf"};
OGRSpatialReference *m_poSRS = nullptr;
double m_dfTileDim0 = 0.0;
double m_dfTopXOrigin = 0.0;
double m_dfTopYOrigin = 0.0;
double m_dfTileDim0 =
0.0; // Extent (in CRS units) of a tile at zoom level 0
double m_dfTopXOrigin = 0.0; // top-left X of tile matrix scheme
double m_dfTopYOrigin = 0.0; // top-left Y of tile matrix scheme
int m_nTileMatrixWidth0 =
1; // Number of tiles along X axis at zoom level 0
int m_nTileMatrixHeight0 =
1; // Number of tiles along Y axis at zoom level 0

static GDALDataset *OpenDirectory(GDALOpenInfo *);

Expand Down Expand Up @@ -322,20 +327,30 @@ class OGRMVTDataset final : public GDALDataset
return m_poSRS;
}

double GetTileDim0() const
inline double GetTileDim0() const
{
return m_dfTileDim0;
}

double GetTopXOrigin() const
inline double GetTopXOrigin() const
{
return m_dfTopXOrigin;
}

double GetTopYOrigin() const
inline double GetTopYOrigin() const
{
return m_dfTopYOrigin;
}

inline int GetTileMatrixWidth0() const
{
return m_nTileMatrixWidth0;
}

inline int GetTileMatrixHeight0() const
{
return m_nTileMatrixHeight0;
}
};

/************************************************************************/
Expand Down Expand Up @@ -1746,8 +1761,10 @@ void OGRMVTDirectoryLayer::SetSpatialFilter(OGRGeometry *poGeomIn)

if (sEnvelope.IsInit() && sEnvelope.MinX >= -10 * m_poDS->GetTileDim0() &&
sEnvelope.MinY >= -10 * m_poDS->GetTileDim0() &&
sEnvelope.MaxX <= 10 * m_poDS->GetTileDim0() &&
sEnvelope.MaxY <= 10 * m_poDS->GetTileDim0())
sEnvelope.MaxX <=
10 * m_poDS->GetTileDim0() * m_poDS->GetTileMatrixWidth0() &&
sEnvelope.MaxY <=
10 * m_poDS->GetTileDim0() * m_poDS->GetTileMatrixHeight0())
{
const double dfTileDim = m_poDS->GetTileDim0() / (1 << m_nZ);
m_nFilterMinX = std::max(
Expand All @@ -1759,18 +1776,30 @@ void OGRMVTDirectoryLayer::SetSpatialFilter(OGRGeometry *poGeomIn)
m_nFilterMaxX = std::min(
static_cast<int>(
ceil((sEnvelope.MaxX - m_poDS->GetTopXOrigin()) / dfTileDim)),
(1 << m_nZ) - 1);
static_cast<int>(std::min<int64_t>(
INT_MAX, (static_cast<int64_t>(1) << m_nZ) *
m_poDS->GetTileMatrixWidth0() -
1)));
m_nFilterMaxY = std::min(
static_cast<int>(
ceil((m_poDS->GetTopYOrigin() - sEnvelope.MinY) / dfTileDim)),
(1 << m_nZ) - 1);
static_cast<int>(std::min<int64_t>(
INT_MAX, (static_cast<int64_t>(1) << m_nZ) *
m_poDS->GetTileMatrixHeight0() -
1)));
}
else
{
m_nFilterMinX = 0;
m_nFilterMinY = 0;
m_nFilterMaxX = (1 << m_nZ) - 1;
m_nFilterMaxY = (1 << m_nZ) - 1;
m_nFilterMaxX = static_cast<int>(
std::min<int64_t>(INT_MAX, (static_cast<int64_t>(1) << m_nZ) *
m_poDS->GetTileMatrixWidth0() -
1));
m_nFilterMaxY = static_cast<int>(
std::min<int64_t>(INT_MAX, (static_cast<int64_t>(1) << m_nZ) *
m_poDS->GetTileMatrixHeight0() -
1));
}
}

Expand Down Expand Up @@ -2429,6 +2458,7 @@ static bool LoadMetadata(const CPLString &osMetadataFile,
CPLJSONArray &oTileStatLayers, CPLJSONObject &oBounds,
OGRSpatialReference *poSRS, double &dfTopX,
double &dfTopY, double &dfTileDim0,
int &nTileMatrixWidth0, int &nTileMatrixHeight0,
const CPLString &osMetadataMemFilename)

{
Expand All @@ -2451,17 +2481,36 @@ static bool LoadMetadata(const CPLString &osMetadataFile,
if (!bLoadOK)
return false;

CPLJSONObject oCrs(oDoc.GetRoot().GetObj("crs"));
CPLJSONObject oTopX(oDoc.GetRoot().GetObj("tile_origin_upper_left_x"));
CPLJSONObject oTopY(oDoc.GetRoot().GetObj("tile_origin_upper_left_y"));
CPLJSONObject oTileDim0(oDoc.GetRoot().GetObj("tile_dimension_zoom_0"));
const CPLJSONObject oCrs(oDoc.GetRoot().GetObj("crs"));
const CPLJSONObject oTopX(
oDoc.GetRoot().GetObj("tile_origin_upper_left_x"));
const CPLJSONObject oTopY(
oDoc.GetRoot().GetObj("tile_origin_upper_left_y"));
const CPLJSONObject oTileDim0(
oDoc.GetRoot().GetObj("tile_dimension_zoom_0"));
nTileMatrixWidth0 = 1;
nTileMatrixHeight0 = 1;
if (oCrs.IsValid() && oTopX.IsValid() && oTopY.IsValid() &&
oTileDim0.IsValid())
{
poSRS->SetFromUserInput(oCrs.ToString().c_str());
dfTopX = oTopX.ToDouble();
dfTopY = oTopY.ToDouble();
dfTileDim0 = oTileDim0.ToDouble();
const CPLJSONObject oTMWidth0(
oDoc.GetRoot().GetObj("tile_matrix_width_zoom_0"));
if (oTMWidth0.GetType() == CPLJSONObject::Type::Integer)
nTileMatrixWidth0 = std::max(1, oTMWidth0.ToInteger());

const CPLJSONObject oTMHeight0(
oDoc.GetRoot().GetObj("tile_matrix_height_zoom_0"));
if (oTMHeight0.GetType() == CPLJSONObject::Type::Integer)
nTileMatrixHeight0 = std::max(1, oTMHeight0.ToInteger());

// Assumes WorldCRS84Quad with 2 tiles in width
// cf https://github.com/OSGeo/gdal/issues/11749
if (!oTMWidth0.IsValid() && dfTopX == -180 && dfTileDim0 == 180)
nTileMatrixWidth0 = 2;
}

oVectorLayers.Deinit();
Expand Down Expand Up @@ -2796,7 +2845,8 @@ GDALDataset *OGRMVTDataset::OpenDirectory(GDALOpenInfo *poOpenInfo)
if (!LoadMetadata(osMetadataFile, osMetadataContent, oVectorLayers,
oTileStatLayers, oBounds, poDS->m_poSRS,
poDS->m_dfTopXOrigin, poDS->m_dfTopYOrigin,
poDS->m_dfTileDim0, osMetadataMemFilename))
poDS->m_dfTileDim0, poDS->m_nTileMatrixWidth0,
poDS->m_nTileMatrixHeight0, osMetadataMemFilename))
{
delete poDS;
return nullptr;
Expand Down Expand Up @@ -3126,7 +3176,8 @@ GDALDataset *OGRMVTDataset::Open(GDALOpenInfo *poOpenInfo, bool bRecurseAllowed)
LoadMetadata(osMetadataFile, CPLString(), oVectorLayers,
oTileStatLayers, oBounds, poDS->m_poSRS,
poDS->m_dfTopXOrigin, poDS->m_dfTopYOrigin,
poDS->m_dfTileDim0, CPLString());
poDS->m_dfTileDim0, poDS->m_nTileMatrixWidth0,
poDS->m_nTileMatrixHeight0, CPLString());
}

const char *pszGeorefTopX =
Expand Down Expand Up @@ -3155,8 +3206,10 @@ GDALDataset *OGRMVTDataset::Open(GDALOpenInfo *poOpenInfo, bool bRecurseAllowed)
int nX = atoi(osX);
int nY = atoi(osY);
int nZ = atoi(osZ);
if (nZ >= 0 && nZ < 30 && nX >= 0 && nX < (1 << nZ) && nY >= 0 &&
nY < (1 << nZ))
if (nZ >= 0 && nZ < 30 && nX >= 0 &&
nX < (static_cast<int64_t>(1) << nZ) * poDS->m_nTileMatrixWidth0 &&
nY >= 0 &&
nY < (static_cast<int64_t>(1) << nZ) * poDS->m_nTileMatrixHeight0)
{
poDS->m_bGeoreferenced = true;
poDS->m_dfTileDimX = poDS->m_dfTileDim0 / (1 << nZ);
Expand Down Expand Up @@ -3325,6 +3378,10 @@ class OGRMVTWriterDataset final : public GDALDataset
double m_dfTopX = 0.0;
double m_dfTopY = 0.0;
double m_dfTileDim0 = 0.0;
int m_nTileMatrixWidth0 =
1; // Number of tiles along X axis at zoom level 0
int m_nTileMatrixHeight0 =
1; // Number of tiles along Y axis at zoom level 0
bool m_bReuseTempFile = false; // debug only

OGRErr PreGenerateForTile(
Expand Down Expand Up @@ -5516,6 +5573,10 @@ bool OGRMVTWriterDataset::GenerateMetadata(
oRoot);
WriteMetadataItem("tile_dimension_zoom_0", m_dfTileDim0, m_hDBMBTILES,
oRoot);
WriteMetadataItem("tile_matrix_width_zoom_0", m_nTileMatrixWidth0,
m_hDBMBTILES, oRoot);
WriteMetadataItem("tile_matrix_height_zoom_0", m_nTileMatrixHeight0,
m_hDBMBTILES, oRoot);
}

CPLJSONDocument oJsonDoc;
Expand Down Expand Up @@ -5818,11 +5879,17 @@ OGRErr OGRMVTWriterDataset::WriteFeature(OGRMVTWriterLayer *poLayer,
const int nTileMaxX =
std::min(static_cast<int>((sExtent.MaxX - m_dfTopX + dfBuffer) /
dfTileDim),
(1 << nZ) - 1);
static_cast<int>(std::min<int64_t>(
INT_MAX, (static_cast<int64_t>(1) << nZ) *
m_nTileMatrixWidth0 -
1)));
const int nTileMaxY =
std::min(static_cast<int>((m_dfTopY - sExtent.MinY + dfBuffer) /
dfTileDim),
(1 << nZ) - 1);
static_cast<int>(std::min<int64_t>(
INT_MAX, (static_cast<int64_t>(1) << nZ) *
m_nTileMatrixHeight0 -
1)));
for (int iX = nTileMinX; iX <= nTileMaxX; iX++)
{
for (int iY = nTileMinY; iY <= nTileMaxY; iY++)
Expand Down Expand Up @@ -6168,20 +6235,32 @@ GDALDataset *OGRMVTWriterDataset::Create(const char *pszFilename, int nXSize,
return nullptr;
}

CPLStringList aoList(CSLTokenizeString2(pszTilingScheme, ",", 0));
if (aoList.Count() == 4)
const CPLStringList aoList(CSLTokenizeString2(pszTilingScheme, ",", 0));
if (aoList.Count() >= 4)
{
poDS->m_poSRS->SetFromUserInput(aoList[0]);
poDS->m_dfTopX = CPLAtof(aoList[1]);
poDS->m_dfTopY = CPLAtof(aoList[2]);
poDS->m_dfTileDim0 = CPLAtof(aoList[3]);
if (aoList.Count() == 6)
{
poDS->m_nTileMatrixWidth0 = std::max(1, atoi(aoList[4]));
poDS->m_nTileMatrixHeight0 = std::max(1, atoi(aoList[5]));
}
else if (poDS->m_dfTopX == -180 && poDS->m_dfTileDim0 == 180)
{
// Assumes WorldCRS84Quad with 2 tiles in width
// cf https://github.com/OSGeo/gdal/issues/11749
poDS->m_nTileMatrixWidth0 = 2;
}
}
else
{
CPLError(CE_Failure, CPLE_AppDefined,
"Wrong format for TILING_SCHEME. "
"Expecting EPSG:XXXX,tile_origin_upper_left_x,"
"tile_origin_upper_left_y,tile_dimension_zoom_0");
"tile_origin_upper_left_y,tile_dimension_zoom_0[,tile_"
"matrix_width_zoom_0,tile_matrix_height_zoom_0]");
delete poDS;
return nullptr;
}
Expand Down Expand Up @@ -6334,7 +6413,8 @@ void RegisterOGRMVT()
" <Option name='TILING_SCHEME' type='string' "
"description='Custom tiling scheme with following format "
"\"EPSG:XXXX,tile_origin_upper_left_x,tile_origin_upper_left_y,"
"tile_dimension_zoom_0\"'/>"
"tile_dimension_zoom_0[,tile_matrix_width_zoom_0,tile_matrix_height_"
"zoom_0]\"'/>"
"</CreationOptionList>");
#endif // HAVE_MVT_WRITE_SUPPORT

Expand Down
Loading