Determining insert point of each WMTS tile

Question

I am implementing my own WMTS tile server. I am using tileserverphp.

I can download each tile that I need to my computer.

The next step is to insert the image into my CAD software (BricsCAD), however, I do not understand how to find the insert point for each tile.

(To clarify, I understand HOW to insert images in BricsCAD, but I don't know WHERE to insert them.)

The tile is a .png image.

Can someone help me figure this out?

Answer 1

You can work it out from the origin of the tile matrix + the tile X (or Y) times the tile size.

So given a tile matrix set:

  <TileMatrixSet>
      <ows:Identifier>EPSG:4326</ows:Identifier>
      <ows:SupportedCRS>urn:ogc:def:crs:EPSG::4326</ows:SupportedCRS>
      <TileMatrix>
        <ows:Identifier>EPSG:4326:0</ows:Identifier>
        <ScaleDenominator>2.795411320143589E8</ScaleDenominator>
        <TopLeftCorner>90.0 -180.0</TopLeftCorner>
        <TileWidth>256</TileWidth>
        <TileHeight>256</TileHeight>
        <MatrixWidth>2</MatrixWidth>
        <MatrixHeight>1</MatrixHeight>
      </TileMatrix>
      <TileMatrix>
        <ows:Identifier>EPSG:4326:1</ows:Identifier>
        <ScaleDenominator>1.3977056600717944E8</ScaleDenominator>
        <TopLeftCorner>90.0 -180.0</TopLeftCorner>
        <TileWidth>256</TileWidth>
        <TileHeight>256</TileHeight>
        <MatrixWidth>4</MatrixWidth>
        <MatrixHeight>2</MatrixHeight>
      </TileMatrix>

A tile in the first level EPSG:4326:0 with a X = 1 and Y = 0 will have it's top left corner at (-180 + tilesize*X, 90 - tilesize *Y) and tilesize is 256 * 2.8e-4 / 111319 where 111319 is a magic number for the number of metres in 60 Nautical Miles because we are using degrees here, for a projected CRS it would be the number of metres per unit of the projection.

Here's the same maths in Java:

private static final double PixelSizeMeters = 0.28e-3;
public static ReferencedEnvelope getExtentFromTileName(
        WMTSTileIdentifier tileIdentifier, TileService service) {
    WMTSZoomLevel zl = new WMTSZoomLevel(tileIdentifier.getZ(), (WMTSTileService) service);
    TileMatrix tileMatrix =
            ((WMTSTileService) service).getMatrixSet().getMatrices().get(zl.getZoomLevel());

    CoordinateReferenceSystem crs = tileMatrix.getCrs();
    CoordinateSystem coordinateSystem = crs.getCoordinateSystem();

    double pixelSpan = getPixelSpan(tileMatrix);
    double tileSpanY = (tileMatrix.getTileHeight() * pixelSpan);
    double tileSpanX = (tileMatrix.getTileWidth() * pixelSpan);

    double tileMatrixMinX;
    double tileMatrixMaxY;
    boolean longFirst = coordinateSystem.getAxis(0).getDirection().equals(AxisDirection.EAST);
    if (longFirst) {
        tileMatrixMinX = tileMatrix.getTopLeft().getX();
        tileMatrixMaxY = tileMatrix.getTopLeft().getY();
    } else {
        tileMatrixMaxY = tileMatrix.getTopLeft().getX();
        tileMatrixMinX = tileMatrix.getTopLeft().getY();
    }
    ReferencedEnvelope ret = new ReferencedEnvelope(crs);
    double minX = tileIdentifier.getX() * tileSpanX + tileMatrixMinX;
    double maxY = tileMatrixMaxY - tileIdentifier.getY() * tileSpanY;
    double maxX = minX + tileSpanX;
    double minY = maxY - tileSpanY;
    if (longFirst) {
        ret.expandToInclude(minX, minY);
        ret.expandToInclude(maxX, maxY);
    } else {
        ret.expandToInclude(minY, minX);
        ret.expandToInclude(maxY, maxX);
    }

    return ret;
}

/** */
private static double getPixelSpan(TileMatrix tileMatrix) {
    CoordinateSystem coordinateSystem = tileMatrix.getCrs().getCoordinateSystem();
    @SuppressWarnings("unchecked")
    Unit<Length> unit = (Unit<Length>) coordinateSystem.getAxis(0).getUnit();

    // now divide by meters per unit!
    double pixelSpan = tileMatrix.getDenominator() * PixelSizeMeters;
    if (unit.equals(NonSI.DEGREE_ANGLE)) {
        /*
         * use the length of a degree at the equator = 60 nautical miles!
         * unit = USCustomary.NAUTICAL_MILE; UnitConverter metersperunit =
         * unit.getConverterTo(SI.METRE); pixelSpan /=
         * metersperunit.convert(60.0);
         */

        // constant value from
        // https://msi.nga.mil/MSISiteContent/StaticFiles/Calculators/degree.html
        // apparently - 60.10764611706782 NaMiles
        pixelSpan /= 111319;
    } else {
        UnitConverter metersperunit = unit.getConverterTo(SI.METRE);
        pixelSpan /= metersperunit.convert(1);
    }
    return pixelSpan;
}

Answer 2

I have worked it out and will post here encase it helps someone else.

You use the Tile Column, Tile Row, Tile Span and the Tile Matrix MinX and Max Y values.

Here is psuedo code from http://docs.opengeospatial.org/is/17-083r2/17-083r2.html#107 Annex I.

The upper-left corner (leftX, upperY) of the tile:

leftX = tileCol * tileSpanX + tileMatrixMinX
upperY = tileMatrixMaxY - tileRow * tileSpanY

and the lower-right corner (rightX, lowerY) of the tile:

rightX = (tileCol+1) * tileSpanX + tileMatrixMinX
lowerY = tileMatrixMaxY – (tileRow+1) * tileSpanY

Answer 3

It's a while back when this question was asked. But, if anyone is interested in a completely open AutoLISP solution to import WMTS, WMS, TMS, etc. into AutoCAD and BricsCAD/Lite, give this a try:

https://github.com/cxcandid/GIS-Tools-for-CAD/tree/main/WMS

The solution utilizes GDAL to import raster image tiles.