Python rasterio.Affine() Examples

def test_data():
    test_data_pan = np.array([
        (np.random.rand(60, 60) * 255).astype(np.uint8)
        ])
    test_data_pan = test_data_pan[0]
    test_data_rgb = np.array([
        (np.random.rand(30, 30) * 255).astype(np.uint8)
        for i in range(3)
        ])
    test_data_src_aff = rasterio.Affine(2.0, 0.0, 0.0, 0.0, - 2.0, 0.0)
    test_data_src_crs = {'init': 'EPSG:3857'}
    test_data_dst_aff = rasterio.Affine(1.0, 0.0, 0.0, 0.0, - 1.0, 0.0)
    test_data_dst_crs = {'init': 'EPSG:3857'}

    return test_data_pan, test_data_rgb, test_data_src_aff,\
        test_data_src_crs, test_data_dst_aff, test_data_dst_crs 

def test_pansharp_data():
    b8_path = 'tests/fixtures/tiny_20_tiffs/LC81070352015122LGN00/'\
              'LC81070352015122LGN00_B8.tif'
    b4_path = 'tests/fixtures/tiny_20_tiffs/LC81070352015122LGN00/'\
              'LC81070352015122LGN00_B4.tif'
    b3_path = 'tests/fixtures/tiny_20_tiffs/LC81070352015122LGN00/'\
              'LC81070352015122LGN00_B3.tif'
    b2_path = 'tests/fixtures/tiny_20_tiffs/LC81070352015122LGN00/'\
              'LC81070352015122LGN00_B2.tif'
    band_paths = [b8_path, b4_path, b3_path, b2_path]
    pan_window = ((1536, 1792), (1280, 1536))
    g_args = {'half_window': False,
              'dst_aff': Affine(75.00483870967741, 0.0, 300892.5,
                                0.0, -75.00475285171103, 4107007.5),
              'verb': False, 'weight': 0.2,
              'dst_crs': {'init': u'epsg:32654'},
              'r_crs': {'init': u'epsg:32654'},
              'dst_dtype': np.__dict__['uint16'],
              'r_aff': Affine(150.0193548387097, 0.0, 300885.0,
                              0.0, -150.0190114068441, 4107015.0),
              'src_nodata': 0}

    return [rasterio.open(f) for f in band_paths],\
        pan_window, (6, 5), g_args 

def merge_rasters(folder, percentile=70):
    """Merge a set of monthly rasters keeping the nth percentile value.

    Used to remove transient features from time-series data.

    Parameters
    ----------
    folder : str, Path
        Folder containing rasters to be merged.
    percentile : int, optional (default 70.)
        Percentile value to use when merging using np.nanpercentile.
        Lower values will result in lower values/brightness.

    Returns
    -------
    raster_merged : numpy array
        The merged array.
    affine : affine.Affine
        The affine transformation for the merged raster.
    """

    affine = None
    rasters = []

    for file in os.listdir(folder):
        if file.endswith(".tif"):
            ntl_rd = rasterio.open(os.path.join(folder, file))
            rasters.append(ntl_rd.read(1))

            if not affine:
                affine = ntl_rd.transform

    raster_arr = np.array(rasters)

    raster_merged = np.percentile(raster_arr, percentile, axis=0)

    return raster_merged, affine 

def makeTesting(output, size, windowsize, bands):
    """Construct test fixture"""
    kwargs = {
        "count": bands,
        "crs": {"init": u"epsg:3857"},
        "dtype": "uint8",
        "driver": u"GTiff",
        "transform": Affine(
            4.595839562240513,
            0.0,
            -13550756.3744,
            0.0,
            -4.595839562240513,
            6315533.02503,
        ),
        "height": size,
        "width": size,
        "compress": "lzw",
        "blockxsize": windowsize,
        "blockysize": windowsize,
        "tiled": True,
    }

    randArr = np.array(
        [(np.random.rand(size, size) * 255).astype(np.uint8) for i in range(bands)]
    )

    with rasterio.open(output, "w", **kwargs) as dst:
        dst.write(randArr) 

def union(inputtiles, parsenames):

    tiles = sutils.tile_parser(inputtiles, parsenames)

    xmin, xmax, ymin, ymax = sutils.get_range(tiles)

    zoom = sutils.get_zoom(tiles)

    # make an array of shape (xrange + 3, yrange + 3)
    burn = sutils.burnXYZs(tiles, xmin, xmax, ymin, ymax, 0)

    nw = mercantile.xy(*mercantile.ul(xmin, ymin, zoom))

    se = mercantile.xy(*mercantile.ul(xmax + 1, ymax + 1, zoom))

    aff = Affine(((se[0] - nw[0]) / float(xmax - xmin + 1)), 0.0, nw[0],
        0.0, -((nw[1] - se[1]) / float(ymax - ymin + 1)), nw[1])

    unprojecter = sutils.Unprojecter()

    unionedTiles = [
        {
            'geometry': unprojecter.unproject(feature),
            'properties': {},
            'type': 'Feature'
        } for feature, shapes in features.shapes(np.asarray(np.flipud(np.rot90(burn)).astype(np.uint8), order='C'), transform=aff) if shapes == 1
    ]

    return unionedTiles 

def updateBoundsAffine(inAffine):
    from rasterio import Affine, coords

    bounds = coords.BoundingBox(
        inAffine.c - 180.0 + (inAffine.a / 2.0),
        -inAffine.f,
        -(inAffine.c - 180.0 + (inAffine.a / 2.0)),
        inAffine.f)

    outAffine = Affine(inAffine.a / 2.0, inAffine.b,inAffine.c - 180.0 + (inAffine.a / 2.0),
         inAffine.d,inAffine.e / 2.0, inAffine.f)

    return outAffine 

def make_affine(height, width, ul, lr):
    """
    Create an affine for a tile of a given size
    """
    xCell = (ul[0] - lr[0]) / width
    yCell = (ul[1] - lr[1]) / height
    return Affine(-xCell, 0.0, ul[0],
        0.0, -yCell, ul[1]) 

def affaux(up):
    return Affine(1, 0, 0, 0, -1, 0), Affine(up, 0, 0, 0, -up, 0) 

def resampleAffine(otrans, factor):
    from rasterio import Affine
    return Affine(otrans.a / float(factor),otrans.b,otrans.c,
             otrans.d,otrans.e / float(factor), otrans.f) 

def test_file_writer():
    test_data = np.zeros((3, 256, 256), dtype=np.uint8)

    test_opts = {
        'driver': 'PNG',
        'dtype': 'uint8',
        'height': 512,
        'width': 512,
        'count': 3,
        'crs': 'EPSG:3857'
    }

    test_affine = Affine(1, 0, 0, 0, -1, 0)

    test_bytearray = _encode_as_png(test_data, test_opts, test_affine)

    assert len(test_bytearray) == 842

    test_complex_data = test_data.copy()

    test_complex_data[0] += (np.random.rand(256, 256) * 255).astype(np.uint8)
    test_complex_data[1] += 10

    test_bytearray_complex = _encode_as_png(test_complex_data, test_opts, test_affine)

    assert len(test_bytearray) < len(test_bytearray_complex) 

def export_prediction_to_tif(self, out_file_path, prediction):
        """

        :param out_file_path:
        :param prediction: a np-array where second dim indicate n-channels
        :return:
        """

        #Check if map-info is available
        if self.map_info is None:
            class MissingMapInfoException(Exception): pass
            raise MissingMapInfoException()

        #Compute geo-meta data
        geo = self.map_info['transform']
        transf = rasterio.Affine(geo[1], geo[2], geo[0], geo[4], geo[5], geo[3])
        crs = {'init': self.map_info['cs_code']}

        #Write to file
        with rasterio.open(out_file_path, "w", driver="GTiff", compress="lzw", bigtiff="YES",
                           height=prediction.shape[0], width=prediction.shape[1], count=prediction.shape[2],
                           dtype=prediction.dtype,
                           crs=crs, transform=transf) as out_file:

            for band_no in range(0, prediction.shape[2]):
                out_file.write(prediction[:,:,band_no], band_no + 1)
        print('Exported predictions to', out_file_path) 

def test_reproject():
    from rasterio.warp import reproject
    from rasterio.enums import Resampling

    with rasterio.Env():
        # As source: a 1024 x 1024 raster centered on 0 degrees E and 0
        # degrees N, each pixel covering 15".
        rows, cols = src_shape = (1024, 1024)
        # decimal degrees per pixel
        d = 1.0 / 240

        # The following is equivalent to
        # A(d, 0, -cols*d/2, 0, -d, rows*d/2).
        src_transform = rasterio.Affine.translation(
                    -cols*d/2,
                    rows*d/2) * rasterio.Affine.scale(d, -d)
        src_crs = {'init': 'EPSG:4326'}
        source = np.ones(src_shape, np.uint8) * 255

        # Destination: a 2048 x 2048 dataset in Web Mercator (EPSG:3857)
        # with origin at 0.0, 0.0.
        dst_shape = (2048, 2048)
        dst_transform = Affine.from_gdal(
            -237481.5, 425.0, 0.0, 237536.4, 0.0, -425.0)
        dst_crs = {'init': 'EPSG:3857'}
        destination = np.zeros(dst_shape, np.uint8)

        reproject(
            source,
            destination,
            src_transform=src_transform,
            src_crs=src_crs,
            dst_transform=dst_transform,
            dst_crs=dst_crs,
            resampling=Resampling.nearest)

        # Assert that the destination is only partly filled.
        assert destination.any()
        assert not destination.all()


# Testing upsample function