Python源码示例:affine.Affine()

示例1
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) 
示例2
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) 
示例3
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) 
示例4
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 
示例5
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()
    ) 
示例6
def _morpho(self, scount):
        aff = self._aff * affine.Affine.translation(-scount, -scount)
        return Footprint(
            gt=aff.to_gdal(),
            rsize=(self.rsize + 2 * scount),
        ) 
示例7
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) 
示例8
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 ******************************************************** ** 
示例9
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) 
示例10
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) 
示例11
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) 
示例12
def affine(self, new_affine):
        assert(isinstance(new_affine, Affine))
        self._affine = new_affine 
示例13
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 = [] 
示例14
def defaults(self):
        props = {
            'affine' : Affine(0,0,0,0,0,0),
            'shape' : (1,1),
            'nodata' : 0,
            'crs' : pyproj.Proj(_pyproj_init),
        }
        return props 
示例15
def nearest_cell(self, x, y, affine=None, snap='corner'):
        """
        Returns the index of the cell (column, row) closest
        to a given geographical coordinate.
 
        Parameters
        ----------
        x : int or float
            x coordinate.
        y : int or float
            y coordinate.
        affine : affine.Affine
                 Affine transformation that defines the translation between
                 geographic x/y coordinate and array row/column coordinate.
                 Defaults to self.affine.
        snap : str
               Indicates the cell indexing method. If "corner", will resolve to 
               snapping the (x,y) geometry to the index of the nearest top-left 
               cell corner. If "center", will return the index of the cell that 
               the geometry falls within.
        Returns
        -------
        x_i, y_i : tuple of ints
                   Column index and row index
        """
        if not affine:
            affine = self.affine
        try:
            assert isinstance(affine, Affine)
        except:
            raise TypeError('affine must be an Affine instance.')
        snap_dict = {'corner': np.around, 'center': np.floor}
        col, row = snap_dict[snap](~affine * (x, y)).astype(int)
        return col, row 
示例16
def affine(self, new_affine):
        assert isinstance(new_affine, Affine)
        self._affine = new_affine 
示例17
def get_geometry_mask(width: int, height: int,
                      geometry: GeometryLike,
                      x_min: float, y_min: float, res: float) -> np.ndarray:
    geometry = convert_geometry(geometry)
    # noinspection PyTypeChecker
    transform = affine.Affine(res, 0.0, x_min,
                              0.0, -res, y_min + res * height)
    return rasterio.features.geometry_mask([geometry],
                                           out_shape=(height, width),
                                           transform=transform,
                                           all_touched=True,
                                           invert=True) 
示例18
def make_transform(tilerange, zoom):
    ulx, uly = mercantile.xy(*mercantile.ul(tilerange['x']['min'], tilerange['y']['min'], zoom))
    lrx, lry = mercantile.xy(*mercantile.ul(tilerange['x']['max'], tilerange['y']['max'], zoom))
    xcell = (lrx - ulx) / float(tilerange['x']['max'] - tilerange['x']['min'])
    ycell = (uly - lry) / float(tilerange['y']['max'] - tilerange['y']['min'])
    return Affine(xcell, 0, ulx,
        0, -ycell, uly) 
示例19
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 
示例20
def test_get_area_def_from_raster(self):
        from pyresample import utils
        from rasterio.crs import CRS
        from affine import Affine
        x_size = 791
        y_size = 718
        transform = Affine(300.0379266750948, 0.0, 101985.0,
                           0.0, -300.041782729805, 2826915.0)
        crs = CRS(init='epsg:3857')
        if utils.is_pyproj2():
            # pyproj 2.0+ expands CRS parameters
            from pyproj import CRS
            proj_dict = CRS(3857).to_dict()
        else:
            proj_dict = crs.to_dict()
        source = tmptiff(x_size, y_size, transform, crs=crs)
        area_id = 'area_id'
        proj_id = 'proj_id'
        description = 'name'
        area_def = utils.rasterio.get_area_def_from_raster(
            source, area_id=area_id, name=description, proj_id=proj_id)
        self.assertEqual(area_def.area_id, area_id)
        self.assertEqual(area_def.proj_id, proj_id)
        self.assertEqual(area_def.description, description)
        self.assertEqual(area_def.width, x_size)
        self.assertEqual(area_def.height, y_size)
        self.assertDictEqual(proj_dict, area_def.proj_dict)
        self.assertTupleEqual(area_def.area_extent, (transform.c, transform.f + transform.e * y_size,
                                                     transform.c + transform.a * x_size, transform.f)) 
示例21
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) 
示例22
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) 
示例23
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) 
示例24
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 
示例25
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 
示例26
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 
示例27
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 
示例28
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") 
示例29
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)),
        )
    ) 
示例30
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)