Source code for buzzard._footprint
""">>> help(Footprint)"""
# pylint: disable=too-many-lines
from __future__ import division, print_function
import logging
import itertools
import shapely
import shapely.geometry as sg
import affine
import numpy as np
import scipy
from osgeo import gdal
from osgeo import ogr
from osgeo import osr
from six.moves import filterfalse
import scipy.ndimage as ndi
import skimage.morphology as skm
from buzzard import _tools
from buzzard._tools import conv
from buzzard._tools import GDALErrorCatcher as Catch
from buzzard._env import env
from buzzard._footprint_tile import TileMixin
from buzzard._footprint_intersection import IntersectionMixin
from buzzard._footprint_move import MoveMixin
LOGGER = logging.getLogger('buzzard')
[docs]class Footprint(TileMixin, IntersectionMixin, MoveMixin):
"""Immutable object representing the location and size of a spatially localized raster. All
methods are thread-safe.
The :code:`Footprint` class:
- is a toolbox class designed to position a rectangle in both image space and geometry space,
- its main purpose is to simplify the manipulation of windows in rasters,
- has many accessors,
- has many algorithms,
- is a constant object,
- is designed to work with any rectangle in space (like non north-up/west-left rasters),
- is independent from projections, units and files,
- uses :code:`affine` library internally for conversions (https://github.com/sgillies/affine).
.. warning::
This class being complex and full python, the constructor is too slow for certain use cases (~0.5ms).
+-------------------------------------------------+--------------------------------------------------------+
| Method category | Method names |
+==============+==================================+========================================================+
| Footprint | from scratch | __init__, of_extent |
| construction +----------------------------------+--------------------------------------------------------+
| | from :code:`Footprint` | __and__, intersection, erode, dilate, ... |
+--------------+----------------------------------+--------------------------------------------------------+
| Conversion | extent, coords, geom, __geo_interface__ |
+--------------+----------------------------------+--------------------------------------------------------+
| Accessors | Spatial - Size and vectors | size, width, height, diagvec, ... |
| +----------------------------------+--------------------------------------------------------+
| | Spatial - Coordinates | tl, bl, br, tr, ... |
| +----------------------------------+--------------------------------------------------------+
| | Spatial - Misc | area, length, semiminoraxis, ...x |
| +----------------------------------+--------------------------------------------------------+
| | Raster - Size | rsize, rwidth, rheight, ... |
| +----------------------------------+--------------------------------------------------------+
| | Raster - Indices | rtl, rbl, rbr, ttr, ... |
| +----------------------------------+--------------------------------------------------------+
| | Raster - Misc | rarea, rlength, rsemiminoraxis, ... |
| +----------------------------------+--------------------------------------------------------+
| | Affine transformations | pxsize, pxvec, angle, ... |
+--------------+----------------------------------+--------------------------------------------------------+
| Binary predicates | __eq__, ... |
+-------------------------------------------------+--------------------------------------------------------+
| Numpy | shape, meshgrid_raster, meshgrid_spatial, slice_in, ...|
+-------------------------------------------------+--------------------------------------------------------+
| Coordinates conversions | spatial_to_raster, raster_to_spatial |
+-------------------------------------------------+--------------------------------------------------------+
| Geometry / Raster conversions | find_polygons, burn_polygons, ... |
+-------------------------------------------------+--------------------------------------------------------+
| Tiling | tile, tile_count, tile_occurrence |
+-------------------------------------------------+--------------------------------------------------------+
| Serialization | __str__, ... |
+-------------------------------------------------+--------------------------------------------------------+
Informations on geo transforms (gt) and affine matrices
-------------------------------------------------------
- http://www.perrygeo.com/python-affine-transforms.html
- https://pypi.python.org/pypi/affine/1.0
GDAL ordering:
+-----+------------------+--------------+-----+-----------------+-------------------+
| c | a | b | f | d | e |
+-----+------------------+--------------+-----+-----------------+-------------------+
| tlx | width of a pixel | row rotation | tly | column rotation | height of a pixel |
+-----+------------------+--------------+-----+-----------------+-------------------+
>>> c, a, b, f, d, e = fp.gt
>>> tlx, dx, rx, tly, ry, dy = fp.gt
Matrix ordering:
+------------------+--------------+-----+-----------------+-------------------+-----+
| a | b | c | d | e | f |
+------------------+--------------+-----+-----------------+-------------------+-----+
| width of a pixel | row rotation | tlx | column rotation | height of a pixel | tly |
+------------------+--------------+-----+-----------------+-------------------+-----+
>>> a, b, c, d, e, f = fp.aff6
>>> dx, rx, tlx, ry, dy, tly = fp.aff6
"""
__slots__ = ['_tl', '_bl', '_br', '_tr', '_aff', '_rsize', '_significant_min']
# Footprint construction ******************************************************************** **
# Footprint construction - from scratch ***************************************************** **
def __init__(self, **kwargs):
"""There are only two ways to construct a Footprint, but several high level constructors
exist, such as `.intersection`.
Usage 1
>>> buzz.Footprint(tl=(0, 10), size=(10, 10), rsize=(100, 100))
Usage 2
>>> buzz.Footprint(gt=(0, .1, 0, 10, 0, -.1), rsize=(100, 100))
Parameters
----------
tl: (nbr, nbr)
raster spatial top left coordinates
gt: (nbr, nbr, nbr, nbr, nbr, nbr)
geotransforms with GDAL ordering
size: (nbr, nbr)
Size of Footprint in space (unsigned)
rsize: (int, int)
Size of raster in pixel (unsigned integers)
"""
if 'rsize' not in kwargs:
raise ValueError('Missing `rsize` parameter')
rsize = np.asarray(kwargs.pop('rsize'), dtype='int32')
if rsize.shape != (2,):
raise ValueError('Invalid rsize shape `%s`' % str(rsize.shape))
if not np.isfinite(rsize).all() or (rsize <= 0).any():
raise ValueError('Invalid rsize value `%s`' % rsize)
if 'gt' in kwargs:
gt = np.asarray(kwargs.pop('gt'), dtype='float64')
if gt.shape != (6,):
raise ValueError('Invalid gt shape `%s`' % gt.shape)
if not np.isfinite(gt).all():
raise ValueError('Invalid gt value `%s`' % gt)
c, a, b, f, d, e = gt
elif 'tl' in kwargs and 'size' in kwargs:
tl = np.asarray(kwargs.pop('tl'), dtype='float64')
if tl.shape != (2,):
raise ValueError('Invalid tl shape `%s`' % tl.shape)
if not np.isfinite(tl).all():
raise ValueError('Invalid tl value `%s`' % tl)
b, d = 0., 0.
c, f = tl
size = np.asarray(kwargs.pop('size'), dtype='float64')
if size.shape != (2,):
raise ValueError('Invalid size shape `%s`' % size.shape)
if not np.isfinite(size).all() or (size <= 0).any():
raise ValueError('Invalid size value `%s`' % size)
a = size[0] / float(rsize[0])
e = -size[1] / float(rsize[1])
else:
raise ValueError('Provide `size & gt` or `rsize & size & tl`')
if kwargs:
raise ValueError('Unknown parameters [{}]'.format(kwargs.keys()))
if a * e - d * b == 0:
raise ValueError('Determinent should not be 0: {}'.format(
a * e - d * b
))
if b != 0 or d != 0 or a <= 0 or e >= 0:
if not env.allow_complex_footprint:
arr = np.asarray([[a, b, c], [d, e, f]])
arr = np.array2string(arr, precision=17).replace('\n', ' ')
s = ('Creating a non north-up/west-left Footprint. Buzzard is designed to handle '
'those Footprints too, but no unit tests are there yet. '
'Use buzz.Env(allow_complex_footprint=True) in a `with statement` to '
'deactivate this error.affine matrix:\n{}').format(arr)
raise ValueError(s)
rsizex, rsizey = rsize
aff = affine.Affine(a, b, c, d, e, f)
self._tl = c, f
self._bl = aff * (0, rsizey)
self._br = aff * (rsizex, rsizey)
self._tr = aff * (rsizex, 0)
self._aff = aff
self._rsize = rsize
self._tl = np.asarray(self._tl, dtype=np.float64)
self._bl = np.asarray(self._bl, dtype=np.float64)
self._br = np.asarray(self._br, dtype=np.float64)
self._tr = np.asarray(self._tr, dtype=np.float64)
self._rsize = np.asarray(self._rsize, dtype=env.default_index_dtype)
rect = _tools.Rect(*self.coords)
self._significant_min = rect.significant_min((rect.size / self._rsize).min())
# Footprint construction - from Footprint *************************************************** **
[docs] def __and__(self, other):
"""Returns Footprint.intersection"""
return self.intersection(other)
[docs] @classmethod
def of_extent(cls, extent, scale):
"""Create a Footprint from a rectangle extent and a scale
Parameters
----------
extent: (nbr, nbr, nbr, nbr)
Spatial coordinates of (minx, maxx, miny, maxy) defining a rectangle
scale: nbr or (nbr, nbr)
Resolution of output Footprint:
* if nbr: resolution = [a, -a]
* if (nbr, nbr): resolution [a, b]
"""
# Check extent parameter
extent = np.asarray(extent, dtype='float64')
if extent.shape != (4,):
raise ValueError('Invalid extent shape `{}`'.format(extent.shape))
if not np.isfinite(extent).all():
raise ValueError('Invalid extent value `{}`'.format(extent))
if extent[0] == extent[1] or extent[2] == extent[3]:
raise ValueError('Empty extent')
# Check scale parameter
scale = np.asarray(scale, dtype='float64')
if scale.ndim == 0:
scale = np.asarray([scale, -scale], dtype='float64')
elif scale.shape == (1,):
scale = np.asarray([scale[0], -scale[0]], dtype='float64')
elif scale.shape != (2,):
raise ValueError('scale has shape {} instead of (2,)'.format(scale.shape))
if (scale == 0).any():
raise ValueError('scale should be greater than 0')
# Work
minx, maxx, miny, maxy = extent
rect = _tools.Rect(
tl=(minx, maxy),
bl=(minx, miny),
br=(maxx, miny),
tr=(maxx, maxy),
)
pxsize = np.abs(scale)
significant_min = rect.significant_min(pxsize.min())
if env.significant <= significant_min:
s = ('This Footprint have large coordinates and small pixels, at least {:.2} '
'significant digits are necessary to perform this operation, but '
'`buzz.env.significant` is set to {}. Increase this value by using '
'buzz.Env(significant={}) in a `with statement`.'
).format(significant_min, env.significant, env.significant + 1)
raise RuntimeError(s)
abstract_grid_density = rect.abstract_grid_density(pxsize.min())
rsize = np.around(rect.size / pxsize * abstract_grid_density, 0) / abstract_grid_density
size = rsize * pxsize
return cls(tl=rect.tl, size=size, rsize=rsize)
[docs] def clip(self, startx, starty, endx, endy):
"""Construct a new Footprint by clipping self using pixel indices
To clip using coordinates see `Footprint.intersection`.
Parameters
----------
startx: int or None
Same rules as regular python slicing
starty: int or None
Same rules as regular python slicing
endx: int or None
Same rules as regular python slicing
endy: int or None
Same rules as regular python slicing
Returns
-------
fp: Footprint
The new clipped :code:`Footprint`
"""
startx, endx, _ = slice(startx, endx).indices(self.rsizex)
starty, endy, _ = slice(starty, endy).indices(self.rsizey)
rsize = np.asarray(
[endx - startx, endy - starty]
)
tl = self.tl + startx * self.pxlrvec + starty * self.pxtbvec
gt = self.gt
gt[0] = tl[0]
gt[3] = tl[1]
return self.__class__(
gt=gt,
rsize=rsize,
)
def _morpho(self, scount):
aff = self._aff * affine.Affine.translation(-scount, -scount)
return Footprint(
gt=aff.to_gdal(),
rsize=(self.rsize + 2 * scount),
)
[docs] def erode(self, count):
"""Construct a new Footprint from self, eroding all edges by :code:`count` pixels"""
assert count >= 0
assert count == int(count)
return self._morpho(-count)
[docs] def dilate(self, count):
"""Construct a new Footprint from self, dilating all edges by :code:`count` pixels"""
assert count >= 0
assert count == int(count)
return self._morpho(count)
[docs] def intersection(self, *objects, **kwargs):
"""intersection(self, *objects, scale='self', rotation='auto', alignment='auto', \
homogeneous=False)
Construct a Footprint bounding the intersection of geometric objects, self being one of the
of input geometry. Inputs' intersection is always within output Footprint.
Parameters
----------
*objects: *object
Any object with a __geo_interface__ attribute defining a geometry, like a Footprint
or a shapely object.
scale: one of {'self', 'highest', 'lowest'} or (nbr, nbr) or nbr
'self': Output Footprint's resolution is the same as self
'highest': Output Footprint's resolution is the highest one among the input Footprints
'lowest': Output Footprint's resolution is the lowest one among the input Footprints
(nbr, nbr): Signed pixel size, aka scale
nbr: Signed pixel width. Signed pixel height is assumed to be -width
rotation: one of {'auto', 'fit'} or nbr
'auto'
If `scale` designate a Footprint object, its rotation is chosen
Else, self's rotation is chosen
'fit'
Output Footprint is the rotated minimum bounding rectangle
nbr
Angle in degree
alignment: {'auto', 'tl', (nbr, nbr)}
'auto'
If `scale` and `rotation` designate the same Footprint object, its alignment
is chosen
Else, 'tl' alignment is chosen
'tl': Ouput Footprint's alignement is the top left most point of the bounding rectangle
of the intersection
(nbr, nbr): Coordinate of a point that lie on the grid.
This point can be anywhere in space.
homogeneous: bool
False: No effect
True: Raise an exception if all input Footprints do not lie on the same grid as self.
Returns
-------
Footprint
"""
# Retrieve *objects
objects = list(objects)
if not objects:
raise ValueError('No other Footprint provided arguments')
def _ensure_shapely(geom):
if not isinstance(geom, sg.base.BaseGeometry):
geom = sg.shape(geom.__geo_interface__)
return geom
def _partition(pred, iterable):
t1, t2 = itertools.tee(iterable)
return list(filterfalse(pred, t1)), list(map(_ensure_shapely, filter(pred, t2)))
def _classify(obj):
if isinstance(obj, self.__class__):
return False
if hasattr(obj, '__geo_interface__'):
return True
raise TypeError("intersection() argument must be a Footprint or a geometry, not %s" % (
type(obj)
))
footprints, geoms = _partition(_classify, [self] + objects)
# Retrieve resolution
resolution = kwargs.pop('scale', 'self')
if isinstance(resolution, str):
if resolution not in self._INTERSECTION_RESOLUTIONS:
raise ValueError('bad resolution parameter')
else:
resolution = np.asarray(resolution, dtype='float64')
if resolution.ndim == 0:
resolution = np.asarray([resolution, -resolution], dtype='float64')
elif resolution.shape == (1,):
resolution = np.asarray([resolution[0], -resolution[0]], dtype='float64')
elif resolution.shape != (2,):
raise ValueError('resolution has shape {}'.format(resolution.shape))
if (resolution == 0).any():
raise ValueError('resolution should be different than zero')
# Retrieve rotation
rotation = kwargs.pop('rotation', 'auto')
if isinstance(rotation, str):
if rotation not in self._INTERSECTION_ROTATIONS:
raise ValueError()
else:
rotation = float(rotation)
# Retrieve alignment
alignment = kwargs.pop('alignment', 'auto')
if isinstance(alignment, str):
if alignment not in self._INTERSECTION_ALIGNMENTS:
raise ValueError('Unknown alignment value')
else:
alignment = np.asarray(alignment, dtype='float64')
if alignment.shape != (2,):
raise ValueError('alignment has shape {}'.format(alignment.shape))
# Retrieve homogeneous
homogeneous = bool(kwargs.pop('homogeneous', False))
if homogeneous:
for fp in footprints:
if not self.same_grid(fp):
raise ValueError(
'{} does not lie on the same grid as self: {}'.format(fp, self)
)
if kwargs:
raise ValueError('Unknown keyword arguments')
return self._intersection_unsafe(
footprints, geoms, resolution, rotation, alignment
)
[docs] def move(self, tl, tr=None, br=None, round_coordinates=False):
"""Create a copy of self moved by an Affine transformation by providing new points.
`rsize` is always conserved
Usage cases
-----------
+-------+-------+-------+-------------------------------------------------------------------+
| tl | tr | br | Affine transformations possible |
+=======+=======+=======+===================================================================+
| coord | None | None | Translation |
+-------+-------+-------+-------------------------------------------------------------------+
| coord | coord | None | Translation, Rotation, Scale x and y uniformly with positive real |
+-------+-------+-------+-------------------------------------------------------------------+
| coord | coord | coord | Translation, Rotation, Scale x and y independently with reals |
+-------+-------+-------+-------------------------------------------------------------------+
Parameters
----------
tl: (nbr, nbr)
New top left coordinates
tr: (nbr, nbr)
New top right coordinates
br: (nbr, nbr)
New bottom right coordinates
round_coordinates: bool
Round the input coordinates with respect to `buzz.env.significant`, so that the output
`Footprint` is as much similar as possible as the input `Footprint` regarding those
properties:
- angle
- pxsize
- pxsizex / pxsizey
This option helps a lot if the input coordinates suffered from floating point
precision loss since it will cancel the noise in the resulting transformation matrix.
.. warning::
Only work when `tr` and `br` are both provided
Returns
-------
Footprint
"""
tl = np.asarray(tl, dtype=np.float64)
if tl.shape != (2,):
raise ValueError('Bad tl shape') # pragma: no cover
if tr is not None:
tr = np.asarray(tr, dtype=np.float64)
if tr.shape != (2,):
raise ValueError('Bad tr shape') # pragma: no cover
if br is not None:
br = np.asarray(br, dtype=np.float64)
if br.shape != (2,):
raise ValueError('Bad br shape') # pragma: no cover
else:
if br is not None:
raise ValueError('If br present, bl should be present too') # pragma: no cover
if round_coordinates:
if br is None: # pragma: no cover
raise ValueError(
'Can only round input coordinates when all parameters are provided'
)
tl, tr, br = self._snap_target_coordinates_before_move(tl, tr, br)
if tr is None:
angle = self.angle
scale = self.scale
elif br is None:
lrvec = tr - tl
angle = float(np.arctan2(lrvec[1], lrvec[0]) * 180. / np.pi)
scale = np.linalg.norm(lrvec) / np.linalg.norm(self.lrvec) * self.scale
else:
rect = _tools.Rect(tl, tl + (br - tr), br, tr)
scale = rect.scale(self.rsize)
angle = rect.angle
significant_min = rect.significant_min(np.abs(scale).min())
if env.significant <= significant_min:
s = ('This Footprint have large coordinates and small pixels, at least {:.2} '
'significant digits are necessary to perform this operation, but '
'`buzz.env.significant` is set to {}. Increase this value by using '
'buzz.Env(significant={}) in a `with statement`.'
).format(significant_min, env.significant, env.significant + 1)
raise RuntimeError(s)
slack_angles = rect.tr_slack_angles
assert slack_angles[0] < slack_angles[1]
if np.prod(np.sign(slack_angles - 90)) != -1:
raise ValueError(
'tl-tr-br angle is between {} and {} degree (should be <90 and >90)'.format(*slack_angles)
)
aff = (
affine.Affine.translation(*tl) *
affine.Affine.rotation(angle) *
affine.Affine.scale(*scale)
)
try:
return self.__class__(
gt=aff.to_gdal(),
rsize=self.rsize,
)
except ValueError as e:
if br is not None and round_coordinates is False and \
len(e.args) > 0 and 'north-up' in e.args[0]:
raise ValueError('Moving Footprint failed. Try using `round_coordinates=True`.')
raise
# Export ************************************************************************************ **
@property
def extent(self):
"""Get the Footprint's extent (`x` then `y`)
Example
-------
>>> minx, maxx, miny, maxy = fp.extent
>>> plt.imshow(arr, extent=fp.extent)
fp.extent from fp.bounds using numpy fancy indexing
>>> minx, maxx, miny, maxy = fp.bounds[[0, 2, 1, 3]]
"""
points = np.r_["1,0,2", self.coords]
return np.asarray([
points[:, 0].min(), points[:, 0].max(),
points[:, 1].min(), points[:, 1].max(),
])
@property
def bounds(self):
"""Get the Footprint's bounds (`min` then `max`)
Example
-------
>>> minx, miny, maxx, maxy = fp.bounds
fp.bounds from fp.extent using numpy fancy indexing
>>> minx, miny, maxx, maxy = fp.extent[[0, 2, 1, 3]]
"""
points = np.r_["1,0,2", self.coords]
return np.asarray([
points[:, 0].min(), points[:, 1].min(),
points[:, 0].max(), points[:, 1].max(),
])
@property
def coords(self):
"""Get corners coordinates
Example
-------
>>> tl, bl, br, tr = fp.coords
"""
return np.asarray(
[self.tl, self.bl, self.br, self.tr]
)
@property
def poly(self):
"""Convert self to shapely.geometry.Polygon"""
return sg.Polygon([
self.tl, self.bl, self.br, self.tr, self.tl
])
@property
def __geo_interface__(self):
return {
'type': 'Polygon',
'coordinates': [[
[self.tlx, self.tly],
[self.blx, self.bly],
[self.brx, self.bry],
[self.trx, self.try_],
[self.tlx, self.tly],
]],
}
# Accessors ********************************************************************************* **
# Accessors - Spatial - Size and vectors **************************************************** **
@property
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)
@property
def sizex(self):
"""Spatial distance: ||raster left - raster right||"""
return float(self.size[0])
@property
def sizey(self):
"""Spatial distance: ||raster top - raster bottom||"""
return float(self.size[1])
@property
def width(self):
"""Spatial distance: ||raster left - raster right||, alias for sizex"""
return float(self.sizex)
@property
def height(self):
"""Spatial distance: ||raster top - raster bottom||, alias for sizey"""
return float(self.sizey)
@property
def w(self):
"""Spatial distance: ||raster left - raster right||, alias for sizex"""
return float(self.sizex)
@property
def h(self):
"""Spatial distance: ||raster top - raster bottom||, alias for sizey"""
return float(self.sizey)
@property
def lrvec(self):
"""Spatial vector: (raster right - raster left)"""
return self.tr - self.tl
@property
def tbvec(self):
"""Spatial vector: (raster bottom - raster top)"""
return self.bl - self.tl
@property
def diagvec(self):
"""Spatial vector: (raster bottom right - raster top left)"""
return self.br - self.tl
# Accessors - Spatial - Coordinates ********************************************************* **
@property
def tl(self):
"""Spatial coordinates: raster top left (x, y)"""
return self._tl.copy()
@property
def tlx(self):
"""Spatial coordinate: raster top left (x)"""
return float(self._tl[0])
@property
def tly(self):
"""Spatial coordinate: raster top left (y)"""
return float(self._tl[1])
@property
def bl(self):
"""Spatial coordinates: raster bottom left (x, y)"""
return self._bl.copy()
@property
def blx(self):
"""Spatial coordinate: raster bottom left (x)"""
return float(self._bl[0])
@property
def bly(self):
"""Spatial coordinate: raster bottom left (y)"""
return float(self._bl[1])
@property
def br(self):
"""Spatial coordinates: raster bottom right (x, y)"""
return self._br.copy()
@property
def brx(self):
"""Spatial coordinate: raster bottom right (x)"""
return float(self._br[0])
@property
def bry(self):
"""Spatial coordinate: raster bottom right (y)"""
return float(self._br[1])
@property
def tr(self):
"""Spatial coordinates: raster top right (x, y)"""
return self._tr.copy()
@property
def trx(self):
"""Spatial coordinate: raster top right (x)"""
return float(self._tr[0])
@property
def try_(self):
"""Spatial coordinate: raster top right (y)
Don't forget the trailing underscore
"""
return float(self._tr[1])
@property
def t(self):
"""Spatial coordinates: raster top center (x, y)"""
return np.array([self.tx, self.ty], dtype=np.float64)
@property
def tx(self):
"""Spatial coordinate: raster top center (x)"""
return (self.tlx + self.trx) / 2.
@property
def ty(self):
"""Spatial coordinate: raster top center (y)"""
return (self.tly + self.try_) / 2.
@property
def l(self):
"""Spatial coordinates: raster center left (x, y)"""
return np.array([self.lx, self.ly], dtype=np.float64)
@property
def lx(self):
"""Spatial coordinate: raster center left (x)"""
return (self.tlx + self.blx) / 2.
@property
def ly(self):
"""Spatial coordinate: raster center left (y)"""
return (self.tly + self.bly) / 2.
@property
def b(self):
"""Spatial coordinates: raster bottom center (x, y)"""
return np.array([self.bx, self.by], dtype=np.float64)
@property
def bx(self):
"""Spatial coordinate: raster bottom center (x)"""
return (self.blx + self.brx) / 2.
@property
def by(self):
"""Spatial coordinate: raster bottom center (y)"""
return (self.bly + self.bry) / 2.
@property
def r(self):
"""Spatial coordinates: raster center right (x, y)"""
return np.array([self.rx, self.ry], dtype=np.float64)
@property
def rx(self):
"""Spatial coordinate: raster center right (x)"""
return (self.brx + self.trx) / 2.
@property
def ry(self):
"""Spatial coordinate: raster center right (y)"""
return (self.bry + self.try_) / 2.
@property
def c(self):
"""Spatial coordinates: raster center (x, y)"""
return np.array([self.cx, self.cy], dtype=np.float64)
@property
def cx(self):
"""Spatial coordinate: raster center (x)"""
return (self.tx + self.bx) / 2
@property
def cy(self):
"""Spatial coordinate: raster center (y)"""
return (self.ly + self.ry) / 2
# Accessors - Spatial - Misc **************************************************************** **
@property
def semiminoraxis(self):
"""Spatial distance: half-size of the smaller side"""
return float(np.min(self.size) / 2.)
@property
def semimajoraxis(self):
"""Spatial distance: half-size of the bigger side"""
return float(np.max(self.size) / 2.)
@property
def area(self):
"""Area: pixel count"""
return np.prod(self.size)
@property
def length(self):
"""Length: circumference of the outer ring"""
return np.sum(self.size) * 2
# Accessors - Raster - Size ***************************************************************** **
@property
def rsize(self):
"""Pixel quantities: (pixel per line, pixel per column)"""
return self._rsize.copy()
@property
def rsizex(self):
"""Pixel quantity: pixel per line"""
return int(self._rsize[0])
@property
def rsizey(self):
"""Pixel quantity: pixel per column"""
return int(self._rsize[1])
@property
def rwidth(self):
"""Pixel quantity: pixel per line, alias for rsizex"""
return int(self.rsizex)
@property
def rheight(self):
"""Pixel quantity: pixel per column, alias for rsizey"""
return int(self.rsizey)
@property
def rw(self):
"""Pixel quantity: pixel per line, alias for rsizex"""
return int(self.rsizex)
@property
def rh(self):
"""Pixel quantity: pixel per column, alias for rsizey"""
return int(self.rsizey)
# Accessors - Raster - Indices ************************************************************** **
@property
def rtl(self):
"""Indices: raster top left pixel (x=0, y=0)"""
return np.array([0, 0], dtype=env.default_index_dtype)
@property
def rtlx(self):
"""Index: raster top left pixel (x=0)"""
return 0
@property
def rtly(self):
"""Index: raster top left pixel (y=0)"""
return 0
@property
def rbl(self):
"""Indices: raster bottom left pixel (x=0, y)"""
return np.array([self.rblx, self.rbly], dtype=env.default_index_dtype)
@property
def rblx(self):
"""Index: raster bottom left pixel (x=0)"""
return 0
@property
def rbly(self):
"""Index: raster bottom left pixel (y)"""
return self.rsizey - 1
@property
def rbr(self):
"""Indices: raster bottom right pixel (x, y)"""
return np.array([self.rbrx, self.rbry], dtype=env.default_index_dtype)
@property
def rbrx(self):
"""Index: raster bottom right pixel (x)"""
return self.rsizex - 1
@property
def rbry(self):
"""Index: raster bottom right pixel (y)"""
return self.rsizey - 1
@property
def rtr(self):
"""Indices: raster top right pixel (x, y=0)"""
return np.array([self.rtrx, self.rtry], dtype=env.default_index_dtype)
@property
def rtrx(self):
"""Index: raster top right pixel (x)"""
return self.rsizex - 1
@property
def rtry(self):
"""Index: raster top right pixel (y=0)"""
return 0
@property
def rt(self):
"""Indices: raster top center pixel (x truncated, y=0)"""
return np.array([self.rtx, self.rty], dtype=env.default_index_dtype)
@property
def rtx(self):
"""Index: raster top center pixel (x truncated)"""
return int((self.rsizex - 1) / 2)
@property
def rty(self):
"""Index: raster top center pixel (y=0)"""
return 0
@property
def rl(self):
"""Indices: raster center left pixel (x=0, y truncated)"""
return np.array([self.rlx, self.rly], dtype=env.default_index_dtype)
@property
def rlx(self):
"""Index: raster center left pixel (x=0)"""
return 0
@property
def rly(self):
"""Index: raster center left pixel (y truncated)"""
return int((self.rsizey - 1) / 2)
@property
def rb(self):
"""Indices: raster bottom center pixel (x truncated, y)"""
return np.array([self.rbx, self.rby], dtype=env.default_index_dtype)
@property
def rbx(self):
"""Index: raster bottom center pixel (x truncated)"""
return int((self.rsizex - 1) / 2)
@property
def rby(self):
"""Index: raster bottom center pixel (y)"""
return self.rsizey - 1
@property
def rr(self):
"""Indices: raster center right pixel (x, y truncated)"""
return np.array([self.rrx, self.rry], dtype=env.default_index_dtype)
@property
def rrx(self):
"""Index: raster center right pixel (x)"""
return self.rsizex - 1
@property
def rry(self):
"""Index: raster center right pixel (y truncated)"""
return int((self.rsizey - 1) / 2)
@property
def rc(self):
"""Indices: raster center pixel (x truncated, y truncated)"""
return np.array([self.rcx, self.rcy], dtype=env.default_index_dtype)
@property
def rcx(self):
"""Index: raster center pixel (x truncated)"""
return int((self.rsizex - 1) / 2)
@property
def rcy(self):
"""Index: raster center pixel (y truncated)"""
return int((self.rsizey - 1) / 2)
# Accessors - Raster - Misc ***************************************************************** **
@property
def rsemiminoraxis(self):
"""Pixel quantity: half pixel count (truncated) of the smaller side"""
return int(np.min(self.rsize) / 2.)
@property
def rsemimajoraxis(self):
"""Pixel quantity: half pixel count (truncated) of the bigger side"""
return int(np.max(self.rsize) / 2.)
@property
def rarea(self):
"""Pixel quantity: pixel count"""
rx, ry = self.rsize
# Convert to int before multiplication to avoid overflow
return int(rx) * int(ry)
@property
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 ******************************************************** **
@property
def gt(self):
"""First 6 numbers of the affine transformation matrix, GDAL ordering"""
return np.array(self._aff.to_gdal(), dtype=np.float64)
@property
def aff33(self):
"""The affine transformation matrix"""
return np.asarray(self._aff, dtype=np.float64).reshape(3, 3)
@property
def aff23(self):
"""Top two rows of the affine transformation matrix"""
return self.aff33[:2]
@property
def aff6(self):
"""First 6 numbers of the affine transformation matrix, left-right/top-bottom ordering"""
return self.aff23.reshape(6)
@property
def affine(self):
"""Underlying affine object"""
# TODO: Check if _aff is mutable
return self._aff
@property
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)
@property
def angle(self):
"""Angle in degree: rotation used in the affine transformation, (0 is north-up)"""
lrvec = self.lrvec
return float(np.arctan2(lrvec[1], lrvec[0]) * 180. / np.pi)
@property
def pxsize(self):
"""Spatial distance: ||pixel bottom right - pixel top left|| (x, y)"""
return self.size / self.rsize
@property
def pxsizex(self):
"""Spatial distance: ||pixel right - pixel left|| (x)"""
return self.sizex / self.rsizex
@property
def pxsizey(self):
"""Spatial distance: ||pixel bottom - pixel top|| (y)"""
return self.sizey / self.rsizey
@property
def pxvec(self):
"""Spatial vector: (pixel bottom right - pixel top left)"""
return self.diagvec / self.rsize
@property
def pxtbvec(self):
"""Spatial vector: (pixel bottom left - pixel top left)"""
return self.tbvec / self.rheight
@property
def pxlrvec(self):
"""Spatial vector: (pixel top right - pixel top left)"""
return self.lrvec / self.rwidth
# Binary predicates ************************************************************************* **
# Binary predicates - geometry ************************************************************** **
[docs] def equals(self, other):
"""Binary predicate: Is other Footprint exactly equal to self
Parameters
----------
other: Footprint
..
Returns
-------
bool
"""
if (self.gt != other.gt).any():
return False
if (self.rsize != other.rsize).any():
return False
return True
[docs] def almost_equals(self, other):
"""Binary predicate: Is other Footprint almost equal to self with regard to
buzz.env.significant.
Parameters
----------
other: Footprint
..
Returns
-------
bool
"""
if env.significant <= self._significant_min:
s = ('This Footprint have large coordinates and small pixels, at least {:.2} '
'significant digits are necessary to perform this operation, but '
'`buzz.env.significant` is set to {}. Increase this value by using '
'buzz.Env(significant={}) in a `with statement`.'
).format(self._significant_min, env.significant, env.significant + 1)
raise RuntimeError(s)
if env.significant <= other._significant_min:
s = ('This Footprint have large coordinates and small pixels, at least {:.2} '
'significant digits are necessary to perform this operation, but '
'`buzz.env.significant` is set to {}. Increase this value by using '
'buzz.Env(significant={}) in a `with statement`.'
).format(self._significant_min, env.significant, env.significant + 1)
raise RuntimeError(s)
if (self.rsize != other.rsize).any():
return False
largest_coord = np.abs(np.r_[self.coords, other.coords]).max()
spatial_precision = largest_coord * 10 ** -env.significant
return (np.abs(self.coords - other.coords) < spatial_precision).all()
[docs] def same_grid(self, other):
"""Binary predicate: Does other Footprint lie on the same grid as self
Parameters
----------
other: Footprint
..
Returns
-------
bool
"""
if env.significant <= self._significant_min:
s = ('This Footprint have large coordinates and small pixels, at least {:.2} '
'significant digits are necessary to perform this operation, but '
'`buzz.env.significant` is set to {}. Increase this value by using '
'buzz.Env(significant={}) in a `with statement`.'
).format(self._significant_min, env.significant, env.significant + 1)
raise RuntimeError(s)
if env.significant <= other._significant_min:
s = ('This Footprint have large coordinates and small pixels, at least {:.2} '
'significant digits are necessary to perform this operation, but '
'`buzz.env.significant` is set to {}. Increase this value by using '
'buzz.Env(significant={}) in a `with statement`.'
).format(self._significant_min, env.significant, env.significant + 1)
raise RuntimeError(s)
largest_coord = np.abs(np.r_[self.coords, other.coords]).max()
spatial_precision = largest_coord * 10 ** -env.significant
rdx, rdy = np.around(~self.affine * other.tl)
errors = other.tl - (self.pxtbvec * rdy + self.pxlrvec * rdx) - self.tl
if (np.abs(errors) >= spatial_precision).any():
return False
errors = self.tl + other.pxtbvec * self.rheight - self.bl
if (np.abs(errors) >= spatial_precision).any():
return False
errors = self.tl + other.pxlrvec * self.rwidth - self.tr
if (np.abs(errors) >= spatial_precision).any():
return False
errors = other.tl + self.pxtbvec * other.rheight - other.bl
if (np.abs(errors) >= spatial_precision).any():
return False
errors = other.tl + self.pxlrvec * other.rwidth - other.tr
if (np.abs(errors) >= spatial_precision).any():
return False
return True
# Numpy ************************************************************************************* **
@property
def shape(self):
"""Pixel quantities: (pixel per column, pixel per line)"""
return np.flipud(self._rsize)
@property
def meshgrid_raster(self):
"""Compute indice matrices
Returns
-------
(x, y): (np.ndarray, np.ndarray)
Raster indices matrices
with shape = self.shape
with dtype = env.default_index_dtype
"""
return np.meshgrid(
np.asarray(range(self.rsizex), dtype=env.default_index_dtype),
np.asarray(range(self.rsizey), dtype=env.default_index_dtype),
copy=False,
)
@property
def meshgrid_spatial(self):
"""Compute coordinate matrices
Returns
-------
(x, y): (np.ndarray, np.ndarray)
Spatial coordinate matrices
with shape = self.shape
with dtype = float32
"""
x, y = self.meshgrid_raster
return (
(x * self._aff.a + y * self._aff.b + self._aff.c).astype(np.float64, copy=False),
(x * self._aff.d + y * self._aff.e + self._aff.f).astype(np.float64, copy=False),
)
[docs] def meshgrid_raster_in(self, other, dtype=None, op=np.floor):
"""Compute raster coordinate matrices of `self` in `other` referential
Parameters
----------
other: Footprint
..
dtype: None or convertible to np.dtype
Output dtype
If None: Use buzz.env.default_index_dtype
op: None or function operating on a vector
Function to apply before casting output to dtype
If None: Do not transform data before casting
Returns
-------
(x, y): (np.ndarray, np.ndarray)
Raster coordinate matrices
with shape = self.shape
with dtype = dtype
"""
# Check other parameter
if not isinstance(other, self.__class__):
raise TypeError('other should be a Footprint') # pragma: no cover
# Check dtype parameter
if dtype is None:
dtype = env.default_index_dtype
else:
dtype = conv.dtype_of_any_downcast(dtype)
# Check op parameter
if not np.issubdtype(dtype, np.integer):
op = None
xy = other.spatial_to_raster(np.dstack(self.meshgrid_spatial), dtype=dtype, op=op)
return xy[..., 0], xy[..., 1]
[docs] def slice_in(self, other, clip=False):
"""Compute location of `self` inside `other` with slice objects.
If other and self do not have the same rotation, operation is undefined
Parameters
----------
other: Footprint
..
clip: bool
False
Does nothing
True
Clip the slices to other bounds. If other and self do not share area,
at least one of the returned slice will have `slice.start == slice.stop`
Returns
-------
(yslice, xslice): (slice, slice)
Example
-------
Burn `small` into `big` if `small` is within `big`
>>> big_data[small.slice_in(big)] = small_data
Burn `small` into `big` where overlapping
>>> big_data[small.slice_in(big, clip=True)] = small_data[big.slice_in(small, clip=True)]
"""
if not isinstance(other, self.__class__):
raise TypeError('other should be a Footprint') # pragma: no cover
startx, starty = other.spatial_to_raster(self.tl)
endx, endy = other.spatial_to_raster(self.br)
if clip:
startx = startx.clip(0, other.rsizex)
endx = endx.clip(0, other.rsizex)
starty = starty.clip(0, other.rsizey)
endy = endy.clip(0, other.rsizey)
return slice(starty, endy), slice(startx, endx)
# Coordinates conversions ******************************************************************* **
[docs] def spatial_to_raster(self, xy, dtype=None, op=np.floor):
"""Convert xy spatial coordinates to raster xy indices
Parameters
----------
xy: sequence of numbers of shape (..., 2)
Spatial coordinates
dtype: None or convertible to np.dtype
Output dtype
If None: Use buzz.env.default_index_dtype
op: None or vectorized function
Function to apply before casting output to dtype
If None: Do not transform data before casting
Returns
-------
out_xy: np.ndarray
Raster indices
with shape = np.asarray(xy).shape
with dtype = dtype
Prototype inspired from
https://mapbox.github.io/rasterio/api/rasterio.io.html#rasterio.io.TransformMethodsMixin.index
"""
# Check xy parameter
xy = np.asarray(xy)
if xy.shape[-1] != 2:
raise ValueError('An array of shape (..., 2) was expected') # pragma: no cover
# Check dtype parameter
if dtype is None:
dtype = env.default_index_dtype
else:
dtype = conv.dtype_of_any_downcast(dtype)
# Check op parameter
if not np.issubdtype(dtype, np.integer):
op = None
if env.significant <= self._significant_min:
s = ('This Footprint have large coordinates and small pixels, at least {:.2} '
'significant digits are necessary to perform this operation, but '
'`buzz.env.significant` is set to {}. Increase this value by using '
'buzz.Env(significant={}) in a `with statement`.'
).format(self._significant_min, env.significant, env.significant + 1)
raise RuntimeError(s)
largest_coord = np.abs(self.coords).max()
spatial_precision = largest_coord * 10 ** -env.significant
smallest_reso = self.pxsize.min()
pixel_precision = spatial_precision / smallest_reso
abstract_grid_density = np.floor(1 / pixel_precision)
workshape = int(xy.size / 2), 2
xy2 = np.empty(workshape, 'float64')
xy2[:, :] = xy.reshape(workshape)
aff = ~self._aff
xy2[:, 0], xy2[:, 1] = (
xy2[:, 0] * aff.a + xy2[:, 1] * aff.b + aff.c,
xy2[:, 0] * aff.d + xy2[:, 1] * aff.e + aff.f,
)
xy2 = np.around(xy2 * abstract_grid_density, 0) / abstract_grid_density # Should move this line in if?
if op is not None:
xy2 = op(xy2)
return xy2.astype(dtype, copy=False).reshape(xy.shape)
[docs] def raster_to_spatial(self, xy):
"""Convert xy raster coordinates to spatial coordinates
Parameters
----------
xy: sequence of numbers of shape (..., 2)
Raster coordinages
Returns
-------
out_xy: np.ndarray
Spatial coordinates
with shape = np.asarray(xy).shape
with dtype = dtype
"""
# Check xy parameter
xy = np.asarray(xy)
if xy.shape[-1] != 2:
raise ValueError('An array of shape (..., 2) was expected') # pragma: no cover
workshape = int(xy.size / 2), 2
xy2 = np.empty(workshape, 'float64')
xy2[:, :] = xy.reshape(workshape)
aff = self._aff
xy2[:, 0], xy2[:, 1] = (
xy2[:, 0] * aff.a + xy2[:, 1] * aff.b + aff.c,
xy2[:, 0] * aff.d + xy2[:, 1] * aff.e + aff.f,
)
return xy2.reshape(xy.shape)
# Geometry / Raster conversions ************************************************************* **
[docs] def find_lines(self, arr, output_offset='middle', merge=True):
"""Create a list of line-strings from a mask. Works with connectivity 4 and 8. The input
raster is preprocessed using `skimage.morphology.thin`. The output linestrings are
postprocessed using `shapely.ops.linemerge`.
.. warning::
All standalone pixels contained in arr will be ignored.
Parameters
----------
arr: np.ndarray of bool of shape (self.shape)
..
output_offset: 'middle' or (nbr, nbr)
Coordinate offset in meter
if `middle`: substituted by `self.pxvec / 2`
Returns
-------
list of shapely.geometry.LineString
Exemple
-------
>>> import buzzard as buzz
>>> import numpy as np
>>> import networkx as nx
>>> with buzz.Env(allow_complex_footprint=1):
... a = np.asarray([
... [0, 1, 1, 1, 0],
... [0, 1, 0, 0, 0],
... [0, 1, 1, 1, 0],
... [0, 1, 0, 0, 0],
... [0, 1, 1, 0, 0],
...
... ])
... fp = buzz.Footprint(gt=(0, 1, 0, 0, 0, 1), rsize=(a.shape))
... lines = fp.find_lines(a, (0, 0))
...
... # Display input / output
... print(fp)
... print(a.astype(int))
... for i, l in enumerate(lines, 1):
... print(f'edge-id:{i} of type:{type(l)} and length:{l.length}')
... print(fp.burn_lines(l).astype(int) * i)
...
... # Build a networkx graph
... g = nx.Graph([(l.coords[0], l.coords[-1]) for l in lines])
... print(repr(g.degree))
...
Footprint(tl=(0.000000, 0.000000), scale=(1.000000, 1.000000), angle=0.000000, rsize=(5, 5))
[[0 1 1 1 0]
[0 1 0 0 0]
[0 1 1 1 0]
[0 1 0 0 0]
[0 1 1 0 0]]
edge-id:1 of type:<class 'shapely.geometry.linestring.LineString'> and length:2.0
[[0 0 0 0 0]
[0 0 0 0 0]
[0 1 1 1 0]
[0 0 0 0 0]
[0 0 0 0 0]]
edge-id:2 of type:<class 'shapely.geometry.linestring.LineString'> and length:3.0
[[0 0 0 0 0]
[0 0 0 0 0]
[0 2 0 0 0]
[0 2 0 0 0]
[0 2 2 0 0]]
edge-id:3 of type:<class 'shapely.geometry.linestring.LineString'> and length:4.0
[[0 3 3 3 0]
[0 3 0 0 0]
[0 3 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]]
DegreeView({(3.0, 2.0): 1, (1.0, 2.0): 3, (2.0, 4.0): 1, (3.0, 0.0): 1})
"""
# Step 1: Parameter checking ************************************************************ **
if arr.shape != tuple(self.shape):
raise ValueError('Incompatible shape between array:%s and self:%s' % (
arr.shape, self.shape
)) # pragma: no cover
if output_offset == 'middle':
output_offset = self.pxvec / 2
else:
output_offset = np.asarray(output_offset, 'float64')
if output_offset.shape != (2,):
raise ValueError(
'`output_offset` should be "middle" or a sequence of float of length 2'
) # pragma: no cover
# Step 2: Array normalization *********************************************************** **
arr = arr.astype(bool, copy=False)
arr = skm.thin(arr)
arr = arr.astype('uint8', copy=False)
# Step 3: Prepare to collapse 2x2 squares to single points ****************************** **
squares2x2_topleft_mask = ndi.convolve(
arr,[
[1, 1, 0],
[1, 1, 0],
[0, 0, 0],
],
mode='constant',
) == 4
squares2x2_yx_links = {
(int(y + dy), int(x + dx)): (int(y), (x))
for y, x in zip(*squares2x2_topleft_mask.nonzero())
for dy in range(2)
for dx in range(2)
# if dy + dx != 0
}
# Step 4: Retrieve pixel indices ******************************************************** **
count = np.sum(arr)
yx_lst = np.stack(arr.nonzero(), -1)
index_lst = np.arange(count)
index = np.empty(self.shape, dtype=int)
index[arr != 0] = index_lst
# Step 5: Retrieve edge indices ********************************************************* **
convolve = lambda arr, kernel: scipy.ndimage.convolve(arr, kernel, mode='constant', cval=0)
has_top = convolve(arr, [[0, 0, 0], [0, 0, 0], [0, 1, 0]]) * arr
has_right = convolve(arr, [[0, 0, 0], [1, 0, 0], [0, 0, 0]]) * arr
has_left = convolve(arr, [[0, 0, 0], [0, 0, 1], [0, 0, 0]]) * arr
has_topright = convolve(arr, [[0, 0, 0], [0, 0, 0], [1, 0, 0]]) * arr
has_topleft = convolve(arr, [[0, 0, 0], [0, 0, 0], [0, 0, 1]]) * arr
has_topright = (has_topright.astype('i1', copy=False) - has_right - has_top).clip(0, 1).astype('u1', copy=False)
has_topleft = (has_topleft.astype('i1', copy=False) - has_left - has_top).clip(0, 1).astype('u1', copy=False)
def _build_neighbors_in_direction(mask, yx_vector):
has_indices = mask[yx_lst[:, 0], yx_lst[:, 1]].nonzero()[0]
neig_yx = yx_lst[has_indices] + yx_vector
neig = index[neig_yx[:, 0], neig_yx[:, 1]]
return np.c_[neig, has_indices]
edges_indices = np.vstack([
_build_neighbors_in_direction(has_top, (-1, 0)),
_build_neighbors_in_direction(has_right, (0, 1)),
_build_neighbors_in_direction(has_topright, (-1, 1)),
_build_neighbors_in_direction(has_topleft, (-1, -1)),
])
if edges_indices.size == 0:
return []
# Step 6: Convert segments made of indices to segments made of shapely objets *********** **
lines = []
for (n1, n2) in edges_indices:
yx1 = tuple(yx_lst[n1].tolist())
yx2 = tuple(yx_lst[n2].tolist())
if yx1 in squares2x2_yx_links and yx2 in squares2x2_yx_links:
# Drop segments inside 2x2 squares
continue
points = [
yx1,
yx2,
]
if points[0] in squares2x2_yx_links:
points.insert(0, squares2x2_yx_links[points[0]])
if points[-1] in squares2x2_yx_links:
points.append(squares2x2_yx_links[points[-1]])
l = shapely.geometry.LineString(self.raster_to_spatial(
list(map(np.flipud, points))
) + output_offset)
lines.append(l)
# Step 7: Merge segments to polylines *************************************************** **
if merge:
mline = shapely.ops.linemerge(lines)
if isinstance(mline, shapely.geometry.LineString):
mline = sg.MultiLineString([mline])
# Step 8: Return ************************************************************************ **
if isinstance(mline, shapely.geometry.LineString):
return [mline]
if isinstance(mline, shapely.geometry.MultiLineString):
return list(mline.geoms)
assert False # pragma: no cover
[docs] def burn_lines(self, obj, all_touched=False, labelize=False):
"""Creates a 2d image from lines. Uses gdal.Polygonize.
Parameters
----------
obj: shapely line or nested iterators over shapely lines
..
labelize: bool
- if `False`: Create a boolean mask
- if `True`: Create an integer matrix containing lines indices from order in input
Returns
----------
np.ndarray
- of bool or uint8 or int
- of shape (self.shape)
"""
lines = list(_line_iterator(obj))
# https://svn.osgeo.org/gdal/trunk/autotest/alg/rasterize.py
sr_wkt = 'LOCAL_CS["arbitrary"]'
sr = osr.SpatialReference(sr_wkt)
if labelize:
if len(lines) > 255:
dtype = conv.dtype_of_any_downcast('int')
else:
dtype = conv.dtype_of_any_downcast('uint8')
else:
dtype = conv.dtype_of_any_downcast('bool')
gdt = conv.gdt_of_any_equiv(dtype) # Set to downcast
target_ds = gdal.GetDriverByName('MEM').Create(
'', int(self.rsizex), int(self.rsizey), 1, gdt
)
target_ds.SetGeoTransform(self.gt)
target_ds.SetProjection(sr_wkt)
rast_ogr_ds = ogr.GetDriverByName('Memory').CreateDataSource('wrk')
rast_mem_lyr = rast_ogr_ds.CreateLayer('line', srs=sr)
val_field = ogr.FieldDefn('val', ogr.OFTInteger64)
rast_mem_lyr.CreateField(val_field)
for i, line in enumerate(lines, 1):
feat = ogr.Feature(rast_mem_lyr.GetLayerDefn())
wkt_geom = line.wkt
feat.SetGeometryDirectly(ogr.Geometry(wkt=wkt_geom))
feat.SetFieldInteger64(0, i)
rast_mem_lyr.CreateFeature(feat)
if all_touched:
options = ["ALL_TOUCHED=TRUE, ATTRIBUTE=val"]
else:
options = ["ATTRIBUTE=val"]
success, payload = Catch(gdal.RasterizeLayer, nonzero_int_is_error=True)(
target_ds, [1], rast_mem_lyr, options=options
)
if not success:
raise ValueError('Could not rasterize (gdal error: `{}`)'.format(payload[1]))
arr = target_ds.GetRasterBand(1).ReadAsArray()
return arr.astype(dtype, copy=False)
[docs] def find_polygons(self, mask):
"""Creates a list of polygons from a mask. Uses gdal.Polygonize.
.. warning::
This method is not equivalent to `cv2.findContours` that considers that pixels are
points and therefore returns the indices of the pixels of the contours of the features.
This method consider that the pixels are areas and therefore returns the coordinates of
the points that surrounds the features.
.. warning::
Some inputs that may produce invalid polygons (see below) are fixed with the \
`shapely.geometry.Polygon.buffer` method.
Shapely will issue several warnings while buzzard fixes the polygons.
>>> # 0 0 0 0 0 0 0
... # 0 1 1 1 0 0 0
... # 0 1 1 1 1 0 0
... # 0 1 1 1 0 1 0 <- This feature has a hole near an edge. GDAL produces a self
... # 0 1 1 1 1 1 1 touching polygon without holes. A polygon with one hole is
... # 0 1 1 1 1 1 1 returned with this method.
... # 0 0 0 0 0 0 0
Parameters
----------
arr: np.ndarray of bool of shape (self.shape)
Returns
-------
list of shapely.geometry.Polygon
"""
if mask.shape != tuple(self.shape):
raise ValueError('Mask shape%s incompatible with self shape%s' % (
mask.shape, tuple(self.shape)
)) # pragma: no cover
mask = mask.astype('uint8', copy=False).clip(0, 1)
sr_wkt = 'LOCAL_CS["arbitrary"]'
sr = osr.SpatialReference(sr_wkt)
source_ds = gdal.GetDriverByName('MEM').Create(
'', int(self.rsizex), int(self.rsizey), 1, gdal.GDT_Byte
)
source_ds.SetGeoTransform(self.gt)
source_ds.SetProjection(sr_wkt)
source_ds.GetRasterBand(1).WriteArray(mask, 0, 0)
ogr_ds = ogr.GetDriverByName('Memory').CreateDataSource('wrk')
ogr_lyr = ogr_ds.CreateLayer('poly', srs=sr)
field_defn = ogr.FieldDefn('elev', ogr.OFTReal)
ogr_lyr.CreateField(field_defn)
success, payload = Catch(gdal.Polygonize, nonzero_int_is_error=True)(
srcBand=source_ds.GetRasterBand(1),
maskBand=source_ds.GetRasterBand(1),
outLayer=ogr_lyr,
iPixValField=0,
)
if not success:
raise ValueError('Could not polygonize (gdal error: `{}`)'.format(payload[1]))
del source_ds
def _polygon_iterator():
feat = ogr_lyr.GetNextFeature()
while feat is not None:
geometry = feat.geometry()
geometry = conv.shapely_of_ogr(geometry)
if not geometry.is_valid:
geometry = geometry.buffer(0)
yield geometry
feat = ogr_lyr.GetNextFeature()
return list(_polygon_iterator())
[docs] def burn_polygons(self, obj, all_touched=False, labelize=False):
"""Creates a 2d image from polygons. Uses gdal.RasterizeLayer.
.. warning::
This method is not equivalent to `cv2.drawContours` that considers that pixels are
points and therefore expect as input the indices of the outer pixels of each feature.
This method consider that the pixels are areas and therefore expect as input the
coordinates of the points surrounding the features.
Parameters
----------
obj: shapely polygon or nested iterators over shapely polygons
..
all_touched: bool
Burn all polygons touched
Returns
----------
np.ndarray
of bool or uint8 or int
of shape (self.shape)
Examples
--------
>>> burn_polygons(poly)
>>> burn_polygons([poly, poly])
>>> burn_polygons([poly, poly, [poly, poly], multipoly, poly])
"""
polys = list(_poly_iterator(obj))
# https://svn.osgeo.org/gdal/trunk/autotest/alg/rasterize.py
sr_wkt = 'LOCAL_CS["arbitrary"]'
sr = osr.SpatialReference(sr_wkt)
if labelize:
if len(polys) >= 65535:
dtype = conv.dtype_of_any_downcast('uint32')
elif len(polys) >= 255:
dtype = conv.dtype_of_any_downcast('uint16')
else:
dtype = conv.dtype_of_any_downcast('uint8')
else:
dtype = conv.dtype_of_any_downcast('bool')
gdt = conv.gdt_of_any_equiv(dtype) # Set to downcast
target_ds = gdal.GetDriverByName('MEM').Create(
'', int(self.rsizex), int(self.rsizey), 1, gdt
)
target_ds.SetGeoTransform(self.gt)
target_ds.SetProjection(sr_wkt)
rast_ogr_ds = ogr.GetDriverByName('Memory').CreateDataSource('wrk')
rast_mem_lyr = rast_ogr_ds.CreateLayer('polygon', srs=sr)
val_field = ogr.FieldDefn('val', ogr.OFTInteger64)
rast_mem_lyr.CreateField(val_field)
for i, poly in enumerate(polys, 1):
feat = ogr.Feature(rast_mem_lyr.GetLayerDefn())
wkt_geom = poly.wkt
feat.SetGeometryDirectly(ogr.Geometry(wkt=wkt_geom))
feat.SetFieldInteger64(0, i)
rast_mem_lyr.CreateFeature(feat)
if all_touched:
options = ["ALL_TOUCHED=TRUE, ATTRIBUTE=val"]
else:
options = ["ATTRIBUTE=val"]
success, payload = Catch(gdal.RasterizeLayer, nonzero_int_is_error=True)(
target_ds, [1], rast_mem_lyr, options=options
)
if not success:
raise ValueError('Could not rasterize (gdal error: `{}`)'.format(payload[1]))
arr = target_ds.GetRasterBand(1).ReadAsArray()
return arr.astype(dtype, copy=False)
# Tiling ************************************************************************************ **
[docs] def tile(self, size, overlapx=0, overlapy=0,
boundary_effect='extend', boundary_effect_locus='br'):
"""Tile a Footprint to a matrix of Footprint
Parameters
----------
size: (int, int)
Tile width and tile height, in pixel
overlapx: int
Width of a tile overlapping with each direct horizontal neighbors, in pixel
overlapy: int
Height of a tile overlapping with each direct vertical neighbors, in pixel
boundary_effect: {'extend', 'exclude', 'overlap', 'shrink', 'exception'}
Behevior at boundary effect locus
- 'extend'
- Preserve tile size
- Preserve overlapx and overlapy
- Sacrifice global bounds, results in tiles partially outside bounds at locus (if necessary)
- Preserve tile count
- Preserve boundary pixels coverage
- 'overlap'
- Preserve tile size
- Sacrifice overlapx and overlapy, results in tiles overlapping more at locus (if necessary)
- Preserve global bounds
- Preserve tile count
- Preserve boundary pixels coverage
- 'exclude'
- Preserve tile size
- Preserve overlapx and overlapy
- Preserve global bounds
- Sacrifice tile count, results in tiles excluded at locus (if necessary)
- Sacrifice boundary pixels coverage at locus (if necessary)
- 'shrink'
- Sacrifice tile size, results in tiles shrinked at locus (if necessary)
- Preserve overlapx and overlapy
- Preserve global bounds
- Preserve tile count
- Preserve boundary pixels coverage
- 'exception'
- Raise an exception if tiles at locus do not lie inside the global bounds
boundary_effect_locus: {'br', 'tr', 'tl', 'bl'}
Locus of the boundary effects
- 'br' : Boundary effect occurs at the bottom right corner of the raster, \
top left coordinates are preserved
- 'tr' : Boundary effect occurs at the top right corner of the raster, \
bottom left coordinates are preserved
- 'tl' : Boundary effect occurs at the top left corner of the raster, \
bottom right coordinates are preserved
- 'bl' : Boundary effect occurs at the bottom left corner of the raster, \
top right coordinates are preserved
Returns
-------
np.ndarray
- of dtype=object (Footprint)
- of shape (M, N)
- with M the line count
- with N the column count
"""
size = np.asarray(size, dtype=int)
overlapx = int(overlapx)
overlapy = int(overlapy)
if size.shape != (2,):
raise ValueError('size.shape(%s) should be (2,)' % str(size.shape))
if (size <= 0).any():
raise ValueError('size(%s) values should satisfy value > 0' % str(tuple(size)))
if not 0 <= overlapx < size[0]:
raise ValueError('overlapx(%d) should satisfy 0 <= overlapx < size[0](%d)' % (
overlapx, size[0]
))
if not 0 <= overlapy < size[1]:
raise ValueError('overlapy(%d) should satisfy 0 <= overlapy < size[1](%d)' % (
overlapy, size[1]
))
if boundary_effect not in self._TILE_BOUNDARY_EFFECTS:
raise ValueError('boundary_effect(%s) should be one of %s' % (
boundary_effect, self._TILE_BOUNDARY_EFFECTS
))
if boundary_effect_locus not in self._TILE_BOUNDARY_EFFECT_LOCI:
raise ValueError('boundary_effect_locus(%s) should be one of %s' % (
boundary_effect_locus, self._TILE_BOUNDARY_EFFECT_LOCI
))
return self._tile_unsafe(size, overlapx, overlapy, boundary_effect, boundary_effect_locus)
[docs] def tile_count(self, rowcount, colcount, overlapx=0, overlapy=0,
boundary_effect='extend', boundary_effect_locus='br'):
"""Tile a Footprint to a matrix of Footprint
Parameters
----------
rowcount: int
Tile count per row
colcount: int
Tile count per column
overlapx: int
Width of a tile overlapping with each direct horizontal neighbors, in pixel
overlapy: int
Height of a tile overlapping with each direct vertical neighbors, in pixel
boundary_effect: {'extend', 'exclude', 'overlap', 'shrink', 'exception'}
Behevior at boundary effect locus
- 'extend'
- Preserve tile size
- Preserve overlapx and overlapy
- Sacrifice global bounds, results in tiles partially outside bounds at locus (if necessary)
- Preserve tile count
- Preserve boundary pixels coverage
- 'overlap'
- Preserve tile size
- Sacrifice overlapx and overlapy, results in tiles overlapping more at locus (if necessary)
- Preserve global bounds
- Preserve tile count
- Preserve boundary pixels coverage
- 'exclude'
- Preserve tile size
- Preserve overlapx and overlapy
- Preserve global bounds
- Preserve tile count
- Sacrifice boundary pixels coverage at locus (if necessary)
- 'shrink'
- Sacrifice tile size, results in tiles shrinked at locus (if necessary)
- Preserve overlapx and overlapy
- Preserve global bounds
- Preserve tile count
- Preserve boundary pixels coverage
- 'exception'
- Raise an exception if tiles at locus do not lie inside the global bounds
boundary_effect_locus: {'br', 'tr', 'tl', 'bl'}
Locus of the boundary effects
- 'br' : Boundary effect occurs at the bottom right corner of the raster, \
top left coordinates are preserved
- 'tr' : Boundary effect occurs at the top right corner of the raster, \
bottom left coordinates are preserved
"""
rowcount = int(rowcount)
colcount = int(colcount)
overlapx = int(overlapx)
overlapy = int(overlapy)
if rowcount <= 0:
raise ValueError('rowcount(%s) should satisfy rowcount > 0' % rowcount)
if colcount <= 0:
raise ValueError('colcount(%s) should satisfy colcount > 0' % colcount)
if overlapx < 0:
raise ValueError('overlapx(%s) should satisfy overlapx >= 0' % overlapx)
if overlapy < 0:
raise ValueError('overlapy(%s) should satisfy overlapy >= 0' % overlapy)
if boundary_effect not in self._TILE_BOUNDARY_EFFECTS:
raise ValueError('boundary_effect(%s) should be one of %s' % (
boundary_effect, self._TILE_BOUNDARY_EFFECTS
))
if boundary_effect_locus not in self._TILE_BOUNDARY_EFFECT_LOCI:
raise ValueError('boundary_effect_locus(%s) should be one of %s' % (
boundary_effect_locus, self._TILE_BOUNDARY_EFFECT_LOCI
))
sizex_float = (self.rsizex + overlapx * (rowcount - 1)) / rowcount
sizey_float = (self.rsizey + overlapy * (colcount - 1)) / colcount
if boundary_effect in ['extend', 'overlap', 'shrink']:
sizex = int(np.ceil(sizex_float))
sizey = int(np.ceil(sizey_float))
elif boundary_effect == 'exclude':
sizex = int(np.floor(sizex_float))
sizey = int(np.floor(sizey_float))
elif boundary_effect == 'exception':
sizex = int(np.floor(sizex_float))
if sizex != sizex_float:
gap = int((sizex_float - sizex) * rowcount)
raise ValueError(
('There is a gap of %d pixel in the x direction, ' +
'`gap:%d %% (sizex:%d - overlapx:%d) == 0` was required') % (
(gap, gap, sizex, overlapx)))
sizey = int(np.floor(sizey_float))
if sizey != sizey_float:
gap = int((sizey_float - sizey) * colcount)
raise ValueError(
('There is a gap of %d pixel in the y direction, ' +
'`gap:%d %% (sizey:%d - overlapy:%d) == 0` was required') % (
(gap, gap, sizey, overlapy)))
else:
assert False # pragma: no cover
if sizex <= overlapx:
raise ValueError('rowcount(%d) with overlapx(%d) would not fit in %d pixels' % (
rowcount, overlapx, self.rw,
))
if sizey <= overlapy:
raise ValueError('colcount(%d) with overlapy(%d) would not fit in %d pixels' % (
colcount, overlapy, self.rw,
))
outsidex = sizex + (rowcount - 1) * (sizex - overlapx) - self.rsizex
if outsidex >= (sizex - overlapx):
raise ValueError('rowcount(%d) with overlapx(%d) would not fit in %d pixels' % (
rowcount, overlapx, self.rw,
))
outsidey = sizey + (colcount - 1) * (sizey - overlapy) - self.rsizey
if outsidey >= (sizey - overlapy):
raise ValueError('colcount(%d) with overlapy(%d) would not fit in %d pixels' % (
colcount, overlapy, self.rw,
))
size = np.asarray((sizex, sizey), dtype=int)
tiles = self._tile_unsafe(size, overlapx, overlapy, boundary_effect, boundary_effect_locus)
if boundary_effect == 'exclude':
if boundary_effect_locus == 'br':
tiles = tiles[0:colcount, 0:rowcount]
elif boundary_effect_locus == 'tl':
tiles = tiles[-colcount:, -rowcount:]
elif boundary_effect_locus == 'tr':
tiles = tiles[-colcount:, 0:rowcount]
elif boundary_effect_locus == 'bl':
tiles = tiles[0:colcount, -rowcount:]
else:
assert False # pragma: no cover
return tiles
[docs] def tile_occurrence(self, size, pixel_occurrencex, pixel_occurrencey,
boundary_effect='extend', boundary_effect_locus='br'):
"""Tile a Footprint to a matrix of Footprint
Each pixel occur `pixel_occurrencex * pixel_occurrencey` times overall in the output
Parameters
----------
size: (int, int)
Tile width and tile height, in pixel
pixel_occurrencex: int
Number of occurence of each pixel in a line of tile
pixel_occurrencey: int
Number of occurence of each pixel in a column of tile
boundary_effect: {'extend', 'exclude', 'overlap', 'shrink', 'exception'}
Behevior at boundary effect locus
- 'extend'
- Preserve tile size
- Preserve overlapx and overlapy
- Sacrifice global bounds, \
results in tiles partially outside bounds at locus (if necessary)
- Preserve tile count
- Preserve boundary pixels coverage
- 'overlap'
- Preserve tile size
- Sacrifice overlapx and overlapy \
results in tiles overlapping more at locus (if necessary)
- Preserve global bounds
- Preserve tile count
- Preserve boundary pixels coverage
- 'exclude'
- Preserve tile size
- Preserve overlapx and overlapy
- Preserve global bounds
- Sacrifice tile count, results in tiles excluded at locus (if necessary)
- Sacrifice boundary pixels coverage at locus (if necessary)
- 'shrink'
- Sacrifice tile size, results in tiles shrinked at locus (if necessary)
- Preserve overlapx and overlapy
- Preserve global bounds
- Preserve tile count
- Preserve boundary pixels coverage
- 'exception'
Raise an exception if tiles at locus do not lie inside the global bounds
boundary_effect_locus: {'br', 'tr', 'tl', 'bl'}
Locus of the boundary effects
- 'br' : Boundary effect occurs at the bottom right corner of the raster \
top left coordinates are preserved
- 'tr' : Boundary effect occurs at the top right corner of the raster, \
bottom left coordinates are preserved
- 'tl' : Boundary effect occurs at the top left corner of the raster, \
bottom right coordinates are preserved
- 'bl' : Boundary effect occurs at the bottom left corner of the raster, \
top right coordinates are preserved
Returns
-------
np.ndarray
- of dtype=object (Footpr
int)
- of shape (M, N)
- with M the line count
- with N the column count
"""
size = np.asarray(size, dtype=int)
pixel_occurrencex = int(pixel_occurrencex)
pixel_occurrencey = int(pixel_occurrencey)
if size.shape != (2,):
raise ValueError('size.shape(%s) should be (2,)' % str(size.shape))
if (size <= 0).any():
raise ValueError('size(%s) values should satisfy value > 0' % str(tuple(size)))
if not pixel_occurrencex > 0:
raise ValueError(
'pixel_occurrencex should satisfy pixel_occurrencex(%s) > 0' % pixel_occurrencex
)
elif not size[0] % pixel_occurrencex == 0:
raise ValueError(
'pixel_occurrencex should satisty size[0](%s) % pixel_occurrencex(%s) == 0'
)
if not pixel_occurrencey > 0:
raise ValueError(
'pixel_occurrencey should satisfy pixel_occurrencey(%s) > 0' % pixel_occurrencey
)
elif not size[1] % pixel_occurrencey == 0:
raise ValueError(
'pixel_occurrencey should satisty size[1](%s) % pixel_occurrencey(%s) == 0'
)
if boundary_effect not in self._TILE_OCCURRENCE_BOUNDARY_EFFECTS:
raise ValueError('boundary_effect(%s) should be one of %s' % (
boundary_effect, self._TILE_OCCURRENCE_BOUNDARY_EFFECTS
))
if boundary_effect_locus not in self._TILE_BOUNDARY_EFFECT_LOCI:
raise ValueError('boundary_effect_locus(%s) should be one of %s' % (
boundary_effect_locus, self._TILE_BOUNDARY_EFFECT_LOCI
))
occurrence = np.asarray([pixel_occurrencex, pixel_occurrencey], dtype=int)
stride = (size / occurrence).astype(int, copy=False)
overlap = size - stride
big_tl = self.tl - self.pxvec * overlap
big_rsize = self.rsize + np.asarray(overlap) * 2
big_size = big_rsize * self.pxsize
big_fp = self.__class__(tl=big_tl, size=big_size, rsize=big_rsize)
tiles = big_fp._tile_unsafe(
size, overlap[0], overlap[1], boundary_effect, boundary_effect_locus
)
return tiles
# Serialization ***************************************************************************** **
[docs] def __str__(self):
if self.angle == 0 and (self.scale * (1, -1) > 0).all():
tup = tuple(np.r_[
self.tl, self.br, self.size, self.rsize,
])
return "Footprint(tl=(%f, %f), br=(%f, %f), size=(%f, %f), rsize=(%d, %d))" % tup
else:
tup = tuple(np.r_[
self.tl, self.scale, self.angle, self.rsize,
])
return "Footprint(tl=(%f, %f), scale=(%f, %f), angle=%f, rsize=(%d, %d))" % tup
[docs] def __repr__(self):
return "Footprint(gt=%s, rsize=(%d, %d))" % (tuple(self.gt), self.rsize[0], self.rsize[1])
[docs] def __hash__(self):
# TODO: Speed test and optimize
return hash((
self._aff.to_gdal(),
tuple(self._rsize.tolist()),
))
# The end *********************************************************************************** **
# ******************************************************************************************* **
def _line_iterator(obj):
if isinstance(obj, (sg.LineString)):
yield obj
elif isinstance(obj, (sg.MultiLineString)):
for obj2 in obj.geoms:
yield obj2
elif isinstance(obj, (sg.Polygon)):
yield sg.LineString(obj.exterior)
for obj2 in obj.interiors:
yield sg.LineString(obj2)
elif isinstance(obj, (sg.MultiPolygon)):
for obj2 in obj.geoms:
yield sg.LineString(obj2.exterior)
for obj3 in obj2.interiors:
yield sg.LineString(obj3)
else:
try:
tup = tuple(obj)
except TypeError:
raise TypeError('Could not use type %s' % type(obj))
else:
for obj2 in tup:
for line in _line_iterator(obj2):
yield line
def _poly_iterator(obj):
if isinstance(obj, (sg.Polygon)):
yield obj
elif isinstance(obj, (sg.MultiPolygon)):
for obj2 in obj.geoms:
yield obj2
else:
try:
tup = tuple(obj)
except TypeError:
raise TypeError('Could not use type %s' % type(obj))
else:
for obj2 in tup:
for poly in _poly_iterator(obj2):
yield poly
def _restore(gt, rsize):
return Footprint(gt=gt, rsize=rsize)
def _angle_between(a, b, c):
return np.arccos(np.dot(
(a - b) / np.linalg.norm(a - b),
(c - b) / np.linalg.norm(c - b),
)) / np.pi * 180.