Python affine.Affine() Examples

def ResampleRaster(InputRasterFile,OutputRasterFile,XResolution,YResolution=None,Format="ENVI"):

    """
    Description goes here...

    MDH

    """

    # import modules
    import rasterio, affine
    from rasterio.warp import reproject, Resampling

    # read the source raster
    with rasterio.open(InputRasterFile) as src:
        Array = src.read()
        OldResolution = src.res

        #setup output resolution
        if YResolution == None:
            YResolution = XResolution
        NewResolution = (XResolution,YResolution)


        # setup the transform to change the resolution
        XResRatio = OldResolution[0]/NewResolution[0]
        YResRatio = OldResolution[1]/NewResolution[1]
        NewArray = np.empty(shape=(Array.shape[0], int(round(Array.shape[1] * XResRatio)), int(round(Array.shape[2] * YResRatio))))
        Aff = src.affine
        NewAff = affine.Affine(Aff.a/XResRatio, Aff.b, Aff.c, Aff.d, Aff.e/YResRatio, Aff.f)

        # reproject the raster
        reproject(Array, NewArray, src_transform=Aff, dst_transform=NewAff, src_crs = src.crs, dst_crs = src.crs, resample=Resampling.bilinear)

        # write results to file
        with rasterio.open(OutputRasterFile, 'w', driver=src.driver, \
                            height=NewArray.shape[1],width=NewArray.shape[2], \
                            nodata=src.nodata,dtype=str(NewArray.dtype), \
                            count=src.count,crs=src.crs,transform=NewAff) as dst:
            dst.write(NewArray) 

def ConvertRaster2LatLong(InputRasterFile,OutputRasterFile):

    """
    Convert a raster to lat long WGS1984 EPSG:4326 coordinates for global plotting

    MDH

    """

    # import modules
    import rasterio
    from rasterio.warp import reproject, calculate_default_transform as cdt, Resampling

    # read the source raster
    with rasterio.open(InputRasterFile) as src:
        #get input coordinate system
        Input_CRS = src.crs
        # define the output coordinate system
        Output_CRS = {'init': "epsg:4326"}
        # set up the transform
        Affine, Width, Height = cdt(Input_CRS,Output_CRS,src.width,src.height,*src.bounds)
        kwargs = src.meta.copy()
        kwargs.update({
            'crs': Output_CRS,
            'transform': Affine,
            'affine': Affine,
            'width': Width,
            'height': Height
        })

        with rasterio.open(OutputRasterFile, 'w', **kwargs) as dst:
            for i in range(1, src.count+1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.affine,
                    src_crs=src.crs,
                    dst_transform=Affine,
                    dst_crs=Output_CRS,
                    resampling=Resampling.bilinear) 

def global_reader(value=None):
    array = np.ma.masked_array(
        np.empty(shape=(360, 720), dtype=np.float32, order="C"),
        np.empty(shape=(360, 720), dtype=np.bool, order="C"),
    )
    if value is None:
        array.mask.fill(True)
    else:
        array.fill(value)
        array.mask.fill(False)
    transform = Affine(-180.0, 0.5, 0.0, 90.0, 0.0, -0.5)
    reader = Mock(spec=DatasetReader)
    reader.read.return_value = array
    reader.shape = array.shape
    reader.transform = transform
    reader.bounds = (-180.0, -90.0, 180.0, 90.0)
    reader.window.return_value = ((0, 359), (0, 719))
    reader.window_transform.return_value = transform
    return reader 

def bounds_to_ranges(out_bounds=None, in_affine=None, in_shape=None):
    """
    Return bounds range values from geolocated input.

    Parameters
    ----------
    out_bounds : tuple
        left, bottom, right, top
    in_affine : Affine
        input geolocation
    in_shape : tuple
        input shape

    Returns
    -------
    minrow, maxrow, mincol, maxcol
    """
    return itertools.chain(
        *from_bounds(
            *out_bounds, transform=in_affine, height=in_shape[-2], width=in_shape[-1]
        ).round_lengths(pixel_precision=0).round_offsets(pixel_precision=0).toranges()
    ) 

def set_indices(self, new_indices):
        """
        Updates self.affine and self.shape to correspond to new indices representing
        a new bounding rectangle. Also resets self.mask.

        Parameters
        ----------
        new_indices : tuple of ints (length 4)
                      (xmin_index, ymin_index, xmax_index, ymax_index)
        """
        affine = self.affine
        assert all((isinstance(ix, int) for ix in new_indices))
        ul = np.asarray((new_indices[0], new_indices[3]))
        lr = np.asarray((new_indices[2], new_indices[1]))
        xmin, ymax = affine * tuple(ul)
        shape = tuple(lr - ul)[::-1]
        new_affine = Affine(affine.a, affine.b, xmin,
                            affine.d, affine.e, ymax)
        self.affine = new_affine
        self.shape = shape
        #TODO: For now, simply reset mask
        self.mask = np.ones(shape, dtype=np.bool) 

def set_bbox(self, new_bbox):
        """
        Sets new bbox while maintaining the same cell dimensions. Updates
        self.affine and self.shape. Also resets self.mask.

        Note that this method rounds the given bbox to match the existing
        cell dimensions.

        Parameters
        ----------
        new_bbox : tuple of floats (length 4)
                   (xmin, ymin, xmax, ymax)
        """
        affine = self.affine
        xmin, ymin, xmax, ymax = new_bbox
        ul = np.around(~affine * (xmin, ymax)).astype(int)
        lr = np.around(~affine * (xmax, ymin)).astype(int)
        xmin, ymax = affine * tuple(ul)
        shape = tuple(lr - ul)[::-1]
        new_affine = Affine(affine.a, affine.b, xmin,
                            affine.d, affine.e, ymax)
        self.affine = new_affine
        self.shape = shape
        #TODO: For now, simply reset mask
        self.mask = np.ones(shape, dtype=np.bool) 

def test_get_area_def_from_raster_rotated_value_err(self):
        from pyresample import utils
        from affine import Affine
        transform = Affine(300.0379266750948, 0.1, 101985.0,
                           0.0, -300.041782729805, 2826915.0)
        source = tmptiff(transform=transform)
        self.assertRaises(ValueError, utils.rasterio.get_area_def_from_raster, source) 

def size(self):
        """Spatial distances: (||raster left - raster right||, ||raster top - raster bottom||)"""
        return np.abs(~affine.Affine.rotation(self.angle) * self.diagvec, dtype=np.float64) 

def test_get_area_def_from_raster_non_georef_value_err(self):
        from pyresample import utils
        from affine import Affine
        transform = Affine(300.0379266750948, 0.0, 101985.0,
                           0.0, -300.041782729805, 2826915.0)
        source = tmptiff(transform=transform)
        self.assertRaises(ValueError, utils.rasterio.get_area_def_from_raster, source) 

def test_get_area_def_from_raster_non_georef_respects_proj_dict(self):
        from pyresample import utils
        from affine import Affine
        transform = Affine(300.0379266750948, 0.0, 101985.0,
                           0.0, -300.041782729805, 2826915.0)
        source = tmptiff(transform=transform)
        proj_dict = {'init': 'epsg:3857'}
        area_def = utils.rasterio.get_area_def_from_raster(source, proj_dict=proj_dict)
        if utils.is_pyproj2():
            from pyproj import CRS
            proj_dict = CRS(3857).to_dict()
        self.assertDictEqual(area_def.proj_dict, proj_dict) 

def raster_file(tmpdir_factory):
    import affine

    raster_data = np.arange(-128 * 256, 128 * 256, dtype='int16').reshape(256, 256)

    # Sprinkle in some more nodata
    raster_data.flat[::5] = 10000

    profile = {
        'driver': 'GTiff',
        'dtype': 'int16',
        'nodata': 10000,
        'width': raster_data.shape[1],
        'height': raster_data.shape[0],
        'count': 1,
        'crs': {'init': 'epsg:32637'},
        'transform': affine.Affine(
            2.0, 0.0, 694920.0,
            0.0, -2.0, 2055666.0
        )
    }

    outpath = tmpdir_factory.mktemp('raster')
    unoptimized_raster = outpath.join('img-raw.tif')
    with rasterio.open(str(unoptimized_raster), 'w', **profile) as dst:
        dst.write(raster_data, 1)

    optimized_raster = outpath.join('img.tif')
    cloud_optimize(unoptimized_raster, optimized_raster)

    return optimized_raster 

def big_raster_file_nodata(tmpdir_factory):
    import affine

    np.random.seed(17)
    raster_data = np.random.randint(0, np.iinfo(np.uint16).max, size=(2048, 2048), dtype='uint16')
    nodata = 10000

    # include some big nodata regions
    ix, iy = np.indices(raster_data.shape)
    circular_mask = np.sqrt((ix - raster_data.shape[0] / 2) ** 2
                            + (iy - raster_data.shape[1] / 2) ** 2) > 1000
    raster_data[circular_mask] = nodata
    raster_data[500:1000, 1000:2000] = nodata
    raster_data[1200, :] = nodata

    profile = {
        'driver': 'GTiff',
        'dtype': 'uint16',
        'nodata': nodata,
        'width': raster_data.shape[1],
        'height': raster_data.shape[0],
        'count': 1,
        'crs': {'init': 'epsg:32637'},
        'transform': affine.Affine(
            10.0, 0.0, 694920.0,
            0.0, -10.0, 2055666.0
        )
    }

    outpath = tmpdir_factory.mktemp('raster')
    unoptimized_raster = outpath.join('img-raw.tif')
    with rasterio.open(str(unoptimized_raster), 'w', **profile) as dst:
        dst.write(raster_data, 1)

    optimized_raster = outpath.join('img-nodata.tif')
    cloud_optimize(unoptimized_raster, optimized_raster)

    return optimized_raster 

def unoptimized_raster_file(tmpdir_factory):
    import affine

    np.random.seed(17)
    raster_data = np.random.randint(0, np.iinfo(np.uint16).max, size=(1024, 1024), dtype='uint16')
    nodata = 10000

    # include some big nodata regions
    ix, iy = np.indices(raster_data.shape)
    circular_mask = np.sqrt((ix - raster_data.shape[0] / 2) ** 2
                            + (iy - raster_data.shape[1] / 2) ** 2) > 400
    raster_data[circular_mask] = nodata
    raster_data[200:600, 400:800] = nodata
    raster_data[500, :] = nodata

    profile = {
        'driver': 'GTiff',
        'dtype': 'uint16',
        'nodata': nodata,
        'width': raster_data.shape[1],
        'height': raster_data.shape[0],
        'count': 1,
        'crs': {'init': 'epsg:32637'},
        'transform': affine.Affine(
            2.0, 0.0, 694920.0,
            0.0, -2.0, 2055666.0
        )
    }

    outpath = tmpdir_factory.mktemp('raster')
    unoptimized_raster = outpath.join('img-raw.tif')
    with rasterio.open(str(unoptimized_raster), 'w', **profile) as dst:
        dst.write(raster_data, 1)

    return unoptimized_raster 

def invalid_raster_file(tmpdir_factory):
    """A raster file that is all nodata"""
    import affine

    raster_data = np.full((256, 256), 0, dtype='uint16')
    profile = {
        'driver': 'GTiff',
        'dtype': 'uint16',
        'nodata': 0,
        'width': raster_data.shape[1],
        'height': raster_data.shape[0],
        'count': 1,
        'crs': {'init': 'epsg:32637'},
        'transform': affine.Affine(
            2.0, 0.0, 694920.0,
            0.0, -2.0, 2055666.0
        )
    }

    outpath = tmpdir_factory.mktemp('raster')
    unoptimized_raster = outpath.join('img-raw.tif')
    with rasterio.open(str(unoptimized_raster), 'w', **profile) as dst:
        dst.write(raster_data, 1)

    optimized_raster = outpath.join('img.tif')
    cloud_optimize(unoptimized_raster, optimized_raster)

    return optimized_raster 

def extract_from_array(in_raster=None, in_affine=None, out_tile=None):
    """
    Extract raster data window array.

    Parameters
    ----------
    in_raster : array or ReferencedRaster
    in_affine : ``Affine`` required if in_raster is an array
    out_tile : ``BufferedTile``

    Returns
    -------
    extracted array : array
    """
    if isinstance(in_raster, ReferencedRaster):
        in_affine, in_raster = in_raster.affine, in_raster.data

    # get range within array
    minrow, maxrow, mincol, maxcol = bounds_to_ranges(
        out_bounds=out_tile.bounds, in_affine=in_affine, in_shape=in_raster.shape
    )
    # if output window is within input window
    if (
        minrow >= 0 and
        mincol >= 0 and
        maxrow <= in_raster.shape[-2] and
        maxcol <= in_raster.shape[-1]
    ):
        return in_raster[..., minrow:maxrow, mincol:maxcol]
    # raise error if output is not fully within input
    else:
        raise ValueError("extraction fails if output shape is not within input") 

def tiles_to_affine_shape(tiles):
    """
    Return Affine and shape of combined tiles.

    Parameters
    ----------
    tiles : iterable
        an iterable containing BufferedTiles

    Returns
    -------
    Affine, Shape
    """
    if not tiles:
        raise TypeError("no tiles provided")
    pixel_size = tiles[0].pixel_x_size
    left, bottom, right, top = (
        min([t.left for t in tiles]),
        min([t.bottom for t in tiles]),
        max([t.right for t in tiles]),
        max([t.top for t in tiles]),
    )
    return (
        Affine(pixel_size, 0, left, 0, -pixel_size, top),
        Shape(
            width=int(round((right - left) / pixel_size, 0)),
            height=int(round((top - bottom) / pixel_size, 0)),
        )
    ) 

def test_transform_from_latlon(lon_start, dlon, lat_start, dlat):

    lon = np.arange(lon_start, 20, dlon)
    lat = np.arange(lat_start, 20, dlat)

    r = _transform_from_latlon(lon, lat)

    assert isinstance(r, Affine)

    expected = np.array(
        [dlon, 0, lon_start - dlon / 2, 0, dlat, lat_start - dlat / 2, 0, 0, 1]
    )

    assert np.allclose(np.array(r), expected) 

def _make_affine(fr_shape, to_shape):
    """Given from and to width and height,
    compute affine transform defining the
    georeferencing of the output array
    """
    fr_window_affine = Affine(
        1, 0, 0,
        0, -1, 0)

    to_window_affine = Affine(
        (fr_shape[1] / float(to_shape[1])), 0, 0,
        0, -(fr_shape[0] / float(to_shape[0])), 0)

    return fr_window_affine, to_window_affine 

def test_rasterio_order_list(self):
        truth_affine_obj = Affine(0.5, 0.0, 733601.0, 0.0, -0.5, 3725139.0)
        affine_list = [0.5, 0.0, 733601.0, 0.0, -0.5, 3725139.0]
        test_affine = list_to_affine(affine_list)

        assert truth_affine_obj == test_affine 

def test_gdal_order_list(self):
        truth_affine_obj = Affine(0.5, 0.0, 733601.0, 0.0, -0.5, 3725139.0)
        gdal_affine_list = [733601.0, 0.5, 0.0, 3725139.0, 0.0, -0.5]
        test_affine = list_to_affine(gdal_affine_list)

        assert truth_affine_obj == test_affine 

def test_get_from_file(self):
        affine_obj = get_geo_transform(os.path.join(data_dir,
                                                    'sample_geotiff.tif'))
        assert affine_obj == Affine(0.5, 0.0, 733601.0, 0.0, -0.5, 3725139.0) 

def test_get_from_opened_raster(self):
        src_obj = sample_load_rasterio()
        affine_obj = get_geo_transform(src_obj)
        assert affine_obj == Affine(0.5, 0.0, 733601.0, 0.0, -0.5, 3725139.0)
        src_obj.close() 

def test_get_from_gdal(self):
        src_obj = sample_load_gdal()
        affine_obj = get_geo_transform(src_obj)
        assert affine_obj == Affine(0.5, 0.0, 733601.0, 0.0, -0.5, 3725139.0) 

def reader_context(**attrs):
    reader = Mock(
        spec=DatasetReader,
        shape=(360, 720),
        transform=Affine(-180., 0.5, 0.0, -90., 0.0, 0.5),
        bounds=BoundingBox(-180., -90., 0., 0.),
        crs={"init": "epsg:4326"},
    )
    reader.configure_mock(**attrs)
    context = Mock()
    context.__enter__ = Mock(return_value=reader)
    context.__exit__ = Mock(return_value=False)
    return context 

def _morpho(self, scount):
        aff = self._aff * affine.Affine.translation(-scount, -scount)
        return Footprint(
            gt=aff.to_gdal(),
            rsize=(self.rsize + 2 * scount),
        ) 

def rlength(self):
        """Pixel quantity: pixel count in the outer ring"""
        rx, ry = self.rsize
        # Convert to int before multiplication to avoid overflow
        inner_area = max(0, int(rx) - 2) * max(0, int(ry) - 2)
        return self.rarea - inner_area

    # Accessors - Affine transformations ******************************************************** ** 

def scale(self):
        """Spatial vector: scale used in the affine transformation, np.abs(scale) == pxsize"""
        aff = ~affine.Affine.rotation(self.angle)
        tl = np.asarray(aff * self.tl)
        br = np.asarray(aff * self.br)
        return np.asarray((br - tl) / self.rsize, dtype=np.float64) 

def __init__(self, affine, shape, mask=None, nodata=None,
                 crs=pyproj.Proj(_pyproj_init),
                 y_coord_ix=0, x_coord_ix=1):
        if affine is not None:
            self.affine = affine
        else:
            self.affine = Affine(0,0,0,0,0,0)
        super().__init__(shape=shape, mask=mask, nodata=nodata, crs=crs,
                         y_coord_ix=y_coord_ix, x_coord_ix=x_coord_ix) 

def affine(self, new_affine):
        assert(isinstance(new_affine, Affine))
        self._affine = new_affine 

def __init__(self, affine=Affine(0,0,0,0,0,0), shape=(1,1), nodata=0,
                 crs=pyproj.Proj(_pyproj_init),
                 mask=None):
        self.affine = affine
        self.shape = shape
        self.nodata = nodata
        self.crs = crs
        # TODO: Mask should be a raster, not an array
        if mask is None:
            self.mask = np.ones(shape)
        self.grids = []