diff options
| -rw-r--r-- | config.yml | 26 | ||||
| -rwxr-xr-x | webmap-import | 45 |
2 files changed, 44 insertions, 27 deletions
@@ -2,19 +2,21 @@ # Spatial Reference System SRS: 'EPSG:3006' # a.k.a. SWEREF99 TM, cf. https://epsg.io/3006 -# (2D) extent +# (2D) extent: minX, minY, maxX, maxY +# +# Lantmäteriet uses a tile-scheme where the origin (upper-left corner) is at N8500000 +# E-1200000 (SWEREF99 TM), where each tile is 256×256 pixels, and where the resolution at +# level 0 is 4096m per pixel. +# +# https://www.lantmateriet.se/globalassets/geodata/geodatatjanster/tb_twk_visning_cache_v1.1.0.pdf +# https://www.lantmateriet.se/globalassets/geodata/geodatatjanster/tb_twk_visning-oversiktlig_v1.0.3.pdf +# +# We set the extent to a 4×4 tiles square at level 2 somehow centered on Norrbotten and +# Västerbotten. This represent a TILEROW (x) offset of 5, and a TILECOL (y) offset of 2. +# A 4×8 tiles rectangle with the same upper-left and upper-right coordinates can be used +# to cover the entire country. +# extent: - # Lantmäteriet uses a tile-scheme where the origin (upper-left corner) is at N8500000 - # E-1200000 (SWEREF99 TM), where each tile is 256×256 pixels, and where the resolution - # at level 0 is 4096m per pixel. - # - # https://www.lantmateriet.se/globalassets/geodata/geodatatjanster/tb_twk_visning_cache_v1.1.0.pdf - # https://www.lantmateriet.se/globalassets/geodata/geodatatjanster/tb_twk_visning-oversiktlig_v1.0.3.pdf - # - # We set the extent to a 4×4 tiles square at level 2 somehow centered on Norrbotten and - # Västerbotten. This represent a TILEROW (x) offset of 5, and a TILECOL (y) offset of 2. - # A 4×8 tiles rectangle with the same upper-left and upper-right coordinates can be used - # to cover the entire country. - 110720 - 6927136 # alternatively 5878560 for the entire country - 1159296 diff --git a/webmap-import b/webmap-import index 0d65db3..3242b3a 100755 --- a/webmap-import +++ b/webmap-import @@ -431,12 +431,18 @@ def getSRS(srs_str): logging.debug('Default SRS: "%s" (%s)', srs.ExportToProj4(), srs.GetName()) return srs -def getExtent(extent): +# Convert extent [minX, minY, maxX, maxY] into a polygon and assign the +# given SRS. Like apps/ogr2ogr_lib.cpp, we segmentize the polygon to +# make sure it is sufficiently densified when transforming to source +# layer SRS for spatial filtering. +def getExtent(extent, srs=None): if extent is None: return if not (isinstance(extent, list) or isinstance(extent, tuple)) or len(extent) != 4: raise Exception(f'Invalid extent {extent}') + elif srs is None: + raise Exception('Configured extent but no SRS') logging.debug('Configured extent in %s: %s', srs.GetName(), ', '.join(map(str, extent))) @@ -451,6 +457,21 @@ def getExtent(extent): polygon = ogr.Geometry(ogr.wkbPolygon) polygon.AddGeometry(ring) + # we expressed extent as minX, minY, maxX, maxY (easting/northing + # ordered, i.e., in traditional GIS order) + srs2 = srs.Clone() + srs2.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) + polygon.AssignSpatialReference(srs2) + polygon.TransformTo(srs) + + segment_distance_metre = 10 * 1000 + if srs.IsGeographic(): + dfMaxLength = segment_distance_metre / math.radians(srs.GetSemiMajor()) + polygon.Segmentize(dfMaxLength) + elif srs.IsProjected(): + dfMaxLength = segment_distance_metre / srs.GetLinearUnits() + polygon.Segmentize(dfMaxLength) + return polygon # Validate the output layer against the provided SRS and creation options @@ -866,6 +887,9 @@ def importSource2(lyr_dst, path, args={}, basedir=None, extent=None): canIgnoreFields = lyr.TestCapability(ogr.OLCIgnoreFields) srs = lyr.GetSpatialRef() + if srs is None: + raise Exception('Source layer has no SRS') + srs_dst = lyr_dst.GetSpatialRef() if srs_dst.IsSame(srs): logging.debug('Both source and destination have the same SRS (%s), skipping coordinate transformation', @@ -914,19 +938,10 @@ def importSource2(lyr_dst, path, args={}, basedir=None, extent=None): count1 = -1 if args.get('spatial-filter', True) and extent is not None: - extent2 = extent.Clone() - extent2.AssignSpatialReference(srs_dst) - - if srs is not None and not srs.IsSame(srs_dst): - # make sure to densify it sufficiently to avoid issues, cf. apps/ogr2ogr_lib.cpp - segment_distance_metre = 10 * 1000 - if srs_dst.IsGeographic(): - dfMaxLength = segment_distance_metre / math.radians(srs_dst.GetSemiMajor()) - extent2.Segmentize(dfMaxLength) - elif srs_dst.IsProjected(): - dfMaxLength = segment_distance_metre / srs_dst.GetLinearUnits() - extent2.Segmentize(dfMaxLength) - + if extent.GetSpatialReference().IsSame(srs): + extent2 = extent + else: + extent2 = extent.Clone() if extent2.TransformTo(srs) != ogr.OGRERR_NONE: raise Exception(f'Could not transform extent {extent.ExportToWkt()} to {srs.GetName()}') @@ -1033,7 +1048,7 @@ if __name__ == '__main__': # get configured Spatial Reference System and extent srs = getSRS(common.config.get('SRS', None)) - extent = getExtent(common.config.get('extent', None)) + extent = getExtent(common.config.get('extent', None), srs=srs) cachedir = Path(args.cachedir) if args.cachedir is not None else None rv = 0 |