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)