Python rasterio.warp.transform() Examples

def feature_to_mercator(feature):
    """Normalize feature and converts coords to 3857.

    Args:
      feature: geojson feature to convert to mercator geometry.
    """
    # Ref: https://gist.github.com/dnomadb/5cbc116aacc352c7126e779c29ab7abe

    src_crs = CRS.from_epsg(4326)
    dst_crs = CRS.from_epsg(3857)

    geometry = feature["geometry"]
    if geometry["type"] == "Polygon":
        xys = (zip(*part) for part in geometry["coordinates"])
        xys = (list(zip(*transform(src_crs, dst_crs, *xy))) for xy in xys)

        yield {"coordinates": list(xys), "type": "Polygon"}

    elif geometry["type"] == "MultiPolygon":
        for component in geometry["coordinates"]:
            xys = (zip(*part) for part in component)
            xys = (list(zip(*transform(src_crs, dst_crs, *xy))) for xy in xys)

            yield {"coordinates": list(xys), "type": "Polygon"} 

def burn(tile, features, size):
    """Burn tile with features.

    Args:
      tile: the mercantile tile to burn.
      features: the geojson features to burn.
      size: the size of burned image.

    Returns:
      image: rasterized file of size with features burned.
    """

    # the value you want in the output raster where a shape exists
    burnval = 1
    shapes = ((geometry, burnval) for feature in features for geometry in feature_to_mercator(feature))

    bounds = mercantile.xy_bounds(tile)
    transform = from_bounds(*bounds, size, size)

    return rasterize(shapes, out_shape=(size, size), transform=transform) 

def transform_coordinates(source_srs, target_srs, x, y):
    return tuple(i[0] for i in transform(source_srs, target_srs, x, y)) 

def get_tile_tif(tile, imagery, folder, kwargs):
    """
    Read a GeoTIFF with a window corresponding to a TMS tile

    The TMS tile bounds are converted to the GeoTIFF source CRS. That bounding
    box is converted to a pixel window which is read from the GeoTIFF. For
    remote files which are internally tiled, this will take advantage of HTTP
    GET Range Requests to avoid downloading the entire file. See more info at:
    http://www.cogeo.org/in-depth.html
    """
    bound = bounds(*[int(t) for t in tile.split('-')])
    imagery_offset = kwargs.get('imagery_offset') or [0, 0]
    with rasterio.open(imagery) as src:
        x_res, y_res = src.transform[0], src.transform[4]

        # offset our imagery in the "destination pixel" space
        offset_bound = dict()
        deg_per_pix_x = (bound.west - bound.east) / 256
        deg_per_pix_y = (bound.north - bound.south) / 256

        offset_bound['west'] = bound.west + imagery_offset[0] * deg_per_pix_x
        offset_bound['east'] = bound.east + imagery_offset[0] * deg_per_pix_x
        offset_bound['north'] = bound.north + imagery_offset[1] * deg_per_pix_y
        offset_bound['south'] = bound.south + imagery_offset[1] * deg_per_pix_y

        # project tile boundaries from lat/lng to source CRS
        x, y = transform(WGS84_CRS, src.crs, [offset_bound['west']], [offset_bound['north']])
        tile_ul_proj = x[0], y[0]

        x, y = transform(WGS84_CRS, src.crs, [offset_bound['east']], [offset_bound['south']])
        tile_lr_proj = x[0], y[0]

        # get origin point from the TIF
        tif_ul_proj = (src.bounds.left, src.bounds.top)

        # use the above information to calculate the pixel indices of the window
        top = int((tile_ul_proj[1] - tif_ul_proj[1]) / y_res)
        left = int((tile_ul_proj[0] - tif_ul_proj[0]) / x_res)
        bottom = int((tile_lr_proj[1] - tif_ul_proj[1]) / y_res)
        right = int((tile_lr_proj[0] - tif_ul_proj[0]) / x_res)

        window = ((top, bottom), (left, right))

        # read the first three bands (assumed RGB) of the TIF into an array
        data = np.empty(shape=(3, 256, 256)).astype(src.profile['dtype'])
        for k in (1, 2, 3):
            src.read(k, window=window, out=data[k - 1], boundless=True)

        # save
        tile_img = op.join(folder, '{}{}'.format(tile, '.jpg'))
        img = Image.fromarray(np.moveaxis(data, 0, -1), mode='RGB')
        img.save(tile_img)

    return tile_img