Python mercantile.bounds() Examples

def pyramid_jobs(x, y, z, job_zoom, max_zoom):
    """
    Generate pyramid jobs for a given job_zoom level
    starting at with the parent tile defined by x, y, z.
    """
    if z == job_zoom:
        bounds = mercantile.bounds(x, y, z)
        yield create_pyramid_job(
            x=x, y=y,
            min_zoom=z, max_zoom=max_zoom,
            bounds=bounds
        )
        return

    tiles = all_descendant_tiles(x, y, z, job_zoom)
    pyramid_zoom_level_tiles = (t for t in tiles if t.z == job_zoom)

    for tile in pyramid_zoom_level_tiles:
        bounds = mercantile.bounds(tile.x, tile.y, tile.z)
        yield create_pyramid_job(tile.x, tile.y, min_zoom=tile.z,
                                 max_zoom=max_zoom, bounds=bounds) 

def test_tile_exists_valid():
    """Should work as expected (return true)."""
    bounds = [-80, 34, -75, 40]
    # Contains
    assert utils.tile_exists(bounds, 7, 36, 50)  # bounds contains tile bounds
    assert utils.tile_exists(bounds, 3, 2, 3)  # tile bounds contains bounds

    # Intersects
    assert utils.tile_exists(bounds, 7, 35, 50)
    assert utils.tile_exists(bounds, 7, 37, 50)
    assert utils.tile_exists(bounds, 7, 36, 51)
    assert utils.tile_exists(bounds, 7, 37, 51)
    assert utils.tile_exists(bounds, 7, 35, 51)
    assert utils.tile_exists(bounds, 7, 35, 48)
    assert utils.tile_exists(bounds, 7, 37, 48)

    # Outside tiles
    assert not utils.tile_exists(bounds, 7, 36, 40)
    assert not utils.tile_exists(bounds, 7, 36, 60)
    assert not utils.tile_exists(bounds, 7, 25, 50)
    assert not utils.tile_exists(bounds, 7, 70, 50) 

def test_tile_read_vrt_option():
    """Should work as expected (read landsat band)."""
    bounds = (
        -8844681.416934313,
        3757032.814272982,
        -8766409.899970293,
        3835304.331237001,
    )
    tilesize = 16
    with rasterio.open(f"{LANDSAT_PATH}_B2.TIF") as src_dst:
        arr, mask = reader.part(
            src_dst,
            bounds,
            tilesize,
            tilesize,
            warp_vrt_option=dict(source_extra=10, num_threads=10),
        )
    assert arr.shape == (1, 16, 16)
    assert mask.shape == (16, 16) 

def pixel_to_location(tile, dx, dy):
    """Converts a pixel in a tile to a coordinate.

    Args:
      tile: the mercantile tile to calculate the location in.
      dx: the relative x offset in range [0, 1].
      dy: the relative y offset in range [0, 1].

    Returns:
      The coordinate for the pixel in the tile.
    """

    assert 0 <= dx <= 1, "x offset is in [0, 1]"
    assert 0 <= dy <= 1, "y offset is in [0, 1]"

    west, south, east, north = mercantile.bounds(tile)

    def lerp(a, b, c):
        return a + c * (b - a)

    lon = lerp(west, east, dx)
    lat = lerp(south, north, dy)

    return lon, lat 

def test_get_vrt_transform_valid():
    """Should return correct transform and size."""
    bounds = (
        -11663507.036777973,
        4715018.0897710975,
        -11663487.927520901,
        4715037.199028169,
    )

    with rasterio.open(S3_PATH) as src:
        vrt_transform, vrt_width, vrt_height = utils.get_vrt_transform(
            src, bounds, 64, 64
        )
        assert vrt_transform[2] == -11663507.036777973
        assert vrt_transform[5] == 4715037.199028169
        assert vrt_width == 100
        assert vrt_height == 100

        vrt_transform, vrt_width, vrt_height = utils.get_vrt_transform(
            src, bounds, 256, 256
        )
        assert vrt_transform[2] == -11663507.036777973
        assert vrt_transform[5] == 4715037.199028169
        assert vrt_width == 256
        assert vrt_height == 256 

def test_get_vrt_transform_valid4326():
    """Should return correct transform and size."""
    bounds = (
        -104.77523803710938,
        38.95353532141205,
        -104.77455139160156,
        38.954069293441066,
    )
    with rasterio.open(S3_PATH) as src:
        vrt_transform, vrt_width, vrt_height = utils.get_vrt_transform(
            src, bounds, 256, 256, dst_crs=constants.WGS84_CRS
        )

    assert vrt_transform[2] == -104.77523803710938
    assert vrt_transform[5] == 38.954069293441066
    assert vrt_width == 420
    assert vrt_height == 327 

def test_aligned_with_internaltile():
    """Check if COG is in WebMercator and aligned with internal tiles."""
    bounds = mercantile.bounds(43, 25, 7)
    with rasterio.open(COG_DST) as src_dst:
        assert not utils._requested_tile_aligned_with_internal_tile(
            src_dst, bounds, 256, 256
        )

    with rasterio.open(NOCOG) as src_dst:
        assert not utils._requested_tile_aligned_with_internal_tile(
            src_dst, bounds, 256, 256
        )

    bounds = mercantile.bounds(147, 182, 9)
    with rasterio.open(COG_NOWEB) as src_dst:
        assert not utils._requested_tile_aligned_with_internal_tile(
            src_dst, bounds, 256, 256
        )

    with rasterio.open(COG_WEB_TILED) as src_dst:
        assert utils._requested_tile_aligned_with_internal_tile(
            src_dst, bounds, 256, 256
        ) 

def get_area_acres(geometry):
    """ Calculate area in acres for a GeoJSON geometry
    :param geometry: A GeoJSON Polygon geometry
    :returns: Area in acres
    """

    shapely_geometry = shape(geometry)
    geom_aea = transform(
        partial(
            pyproj.transform,
            pyproj.Proj(init="EPSG:4326"),
            pyproj.Proj(
                proj="aea",
                lat1=shapely_geometry.bounds[1],
                lat2=shapely_geometry.bounds[3],
            ),
        ),
        shapely_geometry,
    )
    return round(geom_aea.area / 4046.8564224) 

def test_tile_read_extmask():
    """Read masked area."""
    # non-boundless tile covering the masked part
    mercator_tile = mercantile.Tile(x=876431, y=1603669, z=22)
    bounds = mercantile.xy_bounds(mercator_tile)
    with rasterio.Env(GDAL_DISABLE_READDIR_ON_OPEN="TRUE"):
        with rasterio.open(S3_EXTMASK_PATH) as src_dst:
            arr, mask = reader.part(src_dst, bounds, 256, 256)
        assert arr.shape == (3, 256, 256)
        assert mask.shape == (256, 256)
        assert not mask.all()

    # boundless tile covering the masked part
    mercator_tile = mercantile.Tile(x=876431, y=1603668, z=22)
    bounds = mercantile.xy_bounds(mercator_tile)
    with rasterio.Env(GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR"):
        with rasterio.open(S3_MASK_PATH) as src_dst:
            arr, mask = reader.part(src_dst, bounds, 256, 256)
        assert arr.shape == (3, 256, 256)
        assert not mask.all() 

def histogram_match(self, use_bands, blm_source='browse', **kwargs):
        ''' Match the histogram to existing imagery '''
        warnings.warn('Histogram matching has changed due to the Maps API deprecation, see https://github.com/DigitalGlobe/gbdxtools/issues/778')
        assert has_rio, "To match image histograms please install rio_hist"
        data = self._read(self[use_bands,...], **kwargs)
        data = np.rollaxis(data.astype(np.float32), 0, 3)
        if 0 in data:
            data = np.ma.masked_values(data, 0)
        bounds = self._reproject(box(*self.bounds), from_proj=self.proj, to_proj="EPSG:4326").bounds
        ref = BrowseImage(self.cat_id, bbox=bounds).read()
        out = np.dstack([rio_match(data[:,:,idx], ref[:,:,idx].astype(np.double)/255.0)
                        for idx in range(data.shape[-1])])
        if 'stretch' in kwargs or 'gamma' in kwargs:
            return self._histogram_stretch(out, **kwargs)
        else:
            return out 

def test_resampling_returns_different_results():
    bounds = (
        -8844681.416934313,
        3757032.814272982,
        -8766409.899970293,
        3835304.331237001,
    )
    with rasterio.open(f"{LANDSAT_PATH}_B2.TIF") as src_dst:
        arr, _ = reader.part(
            src_dst, bounds, 16, 16, dst_crs=constants.WEB_MERCATOR_CRS
        )
        arr2, _ = reader.part(
            src_dst,
            bounds,
            16,
            16,
            dst_crs=constants.WEB_MERCATOR_CRS,
            resampling_method="bilinear",
        )

    assert not numpy.array_equal(arr, arr2) 

def test_make_tiles_tile_bounds(x, y):
    '''
    Test if children tiles from z10 are created correctly
    '''
    test_bounds = mercantile.bounds(x, y, 10)

    test_bbox = list(mercantile.xy(test_bounds.west, test_bounds.south)) + list(mercantile.xy(test_bounds.east, test_bounds.north))

    test_crs = 'epsg:3857'
    test_minz = 10
    test_maxz = 13

    created_tiles_gen = _make_tiles(test_bbox, test_crs, test_minz, test_maxz)

    assert isinstance(created_tiles_gen, types.GeneratorType)

    created_tiles = list(created_tiles_gen)

    assert len(created_tiles) == 85 

def create_pyramid_job(x, y, min_zoom, max_zoom, bounds):
    pyramid = {
        'tile': {
            'x': x,
            'y': y,
            'min_zoom': min_zoom,
            'max_zoom': max_zoom
        },
        'bounds': {
            'west': bounds.west,
            'south': bounds.south,
            'east': bounds.east,
            'north': bounds.north
        }
    }

    def payload_id():
        hash_obj = json.dumps(pyramid, sort_keys=True).encode('utf-8')
        return hashlib.sha1(hash_obj).hexdigest()

    return {
        'id': payload_id(),
        'type': 'pyramid',
        'pyramid': pyramid
    } 

def get_mvt(zoom,x,y):
    try:								# Sanitize the inputs
        sani_zoom,sani_x,sani_y = float(zoom),float(x),float(y)
        del zoom,x,y
    except:
        print('suspicious')
        return 1

    scale_denom = zoom_to_scale_denom(sani_zoom)
    tilebounds = bounds(sani_zoom,sani_x,sani_y)
    s,w,n,e = str(tilebounds['s']),str(tilebounds['w']),str(tilebounds['n']),str(tilebounds['e'])
    final_query = "EXECUTE gettile(!bbox!, !scale_denominator!, !pixel_width!, !pixel_height!);"
    sent_query = replace_tokens(final_query,s,w,n,e,scale_denom)
    response = list(session.execute(sent_query))
    print(sent_query)
    layers = filter(None,list(itertools.chain.from_iterable(response)))
    final_tile = b''
    for layer in layers:
        final_tile = final_tile + io.BytesIO(layer).getvalue() 
    return final_tile 

def test_src_meta_making(expectedMeta):
    bounds = merc.bounds(10, 10, 10)

    src_meta = untiler.make_src_meta(bounds, 4096)

    for k, e in zip(sorted(src_meta), sorted(expectedMeta)):
        assert k == e
        # assert src_meta[k] == expectedMeta[e]

    print("# OK - %s " % (inspect.stack()[0][3])) 

def make_src_meta(bounds, size, creation_opts={}):
    """
    Create metadata for output tiles
    """

    ul = merc.xy(bounds.west, bounds.north)
    lr = merc.xy(bounds.east, bounds.south)

    aff = make_affine(size, size, ul, lr)

    ## default values
    src_meta = {
        'driver': 'GTiff',
        'height': size,
        'width': size,
        'count': 4,
        'dtype': np.uint8,
        'affine': aff,
        "crs": 'EPSG:3857',
        'compress': 'JPEG',
        'tiled': True,
        'blockxsize': 256,
        'blockysize': 256
    }

    for c in creation_opts.keys():
        src_meta[c] = creation_opts[c]

    return src_meta 

def _calc_tms_zoom(self, scale):
        for z in range(15,20):
            b = mercantile.bounds(0,0,z)
            if scale > math.sqrt((b.north - b.south)*(b.east - b.west) / (256*256)):
                return z 

def get_tile_wms(tile, imagery, folder, kwargs):
    """
    Read a WMS endpoint with query parameters corresponding to a TMS tile

    Converts the tile boundaries to the spatial/coordinate reference system
    (SRS or CRS) specified by the WMS query parameter.
    """
    # retrieve the necessary parameters from the query string
    query_dict = parse_qs(imagery.lower())
    image_format = query_dict.get('format')[0].split('/')[1]
    wms_version = query_dict.get('version')[0]
    if wms_version == '1.3.0':
        wms_srs = query_dict.get('crs')[0]
    else:
        wms_srs = query_dict.get('srs')[0]

    # find our tile bounding box
    bound = bounds(*[int(t) for t in tile.split('-')])
    xmin, ymin, xmax, ymax = transform_bounds(WGS84_CRS, CRS.from_string(wms_srs), *bound, densify_pts=21)

    # project the tile bounding box from lat/lng to WMS SRS
    bbox = (
        [ymin, xmin, ymax, xmax] if wms_version == "1.3.0" else [xmin, ymin, xmax, ymax]
    )

    # request the image with the transformed bounding box and save
    wms_url = imagery.replace('{bbox}', ','.join([str(b) for b in bbox]))
    r = requests.get(wms_url, auth=kwargs.get('http_auth'))
    tile_img = op.join(folder, '{}.{}'.format(tile, image_format))
    with open(tile_img, 'wb') as w:
        w.write(r.content)
    return tile_img 

def getCorners(bounds, boolKey):
    corners = np.array([
        [bounds.west, bounds.south],
        [bounds.east, bounds.south],
        [bounds.east, bounds.north],
        [bounds.west, bounds.north]
        ])

    return [
        corners[coordOrd[boolKey][0]],
        corners[coordOrd[boolKey][1]]
    ] 

def createFacets(tileMin, tileMax, zoom, parentGet):
    for r in range(tileMin.y, tileMax.y + 1):
        for c in range(tileMin.x, tileMax.x + 1):
            quad = tools.quadtree(c, r, zoom)
            boolKey = (r+c) % 2 == 0
            n = parentGet.getParents('n', c, r, zoom)
            s = parentGet.getParents('s', c, r, zoom)
            coords = getCorners(mercantile.bounds(c, r, zoom), boolKey)
            nQT = ''.join(np.dstack((n, quad)).flatten()) + 'n'
            sQT = ''.join(np.dstack((s, quad)).flatten()) + 's'

            yield {
                "type": "Feature",
                "properties": {
                    "qt": nQT
                    },
                "geometry": {
                    "type": "Polygon",
                    "coordinates": [coords[0].tolist()]
                    }
                }
            yield {
                "type": "Feature",
                "properties": {
                    "qt": sQT,
                    },
                "geometry": {
                    "type": "Polygon",
                    "coordinates": [coords[1].tolist()]
                    }
                } 

def test_ovr_level():
    """Should return the correct overview level."""
    with rasterio.open(COG_DST) as src_dst:
        # raw/-1: 2667x2658 0: 1329x1334, 1: 665x667, 2: 333x334, 3: 167x167
        assert (
            utils.get_overview_level(
                src_dst, src_dst.bounds, 100, 100, dst_crs=src_dst.crs
            )
            == 3
        )
        assert (
            utils.get_overview_level(
                src_dst, src_dst.bounds, 200, 200, dst_crs=src_dst.crs
            )
            == 2
        )
        assert (
            utils.get_overview_level(
                src_dst, src_dst.bounds, 500, 500, dst_crs=src_dst.crs
            )
            == 1
        )
        assert (
            utils.get_overview_level(
                src_dst, src_dst.bounds, 800, 800, dst_crs=src_dst.crs
            )
            == 0
        )
        assert (
            utils.get_overview_level(
                src_dst, src_dst.bounds, 1500, 1500, dst_crs=src_dst.crs
            )
            == -1
        )
        assert (
            utils.get_overview_level(
                src_dst, src_dst.bounds, 3000, 3000, dst_crs=src_dst.crs
            )
            == -1
        ) 

def triangulate(zoom, output, bounds, tile, tableid):
    if bounds:
        bounds = np.array(bounds).astype(np.float64)
    elif tile:
        epsilon = 1.0e-10
        tile = np.array(tile).astype(np.uint16)
        tBounds = mercantile.bounds(*tile)
        bounds = np.array([
            tBounds.west + epsilon,
            tBounds.south + epsilon,
            tBounds.east - epsilon,
            tBounds.north - epsilon
            ])
    else:
        sys.exit('Error: A bounds or tile must be specified')

    tileMin = mercantile.tile(bounds[0], bounds[3], zoom)
    tileMax = mercantile.tile(bounds[2], bounds[1], zoom)

    pGet = facetParent()

    if tableid:
        gJSON = createDBinit(tileMin, tileMax, zoom, pGet, tableid)
    else:
        gJSON = createFacets(tileMin, tileMax, zoom, pGet)

    if output:
        with open(output, 'w') as oFile:
            for feat in gJSON:
                oFile.write(json.dumps(feat) + '\n')
    else:
        for feat in gJSON:
            click.echo(json.dumps(feat)) 

def _make_tiles(bbox, src_crs, minz, maxz):
    """
    Given a bounding box, zoom range, and source crs,
    find all tiles that would intersect

    Parameters
    -----------
    bbox: list
        [w, s, e, n] bounds
    src_crs: str
        the source crs of the input bbox
    minz: int
        minumum zoom to find tiles for
    maxz: int
        maximum zoom to find tiles for

    Returns
    --------
    tiles: generator
        generator of [x, y, z] tiles that intersect
        the provided bounding box
    """
    w, s, e, n = transform_bounds(*[src_crs, "epsg:4326"] + bbox, densify_pts=0)

    EPSILON = 1.0e-10

    w += EPSILON
    s += EPSILON
    e -= EPSILON
    n -= EPSILON

    for z in range(minz, maxz + 1):
        for x, y in _tile_range(mercantile.tile(w, n, z), mercantile.tile(e, s, z)):
            yield [x, y, z] 

def test_tile_read_crs():
    """Read tile using different target CRS and bounds CRS."""
    bounds = (
        -11663507.036777973,
        4715018.0897710975,
        -11663487.927520901,
        4715037.199028169,
    )
    tilesize = 16
    with rasterio.open(S3_PATH) as src_dst:
        # Test target CRS with input bounds in bounds_crs
        arr, mask = reader.part(
            src_dst,
            bounds,
            tilesize,
            tilesize,
            indexes=(3, 2, 1),
            dst_crs=constants.WGS84_CRS,
            bounds_crs=constants.WEB_MERCATOR_CRS,
        )
        assert arr.shape == (3, 16, 16)
        assert mask.shape == (16, 16)

        # Test target CRS with input bounds in target CRS
        bounds = (
            -104.7750663757324,
            38.95353532141203,
            -104.77489471435543,
            38.95366881479646,
        )
        arr_crs, _ = reader.part(
            src_dst,
            bounds,
            tilesize,
            tilesize,
            indexes=(3, 2, 1),
            dst_crs=constants.WGS84_CRS,
        )

        assert numpy.array_equal(arr, arr_crs) 

def test_tile_read_not_covering_the_whole_tile():
    """Should raise an error when dataset doesn't cover more than 50% of the tile."""
    bounds = (
        -9079495.967826376,
        3991847.365165044,
        -9001224.450862356,
        4070118.882129065,
    )
    tilesize = 16
    with pytest.raises(TileOutsideBounds):
        with rasterio.open(f"{LANDSAT_PATH}_B2.TIF") as src_dst:
            reader.part(src_dst, bounds, tilesize, tilesize, minimum_overlap=0.6)


# See https://github.com/cogeotiff/rio-tiler/issues/105#issuecomment-492268836 

def bounds(zoom,x,y):
    inProj = pyproj.Proj(init='epsg:4326')
    outProj = pyproj.Proj(init='epsg:3857')
    lnglatbbox = mercantile.bounds(x,y,zoom)
    ws = (pyproj.transform(inProj,outProj,lnglatbbox[0],lnglatbbox[1]))
    en = (pyproj.transform(inProj,outProj,lnglatbbox[2],lnglatbbox[3]))
    return {'w':ws[0],'s':ws[1],'e':en[0],'n':en[1]} 

def linear_rescale(
    image: numpy.ndarray,
    in_range: Tuple[Union[int, float], Union[int, float]] = (0, 1),
    out_range: Tuple[Union[int, float], Union[int, float]] = (1, 255),
) -> numpy.ndarray:
    """
    Linear rescaling.

    Attributes
    ----------
        image : numpy ndarray
            Image array to rescale.
        in_range : list, int, optional, (default: [0,1])
            Image min/max value to rescale.
        out_range : list, int, optional, (default: [1,255])
            output min/max bounds to rescale to.

    Returns
    -------
        out : numpy ndarray
            returns rescaled image array.

    """
    imin, imax = in_range
    omin, omax = out_range
    image = numpy.clip(image, imin, imax) - imin
    image = image / numpy.float(imax - imin)
    return image * (omax - omin) + omin 

def tile_exists(
    bounds: Tuple[float, float, float, float], tile_z: int, tile_x: int, tile_y: int
) -> bool:
    """
    Check if a mercatile tile is inside a given bounds.

    Attributes
    ----------
        bounds : list
            WGS84 bounds (left, bottom, right, top).
        z : int
            Mercator tile ZOOM level.
        y : int
            Mercator tile Y index.
        x : int
            Mercator tile Y index.

    Returns
    -------
        out : boolean
            if True, the z-x-y mercator tile in inside the bounds.

    """
    tile_bounds = mercantile.bounds(mercantile.Tile(tile_x, tile_y, tile_z))
    return (
        (tile_bounds[0] < bounds[2])
        and (tile_bounds[2] > bounds[0])
        and (tile_bounds[3] > bounds[1])
        and (tile_bounds[1] < bounds[3])
    ) 

def _requested_tile_aligned_with_internal_tile(
    src_dst: Union[DatasetReader, DatasetWriter, WarpedVRT],
    bounds: Tuple[float, float, float, float],
    height: int,
    width: int,
) -> bool:
    """Check if tile is aligned with internal tiles."""
    if not src_dst.is_tiled:
        return False

    if src_dst.crs != constants.WEB_MERCATOR_CRS:
        return False

    col_off, row_off, w, h = windows.from_bounds(
        *bounds, height=height, transform=src_dst.transform, width=width
    ).flatten()

    if round(w) % 64 and round(h) % 64:
        return False

    if (src_dst.width - round(col_off)) % 64:
        return False

    if (src_dst.height - round(row_off)) % 64:
        return False

    return True 

def test_tile_read_validResampling():
    """Should return a 1 band array and a mask."""
    bounds = (
        -8844681.416934313,
        3757032.814272982,
        -8766409.899970293,
        3835304.331237001,
    )
    with rasterio.open(f"{LANDSAT_PATH}_B2.TIF") as src_dst:
        arr, mask = reader.part(src_dst, bounds, 16, 16, resampling_method="bilinear")
    assert arr.shape == (1, 16, 16)
    assert mask.shape == (16, 16)