"""
Spatial extents classes for map layer and space time datasets
Usage:
.. code-block:: python
>>> import grass.temporal as tgis
>>> tgis.init()
>>> extent = tgis.RasterSpatialExtent(
... ident="raster@PERMANENT",
... north=90,
... south=90,
... east=180,
... west=180,
... top=100,
... bottom=-20,
... )
>>> extent = tgis.Raster3DSpatialExtent(
... ident="raster3d@PERMANENT",
... north=90,
... south=90,
... east=180,
... west=180,
... top=100,
... bottom=-20,
... )
>>> extent = tgis.VectorSpatialExtent(
... ident="vector@PERMANENT",
... north=90,
... south=90,
... east=180,
... west=180,
... top=100,
... bottom=-20,
... )
>>> extent = tgis.STRDSSpatialExtent(
... ident="strds@PERMANENT",
... north=90,
... south=90,
... east=180,
... west=180,
... top=100,
... bottom=-20,
... )
>>> extent = tgis.STR3DSSpatialExtent(
... ident="str3ds@PERMANENT",
... north=90,
... south=90,
... east=180,
... west=180,
... top=100,
... bottom=-20,
... )
>>> extent = tgis.STVDSSpatialExtent(
... ident="stvds@PERMANENT",
... north=90,
... south=90,
... east=180,
... west=180,
... top=100,
... bottom=-20,
... )
(C) 2012-2013 by the GRASS Development Team
This program is free software under the GNU General Public
License (>=v2). Read the file COPYING that comes with GRASS
for details.
:authors: Soeren Gebbert
"""
from .base import SQLDatabaseInterface
[docs]class SpatialExtent(SQLDatabaseInterface):
"""This is the spatial extent base class for all maps and space time datasets
This class implements a three dimensional axis aligned bounding box
and functions to compute topological relationships
Usage:
.. code-block:: python
>>> init()
>>> extent = SpatialExtent(
... table="raster_spatial_extent",
... ident="soil@PERMANENT",
... north=90,
... south=90,
... east=180,
... west=180,
... top=100,
... bottom=-20,
... )
>>> extent.id
'soil@PERMANENT'
>>> extent.north
90.0
>>> extent.south
90.0
>>> extent.east
180.0
>>> extent.west
180.0
>>> extent.top
100.0
>>> extent.bottom
-20.0
>>> extent.print_info()
+-------------------- Spatial extent ----------------------------------------+
| North:...................... 90.0
| South:...................... 90.0
| East:.. .................... 180.0
| West:....................... 180.0
| Top:........................ 100.0
| Bottom:..................... -20.0
>>> extent.print_shell_info()
north=90.0
south=90.0
east=180.0
west=180.0
top=100.0
bottom=-20.0
"""
def __init__(
self,
table=None,
ident=None,
north=None,
south=None,
east=None,
west=None,
top=None,
bottom=None,
proj="XY",
):
SQLDatabaseInterface.__init__(self, table, ident)
self.set_id(ident)
self.set_spatial_extent_from_values(north, south, east, west, top, bottom)
self.set_projection(proj)
[docs] def overlapping_2d(self, extent):
"""Return True if this (A) and the provided spatial extent (B) overlaps
in two dimensional space.
Code is lend from wind_overlap.c in lib/gis
Overlapping includes the spatial relations:
- contain
- in
- cover
- covered
- equivalent
.. code-block:: python
>>> A = SpatialExtent(north=80, south=20, east=60, west=10)
>>> B = SpatialExtent(north=80, south=20, east=60, west=10)
>>> A.overlapping_2d(B)
True
:param extent: The spatial extent to check overlapping with
:return: True or False
"""
if self.get_projection() != extent.get_projection():
self.msgr.error(
_(
"Projections are different. Unable to compute "
"overlapping_2d for spatial extents"
)
)
return False
N = extent.get_north()
S = extent.get_south()
E = extent.get_east()
W = extent.get_west()
# Adjust the east and west in case of LL projection
if self.get_projection() == "LL":
while E < self.get_west():
E += 360.0
W += 360.0
while W > self.get_east():
E -= 360.0
W -= 360.0
if self.get_north() <= S:
return False
if self.get_south() >= N:
return False
if self.get_east() <= W:
return False
if self.get_west() >= E:
return False
return True
[docs] def overlapping(self, extent):
"""Return True if this (A) and the provided spatial
extent (B) overlaps in three dimensional space.
Overlapping includes the spatial relations:
- contain
- in
- cover
- covered
- equivalent
Usage:
.. code-block:: python
>>> A = SpatialExtent(
... north=80, south=20, east=60, west=10, bottom=-50, top=50
... )
>>> B = SpatialExtent(
... north=80, south=20, east=60, west=10, bottom=-50, top=50
... )
>>> A.overlapping(B)
True
:param extent: The spatial extent to check overlapping with
:return: True or False
"""
if not self.overlapping_2d(extent):
return False
T = extent.get_top()
B = extent.get_bottom()
if self.get_top() <= B:
return False
if self.get_bottom() >= T:
return False
return True
[docs] def intersect_2d(self, extent):
"""Return the two dimensional intersection as spatial_extent
object or None in case no intersection was found.
:param extent: The spatial extent to intersect with
:return: The intersection spatial extent
"""
if not self.overlapping_2d(extent):
return None
eN = extent.get_north()
eS = extent.get_south()
eE = extent.get_east()
eW = extent.get_west()
N = self.get_north()
S = self.get_south()
E = self.get_east()
W = self.get_west()
# Adjust the east and west in case of LL projection
if self.get_projection() == "LL":
while eE < W:
eE += 360.0
eW += 360.0
while eW > E:
eE -= 360.0
eW -= 360.0
# Compute the extent
nN = N
nS = S
nE = E
nW = W
if W < eW:
nW = eW
if E > eE:
nE = eE
if N > eN:
nN = eN
if S < eS:
nS = eS
new = SpatialExtent(
north=nN,
south=nS,
east=nE,
west=nW,
top=0,
bottom=0,
proj=self.get_projection(),
)
return new
[docs] def intersect(self, extent):
"""Return the three dimensional intersection as spatial_extent
object or None in case no intersection was found.
Usage:
.. code-block:: python
>>> A = SpatialExtent(
... north=80, south=20, east=60, west=10, bottom=-50, top=50
... )
>>> B = SpatialExtent(
... north=80, south=20, east=60, west=10, bottom=-50, top=50
... )
>>> C = A.intersect(B)
>>> C.print_info()
+-------------------- Spatial extent ----------------------------------------+
| North:...................... 80.0
| South:...................... 20.0
| East:.. .................... 60.0
| West:....................... 10.0
| Top:........................ 50.0
| Bottom:..................... -50.0
>>> B = SpatialExtent(
... north=40, south=30, east=60, west=10, bottom=-50, top=50
... )
>>> C = A.intersect(B)
>>> C.print_info()
+-------------------- Spatial extent ----------------------------------------+
| North:...................... 40.0
| South:...................... 30.0
| East:.. .................... 60.0
| West:....................... 10.0
| Top:........................ 50.0
| Bottom:..................... -50.0
>>> B = SpatialExtent(
... north=40, south=30, east=60, west=30, bottom=-50, top=50
... )
>>> C = A.intersect(B)
>>> C.print_info()
+-------------------- Spatial extent ----------------------------------------+
| North:...................... 40.0
| South:...................... 30.0
| East:.. .................... 60.0
| West:....................... 30.0
| Top:........................ 50.0
| Bottom:..................... -50.0
>>> B = SpatialExtent(
... north=40, south=30, east=60, west=30, bottom=-30, top=50
... )
>>> C = A.intersect(B)
>>> C.print_info()
+-------------------- Spatial extent ----------------------------------------+
| North:...................... 40.0
| South:...................... 30.0
| East:.. .................... 60.0
| West:....................... 30.0
| Top:........................ 50.0
| Bottom:..................... -30.0
>>> B = SpatialExtent(
... north=40, south=30, east=60, west=30, bottom=-30, top=30
... )
>>> C = A.intersect(B)
>>> C.print_info()
+-------------------- Spatial extent ----------------------------------------+
| North:...................... 40.0
| South:...................... 30.0
| East:.. .................... 60.0
| West:....................... 30.0
| Top:........................ 30.0
| Bottom:..................... -30.0
:param extent: The spatial extent to intersect with
:return: The intersection spatial extent
""" # noqa: E501
if not self.overlapping(extent):
return None
new = self.intersect_2d(extent)
eT = extent.get_top()
eB = extent.get_bottom()
T = self.get_top()
B = self.get_bottom()
nT = T
nB = B
if B < eB:
nB = eB
if T > eT:
nT = eT
new.set_top(nT)
new.set_bottom(nB)
return new
[docs] def union_2d(self, extent):
"""Return the two dimensional union as spatial_extent
object or None in case the extents does not overlap or meet.
:param extent: The spatial extent to create a union with
:return: The union spatial extent
"""
if not self.overlapping_2d(extent) and not self.meet_2d(extent):
return None
return self.disjoint_union_2d(extent)
[docs] def disjoint_union_2d(self, extent):
"""Return the two dimensional union as spatial_extent.
:param extent: The spatial extent to create a union with
:return: The union spatial extent
"""
eN = extent.get_north()
eS = extent.get_south()
eE = extent.get_east()
eW = extent.get_west()
N = self.get_north()
S = self.get_south()
E = self.get_east()
W = self.get_west()
# Adjust the east and west in case of LL projection
if self.get_projection() == "LL":
while eE < W:
eE += 360.0
eW += 360.0
while eW > E:
eE -= 360.0
eW -= 360.0
# Compute the extent
nN = N
nS = S
nE = E
nW = W
if W > eW:
nW = eW
if E < eE:
nE = eE
if N < eN:
nN = eN
if S > eS:
nS = eS
new = SpatialExtent(
north=nN,
south=nS,
east=nE,
west=nW,
top=0,
bottom=0,
proj=self.get_projection(),
)
return new
[docs] def union(self, extent):
"""Return the three dimensional union as spatial_extent
object or None in case the extents does not overlap or meet.
:param extent: The spatial extent to create a union with
:return: The union spatial extent
"""
if not self.overlapping(extent) and not self.meet(extent):
return None
return self.disjoint_union(extent)
[docs] def disjoint_union(self, extent):
"""Return the three dimensional union as spatial_extent .
Usage:
.. code-block:: python
>>> A = SpatialExtent(
... north=80, south=20, east=60, west=10, bottom=-50, top=50
... )
>>> B = SpatialExtent(
... north=80, south=20, east=60, west=10, bottom=-50, top=50
... )
>>> C = A.disjoint_union(B)
>>> C.print_info()
+-------------------- Spatial extent ----------------------------------------+
| North:...................... 80.0
| South:...................... 20.0
| East:.. .................... 60.0
| West:....................... 10.0
| Top:........................ 50.0
| Bottom:..................... -50.0
>>> B = SpatialExtent(
... north=40, south=30, east=60, west=10, bottom=-50, top=50
... )
>>> C = A.disjoint_union(B)
>>> C.print_info()
+-------------------- Spatial extent ----------------------------------------+
| North:...................... 80.0
| South:...................... 20.0
| East:.. .................... 60.0
| West:....................... 10.0
| Top:........................ 50.0
| Bottom:..................... -50.0
>>> B = SpatialExtent(
... north=40, south=30, east=60, west=30, bottom=-50, top=50
... )
>>> C = A.disjoint_union(B)
>>> C.print_info()
+-------------------- Spatial extent ----------------------------------------+
| North:...................... 80.0
| South:...................... 20.0
| East:.. .................... 60.0
| West:....................... 10.0
| Top:........................ 50.0
| Bottom:..................... -50.0
>>> B = SpatialExtent(
... north=40, south=30, east=60, west=30, bottom=-30, top=50
... )
>>> C = A.disjoint_union(B)
>>> C.print_info()
+-------------------- Spatial extent ----------------------------------------+
| North:...................... 80.0
| South:...................... 20.0
| East:.. .................... 60.0
| West:....................... 10.0
| Top:........................ 50.0
| Bottom:..................... -50.0
>>> B = SpatialExtent(
... north=40, south=30, east=60, west=30, bottom=-30, top=30
... )
>>> C = A.disjoint_union(B)
>>> C.print_info()
+-------------------- Spatial extent ----------------------------------------+
| North:...................... 80.0
| South:...................... 20.0
| East:.. .................... 60.0
| West:....................... 10.0
| Top:........................ 50.0
| Bottom:..................... -50.0
>>> A = SpatialExtent(
... north=80, south=20, east=60, west=10, bottom=-50, top=50
... )
>>> B = SpatialExtent(
... north=90, south=80, east=70, west=20, bottom=-30, top=60
... )
>>> C = A.disjoint_union(B)
>>> C.print_info()
+-------------------- Spatial extent ----------------------------------------+
| North:...................... 90.0
| South:...................... 20.0
| East:.. .................... 70.0
| West:....................... 10.0
| Top:........................ 60.0
| Bottom:..................... -50.0
:param extent: The spatial extent to create a disjoint union with
:return: The union spatial extent
""" # noqa: E501
new = self.disjoint_union_2d(extent)
eT = extent.get_top()
eB = extent.get_bottom()
T = self.get_top()
B = self.get_bottom()
nT = T
nB = B
if B > eB:
nB = eB
if T < eT:
nT = eT
new.set_top(nT)
new.set_bottom(nB)
return new
[docs] def is_in_2d(self, extent):
"""Return True if this extent (A) is located in the provided spatial
extent (B) in two dimensions.
::
_____
|A _ |
| |_| |
|_____|B
:param extent: The spatial extent
:return: True or False
"""
if self.get_projection() != extent.get_projection():
self.msgr.error(
_(
"Projections are different. Unable to compute "
"is_in_2d for spatial extents"
)
)
return False
eN = extent.get_north()
eS = extent.get_south()
eE = extent.get_east()
eW = extent.get_west()
N = self.get_north()
S = self.get_south()
E = self.get_east()
W = self.get_west()
# Adjust the east and west in case of LL projection
if self.get_projection() == "LL":
while eE < W:
eE += 360.0
eW += 360.0
while eW > E:
eE -= 360.0
eW -= 360.0
if W <= eW:
return False
if E >= eE:
return False
if N >= eN:
return False
if S <= eS:
return False
return True
[docs] def is_in(self, extent):
"""Return True if this extent (A) is located in the provided spatial
extent (B) in three dimensions.
Usage:
.. code-block:: python
>>> A = SpatialExtent(
... north=79, south=21, east=59, west=11, bottom=-49, top=49
... )
>>> B = SpatialExtent(
... north=80, south=20, east=60, west=10, bottom=-50, top=50
... )
>>> A.is_in(B)
True
>>> B.is_in(A)
False
:param extent: The spatial extent
:return: True or False
"""
if not self.is_in_2d(extent):
return False
eT = extent.get_top()
eB = extent.get_bottom()
T = self.get_top()
B = self.get_bottom()
if B <= eB:
return False
if T >= eT:
return False
return True
[docs] def contain_2d(self, extent):
"""Return True if this extent (A) contains the provided spatial
extent (B) in two dimensions.
Usage:
.. code-block:: python
>>> A = SpatialExtent(north=80, south=20, east=60, west=10)
>>> B = SpatialExtent(north=79, south=21, east=59, west=11)
>>> A.contain_2d(B)
True
>>> B.contain_2d(A)
False
:param extent: The spatial extent
:return: True or False
"""
return extent.is_in_2d(self)
[docs] def contain(self, extent):
"""Return True if this extent (A) contains the provided spatial
extent (B) in three dimensions.
Usage:
.. code-block:: python
>>> A = SpatialExtent(
... north=80, south=20, east=60, west=10, bottom=-50, top=50
... )
>>> B = SpatialExtent(
... north=79, south=21, east=59, west=11, bottom=-49, top=49
... )
>>> A.contain(B)
True
>>> B.contain(A)
False
:param extent: The spatial extent
:return: True or False
"""
return extent.is_in(self)
[docs] def equivalent_2d(self, extent):
"""Return True if this extent (A) is equal to the provided spatial
extent (B) in two dimensions.
Usage:
.. code-block:: python
>>> A = SpatialExtent(north=80, south=20, east=60, west=10)
>>> B = SpatialExtent(north=80, south=20, east=60, west=10)
>>> A.equivalent_2d(B)
True
>>> B.equivalent_2d(A)
True
:param extent: The spatial extent
:return: True or False
"""
if self.get_projection() != extent.get_projection():
self.msgr.error(
_(
"Projections are different. Unable to compute "
"equivalent_2d for spatial extents"
)
)
return False
eN = extent.get_north()
eS = extent.get_south()
eE = extent.get_east()
eW = extent.get_west()
N = self.get_north()
S = self.get_south()
E = self.get_east()
W = self.get_west()
# Adjust the east and west in case of LL projection
if self.get_projection() == "LL":
while eE < W:
eE += 360.0
eW += 360.0
while eW > E:
eE -= 360.0
eW -= 360.0
if W != eW:
return False
if E != eE:
return False
if N != eN:
return False
if S != eS:
return False
return True
[docs] def equivalent(self, extent):
"""Return True if this extent (A) is equal to the provided spatial
extent (B) in three dimensions.
Usage:
.. code-block:: python
>>> A = SpatialExtent(
... north=80, south=20, east=60, west=10, bottom=-50, top=50
... )
>>> B = SpatialExtent(
... north=80, south=20, east=60, west=10, bottom=-50, top=50
... )
>>> A.equivalent(B)
True
>>> B.equivalent(A)
True
:param extent: The spatial extent
:return: True or False
"""
if not self.equivalent_2d(extent):
return False
eT = extent.get_top()
eB = extent.get_bottom()
T = self.get_top()
B = self.get_bottom()
if B != eB:
return False
if T != eT:
return False
return True
[docs] def cover_2d(self, extent):
"""Return True if this extent (A) covers the provided spatial
extent (B) in two dimensions.
::
_____ _____ _____ _____
|A __| |__ A| |A | B| |B | A|
| |B | | B| | | |__| |__| |
|__|__| |__|__| |_____| |_____|
_____ _____ _____ _____
|A|B| | |A __| |A _ | |__ A|
| |_| | | |__|B | |B| | B|__| |
|_____| |_____| |_|_|_| |_____|
_____ _____ _____ _____
|A|B | |_____|A |A|B|A| |_____|A
| | | |B | | | | | |_____|B
|_|___| |_____| |_|_|_| |_____|A
The following cases are excluded:
- contain
- in
- equivalent
:param extent: The spatial extent
:return: True or False
"""
if self.get_projection() != extent.get_projection():
self.msgr.error(
_(
"Projections are different. Unable to compute"
" cover_2d for spatial extents"
)
)
return False
# Exclude equivalent_2d
if self.equivalent_2d(extent):
return False
eN = extent.get_north()
eS = extent.get_south()
eE = extent.get_east()
eW = extent.get_west()
N = self.get_north()
S = self.get_south()
E = self.get_east()
W = self.get_west()
# Adjust the east and west in case of LL projection
if self.get_projection() == "LL":
while eE < W:
eE += 360.0
eW += 360.0
while eW > E:
eE -= 360.0
eW -= 360.0
# Edges of extent located outside of self are not allowed
if E <= eW:
return False
if W >= eE:
return False
if N <= eS:
return False
if S >= eN:
return False
# First we check that at least one edge of extent meets an edge of self
if W != eW and E != eE and N != eN and S != eS:
return False
# We check that at least one edge of extent is located in self
edge_count = 0
if W < eW and E > eW:
edge_count += 1
if E > eE and W < eE:
edge_count += 1
if N > eN and S < eN:
edge_count += 1
if S < eS and N > eS:
edge_count += 1
if edge_count == 0:
return False
return True
[docs] def cover(self, extent):
"""Return True if this extent covers the provided spatial
extent in three dimensions.
The following cases are excluded:
- contain
- in
- equivalent
:param extent: The spatial extent
:return: True or False
"""
if self.get_projection() != extent.get_projection():
self.msgr.error(
_(
"Projections are different. Unable to compute "
"cover for spatial extents"
)
)
return False
# Exclude equivalent_2d
if self.equivalent_2d(extent):
return False
eN = extent.get_north()
eS = extent.get_south()
eE = extent.get_east()
eW = extent.get_west()
eT = extent.get_top()
eB = extent.get_bottom()
N = self.get_north()
S = self.get_south()
E = self.get_east()
W = self.get_west()
T = self.get_top()
B = self.get_bottom()
# Adjust the east and west in case of LL projection
if self.get_projection() == "LL":
while eE < W:
eE += 360.0
eW += 360.0
while eW > E:
eE -= 360.0
eW -= 360.0
# Edges of extent located outside of self are not allowed
if E <= eW:
return False
if W >= eE:
return False
if N <= eS:
return False
if S >= eN:
return False
if T <= eB:
return False
if B >= eT:
return False
# First we check that at least one edge of extent meets an edge of self
if W != eW and E != eE and N != eN and S != eS and B != eB and T != eT:
return False
# We check that at least one edge of extent is located in self
edge_count = 0
if W < eW and E > eW:
edge_count += 1
if E > eE and W < eE:
edge_count += 1
if N > eN and S < eN:
edge_count += 1
if S < eS and N > eS:
edge_count += 1
if N > eN and S < eN:
edge_count += 1
if S < eS and N > eS:
edge_count += 1
if T > eT and B < eT:
edge_count += 1
if B < eB and T > eB:
edge_count += 1
if edge_count == 0:
return False
return True
[docs] def covered_2d(self, extent):
"""Return True if this extent is covered by the provided spatial
extent in two dimensions.
The following cases are excluded:
- contain
- in
- equivalent
:param extent: The spatial extent
:return: True or False
"""
return extent.cover_2d(self)
[docs] def covered(self, extent):
"""Return True if this extent is covered by the provided spatial
extent in three dimensions.
The following cases are excluded:
- contain
- in
- equivalent
:param extent: The spatial extent
:return: True or False
"""
return extent.cover(self)
[docs] def overlap_2d(self, extent):
"""Return True if this extent (A) overlaps with the provided spatial
extent (B) in two dimensions.
Code is lend from wind_overlap.c in lib/gis
::
_____
|A __|__
| | | B|
|__|__| |
|_____|
The following cases are excluded:
- contain
- in
- cover
- covered
- equivalent
:param extent: The spatial extent
:return: True or False
"""
if self.contain_2d(extent):
return False
if self.is_in_2d(extent):
return False
if self.cover_2d(extent):
return False
if self.covered_2d(extent):
return False
if self.equivalent_2d(extent):
return False
N = extent.get_north()
S = extent.get_south()
E = extent.get_east()
W = extent.get_west()
# Adjust the east and west in case of LL projection
if self.get_projection() == "LL":
while E < self.get_west():
E += 360.0
W += 360.0
while W > self.get_east():
E -= 360.0
W -= 360.0
if self.get_north() <= S:
return False
if self.get_south() >= N:
return False
if self.get_east() <= W:
return False
if self.get_west() >= E:
return False
return True
[docs] def overlap(self, extent):
"""Return True if this extent overlaps with the provided spatial
extent in three dimensions.
The following cases are excluded:
- contain
- in
- cover
- covered
- equivalent
:param extent: The spatial extent
:return: True or False
"""
if self.is_in(extent):
return False
if self.contain(extent):
return False
if self.cover(extent):
return False
if self.covered(extent):
return False
if self.equivalent(extent):
return False
N = extent.get_north()
S = extent.get_south()
E = extent.get_east()
W = extent.get_west()
T = extent.get_top()
B = extent.get_bottom()
# Adjust the east and west in case of LL projection
if self.get_projection() == "LL":
while E < self.get_west():
E += 360.0
W += 360.0
while W > self.get_east():
E -= 360.0
W -= 360.0
if self.get_north() <= S:
return False
if self.get_south() >= N:
return False
if self.get_east() <= W:
return False
if self.get_west() >= E:
return False
if self.get_top() <= B:
return False
if self.get_bottom() >= T:
return False
return True
[docs] def meet_2d(self, extent):
"""Return True if this extent (A) meets with the provided spatial
extent (B) in two dimensions.
::
_____ _____
| A | B |
|_____| |
|_____|
_____ _____
| B | A |
| | |
|_____|_____|
___
| A |
| |
|___|
| B |
| |
|_____|
_____
| B |
| |
|_____|_
| A |
| |
|_____|
:param extent: The spatial extent
:return: True or False
"""
eN = extent.get_north()
eS = extent.get_south()
eE = extent.get_east()
eW = extent.get_west()
N = self.get_north()
S = self.get_south()
E = self.get_east()
W = self.get_west()
# Adjust the east and west in case of LL projection
if self.get_projection() == "LL":
while eE < W:
eE += 360.0
eW += 360.0
while eW > E:
eE -= 360.0
eW -= 360.0
edge = None
edge_count = 0
if E == eW:
edge = "E"
edge_count += 1
if W == eE:
edge = "W"
edge_count += 1
if N == eS:
edge = "N"
edge_count += 1
if S == eN:
edge = "S"
edge_count += 1
# Meet a a single edge only
if edge_count != 1:
return False
# Check boundaries of the faces
if edge == "E" or edge == "W":
if N < eS or S > eN:
return False
if edge == "N" or edge == "S":
if E < eW or W > eE:
return False
return True
[docs] def meet(self, extent):
"""Return True if this extent meets with the provided spatial
extent in three dimensions.
:param extent: The spatial extent
:return: True or False
"""
eN = extent.get_north()
eS = extent.get_south()
eE = extent.get_east()
eW = extent.get_west()
eT = extent.get_top()
eB = extent.get_bottom()
N = self.get_north()
S = self.get_south()
E = self.get_east()
W = self.get_west()
T = self.get_top()
B = self.get_bottom()
# Adjust the east and west in case of LL projection
if self.get_projection() == "LL":
while eE < W:
eE += 360.0
eW += 360.0
while eW > E:
eE -= 360.0
eW -= 360.0
edge = None
edge_count = 0
if E == eW:
edge = "E"
edge_count += 1
if W == eE:
edge = "W"
edge_count += 1
if N == eS:
edge = "N"
edge_count += 1
if S == eN:
edge = "S"
edge_count += 1
if T == eB:
edge = "T"
edge_count += 1
if B == eT:
edge = "B"
edge_count += 1
# Meet a single edge only
if edge_count != 1:
return False
# Check boundaries of the faces
if edge == "E" or edge == "W":
if N < eS or S > eN:
return False
if T < eB or B > eT:
return False
if edge == "N" or edge == "S":
if E < eW or W > eE:
return False
if T < eB or B > eT:
return False
if edge == "T" or edge == "B":
if E < eW or W > eE:
return False
if N < eS or S > eN:
return False
return True
[docs] def disjoint_2d(self, extent):
"""Return True if this extent (A) is disjoint with the provided spatial
extent (B) in three dimensions.
::
_____
| A |
|_____|
_______
| B |
|_______|
:param extent: The spatial extent
:return: True or False
"""
if self.is_in_2d(extent):
return False
if self.contain_2d(extent):
return False
if self.cover_2d(extent):
return False
if self.covered_2d(extent):
return False
if self.equivalent_2d(extent):
return False
if self.overlapping_2d(extent):
return False
if self.meet_2d(extent):
return False
return True
[docs] def disjoint(self, extent):
"""Return True if this extent is disjoint with the provided spatial
extent in three dimensions.
:param extent: The spatial extent
:return: True or False
"""
if self.is_in(extent):
return False
if self.contain(extent):
return False
if self.cover(extent):
return False
if self.covered(extent):
return False
if self.equivalent(extent):
return False
if self.overlapping(extent):
return False
if self.meet(extent):
return False
return True
[docs] def spatial_relation_2d(self, extent):
"""Returns the two dimensional spatial relation between this
extent and the provided spatial extent in two dimensions.
Spatial relations are:
- disjoint
- meet
- overlap
- cover
- covered
- in
- contain
- equivalent
Usage: see self.spatial_relation()
"""
if self.equivalent_2d(extent):
return "equivalent"
if self.contain_2d(extent):
return "contain"
if self.is_in_2d(extent):
return "in"
if self.cover_2d(extent):
return "cover"
if self.covered_2d(extent):
return "covered"
if self.overlap_2d(extent):
return "overlap"
if self.meet_2d(extent):
return "meet"
if self.disjoint_2d(extent):
return "disjoint"
return "unknown"
[docs] def spatial_relation(self, extent):
"""Returns the two dimensional spatial relation between this
extent and the provided spatial extent in three dimensions.
Spatial relations are:
- disjoint
- meet
- overlap
- cover
- covered
- in
- contain
- equivalent
Usage:
.. code-block:: python
>>> A = SpatialExtent(
... north=80, south=20, east=60, west=10, bottom=-50, top=50
... )
>>> B = SpatialExtent(
... north=80, south=20, east=60, west=10, bottom=-50, top=50
... )
>>> A.spatial_relation(B)
'equivalent'
>>> B.spatial_relation(A)
'equivalent'
>>> B = SpatialExtent(
... north=70, south=20, east=60, west=10, bottom=-50, top=50
... )
>>> A.spatial_relation_2d(B)
'cover'
>>> A.spatial_relation(B)
'cover'
>>> B = SpatialExtent(
... north=70, south=30, east=60, west=10, bottom=-50, top=50
... )
>>> A.spatial_relation_2d(B)
'cover'
>>> A.spatial_relation(B)
'cover'
>>> B.spatial_relation_2d(A)
'covered'
>>> B.spatial_relation(A)
'covered'
>>> B = SpatialExtent(
... north=70, south=30, east=50, west=10, bottom=-50, top=50
... )
>>> A.spatial_relation_2d(B)
'cover'
>>> B.spatial_relation_2d(A)
'covered'
>>> A.spatial_relation(B)
'cover'
>>> B = SpatialExtent(
... north=70, south=30, east=50, west=20, bottom=-50, top=50
... )
>>> B.spatial_relation(A)
'covered'
>>> B = SpatialExtent(
... north=70, south=30, east=50, west=20, bottom=-50, top=50
... )
>>> A.spatial_relation_2d(B)
'contain'
>>> A.spatial_relation(B)
'cover'
>>> B = SpatialExtent(
... north=70, south=30, east=50, west=20, bottom=-40, top=50
... )
>>> A.spatial_relation(B)
'cover'
>>> B = SpatialExtent(
... north=70, south=30, east=50, west=20, bottom=-40, top=40
... )
>>> A.spatial_relation(B)
'contain'
>>> B.spatial_relation(A)
'in'
>>> B = SpatialExtent(
... north=90, south=30, east=50, west=20, bottom=-40, top=40
... )
>>> A.spatial_relation_2d(B)
'overlap'
>>> A.spatial_relation(B)
'overlap'
>>> B = SpatialExtent(
... north=90, south=5, east=70, west=5, bottom=-40, top=40
... )
>>> A.spatial_relation_2d(B)
'in'
>>> A.spatial_relation(B)
'overlap'
>>> B = SpatialExtent(
... north=90, south=5, east=70, west=5, bottom=-40, top=60
... )
>>> A.spatial_relation(B)
'overlap'
>>> B = SpatialExtent(
... north=90, south=5, east=70, west=5, bottom=-60, top=60
... )
>>> A.spatial_relation(B)
'in'
>>> A = SpatialExtent(
... north=80, south=60, east=60, west=10, bottom=-50, top=50
... )
>>> B = SpatialExtent(
... north=60, south=20, east=60, west=10, bottom=-50, top=50
... )
>>> A.spatial_relation_2d(B)
'meet'
>>> A.spatial_relation(B)
'meet'
>>> A = SpatialExtent(
... north=60, south=40, east=60, west=10, bottom=-50, top=50
... )
>>> B = SpatialExtent(
... north=80, south=60, east=60, west=10, bottom=-50, top=50
... )
>>> A.spatial_relation_2d(B)
'meet'
>>> A.spatial_relation(B)
'meet'
>>> A = SpatialExtent(
... north=80, south=40, east=60, west=40, bottom=-50, top=50
... )
>>> B = SpatialExtent(
... north=80, south=40, east=40, west=20, bottom=-50, top=50
... )
>>> A.spatial_relation_2d(B)
'meet'
>>> A.spatial_relation(B)
'meet'
>>> A = SpatialExtent(
... north=80, south=40, east=40, west=20, bottom=-50, top=50
... )
>>> B = SpatialExtent(
... north=90, south=30, east=60, west=40, bottom=-50, top=50
... )
>>> A.spatial_relation_2d(B)
'meet'
>>> A.spatial_relation(B)
'meet'
>>> A = SpatialExtent(
... north=80, south=40, east=40, west=20, bottom=-50, top=50
... )
>>> B = SpatialExtent(
... north=70, south=50, east=60, west=40, bottom=-50, top=50
... )
>>> A.spatial_relation_2d(B)
'meet'
>>> A.spatial_relation(B)
'meet'
>>> A = SpatialExtent(
... north=80, south=40, east=40, west=20, bottom=-50, top=50
... )
>>> B = SpatialExtent(
... north=60, south=20, east=60, west=40, bottom=-50, top=50
... )
>>> A.spatial_relation_2d(B)
'meet'
>>> A.spatial_relation(B)
'meet'
>>> A = SpatialExtent(
... north=80, south=40, east=40, west=20, bottom=-50, top=50
... )
>>> B = SpatialExtent(
... north=40, south=20, east=60, west=40, bottom=-50, top=50
... )
>>> A.spatial_relation_2d(B)
'disjoint'
>>> A.spatial_relation(B)
'disjoint'
>>> A = SpatialExtent(
... north=80, south=40, east=40, west=20, bottom=-50, top=50
... )
>>> B = SpatialExtent(
... north=60, south=20, east=60, west=40, bottom=-60, top=60
... )
>>> A.spatial_relation(B)
'meet'
>>> A = SpatialExtent(
... north=80, south=40, east=40, west=20, bottom=-50, top=50
... )
>>> B = SpatialExtent(
... north=90, south=30, east=60, west=40, bottom=-40, top=40
... )
>>> A.spatial_relation(B)
'meet'
>>> A = SpatialExtent(
... north=80, south=40, east=60, west=20, bottom=0, top=50
... )
>>> B = SpatialExtent(
... north=80, south=40, east=60, west=20, bottom=-50, top=0
... )
>>> A.spatial_relation(B)
'meet'
>>> A = SpatialExtent(
... north=80, south=40, east=60, west=20, bottom=0, top=50
... )
>>> B = SpatialExtent(
... north=80, south=50, east=60, west=30, bottom=-50, top=0
... )
>>> A.spatial_relation(B)
'meet'
>>> A = SpatialExtent(
... north=80, south=40, east=60, west=20, bottom=0, top=50
... )
>>> B = SpatialExtent(
... north=70, south=50, east=50, west=30, bottom=-50, top=0
... )
>>> A.spatial_relation(B)
'meet'
>>> A = SpatialExtent(
... north=80, south=40, east=60, west=20, bottom=0, top=50
... )
>>> B = SpatialExtent(
... north=90, south=30, east=70, west=10, bottom=-50, top=0
... )
>>> A.spatial_relation(B)
'meet'
>>> A = SpatialExtent(
... north=80, south=40, east=60, west=20, bottom=0, top=50
... )
>>> B = SpatialExtent(
... north=70, south=30, east=50, west=10, bottom=-50, top=0
... )
>>> A.spatial_relation(B)
'meet'
>>> A = SpatialExtent(
... north=80, south=40, east=60, west=20, bottom=-50, top=0
... )
>>> B = SpatialExtent(
... north=80, south=40, east=60, west=20, bottom=0, top=50
... )
>>> A.spatial_relation(B)
'meet'
>>> A = SpatialExtent(
... north=80, south=40, east=60, west=20, bottom=-50, top=0
... )
>>> B = SpatialExtent(
... north=80, south=50, east=60, west=30, bottom=0, top=50
... )
>>> A.spatial_relation(B)
'meet'
>>> A = SpatialExtent(
... north=80, south=40, east=60, west=20, bottom=-50, top=0
... )
>>> B = SpatialExtent(
... north=70, south=50, east=50, west=30, bottom=0, top=50
... )
>>> A.spatial_relation(B)
'meet'
>>> A = SpatialExtent(
... north=80, south=40, east=60, west=20, bottom=-50, top=0
... )
>>> B = SpatialExtent(
... north=90, south=30, east=70, west=10, bottom=0, top=50
... )
>>> A.spatial_relation(B)
'meet'
>>> A = SpatialExtent(
... north=80, south=40, east=60, west=20, bottom=-50, top=0
... )
>>> B = SpatialExtent(
... north=70, south=30, east=50, west=10, bottom=0, top=50
... )
>>> A.spatial_relation(B)
'meet'
>>> A = SpatialExtent(
... north=80, south=20, east=60, west=10, bottom=-50, top=50
... )
>>> B = SpatialExtent(
... north=90, south=81, east=60, west=10, bottom=-50, top=50
... )
>>> A.spatial_relation(B)
'disjoint'
>>> A = SpatialExtent(
... north=80, south=20, east=60, west=10, bottom=-50, top=50
... )
>>> B = SpatialExtent(
... north=90, south=80, east=60, west=10, bottom=-50, top=50
... )
>>> A.spatial_relation(B)
'meet'
"""
if self.equivalent(extent):
return "equivalent"
if self.contain(extent):
return "contain"
if self.is_in(extent):
return "in"
if self.cover(extent):
return "cover"
if self.covered(extent):
return "covered"
if self.overlap(extent):
return "overlap"
if self.meet(extent):
return "meet"
if self.disjoint(extent):
return "disjoint"
return "unknown"
[docs] def set_spatial_extent_from_values(self, north, south, east, west, top, bottom):
"""Set the three dimensional spatial extent
:param north: The northern edge
:param south: The southern edge
:param east: The eastern edge
:param west: The western edge
:param top: The top edge
:param bottom: The bottom edge
"""
self.set_north(north)
self.set_south(south)
self.set_east(east)
self.set_west(west)
self.set_top(top)
self.set_bottom(bottom)
[docs] def set_spatial_extent(self, spatial_extent):
"""Set the three dimensional spatial extent
:param spatial_extent: An object of type SpatialExtent or its
subclasses
"""
self.set_north(spatial_extent.get_north())
self.set_south(spatial_extent.get_south())
self.set_east(spatial_extent.get_east())
self.set_west(spatial_extent.get_west())
self.set_top(spatial_extent.get_top())
self.set_bottom(spatial_extent.get_bottom())
[docs] def set_projection(self, proj):
"""Set the projection of the spatial extent it should be XY or LL.
As default the projection is XY
"""
if proj is None or (proj != "XY" and proj != "LL"):
self.D["proj"] = "XY"
else:
self.D["proj"] = proj
[docs] def set_spatial_extent_from_values_2d(self, north, south, east, west):
"""Set the two dimensional spatial extent from values
:param north: The northern edge
:param south: The southern edge
:param east: The eastern edge
:param west: The western edge
"""
self.set_north(north)
self.set_south(south)
self.set_east(east)
self.set_west(west)
[docs] def set_spatial_extent_2d(self, spatial_extent):
"""Set the three dimensional spatial extent
:param spatial_extent: An object of type SpatialExtent or its
subclasses
"""
self.set_north(spatial_extent.north)
self.set_south(spatial_extent.south)
self.set_east(spatial_extent.east)
self.set_west(spatial_extent.west)
[docs] def set_id(self, ident):
"""Convenient method to set the unique identifier (primary key)"""
self.ident = ident
self.D["id"] = ident
[docs] def set_north(self, north):
"""Set the northern edge of the map"""
if north is not None:
self.D["north"] = float(north)
else:
self.D["north"] = None
[docs] def set_south(self, south):
"""Set the southern edge of the map"""
if south is not None:
self.D["south"] = float(south)
else:
self.D["south"] = None
[docs] def set_west(self, west):
"""Set the western edge of the map"""
if west is not None:
self.D["west"] = float(west)
else:
self.D["west"] = None
[docs] def set_east(self, east):
"""Set the eastern edge of the map"""
if east is not None:
self.D["east"] = float(east)
else:
self.D["east"] = None
[docs] def set_top(self, top):
"""Set the top edge of the map"""
if top is not None:
self.D["top"] = float(top)
else:
self.D["top"] = None
[docs] def set_bottom(self, bottom):
"""Set the bottom edge of the map"""
if bottom is not None:
self.D["bottom"] = float(bottom)
else:
self.D["bottom"] = None
[docs] def get_id(self):
"""Convenient method to get the unique identifier (primary key)
:return: None if not found
"""
if "id" in self.D:
return self.D["id"]
else:
return None
[docs] def get_projection(self):
"""Get the projection of the spatial extent"""
return self.D["proj"]
[docs] def get_volume(self):
"""Compute the volume of the extent, in case z is zero
(top == bottom or top - bottom = 1) the area is returned"""
if self.get_projection() == "LL":
self.msgr.error(
_("Volume computation is not supported " "for LL projections")
)
area = self.get_area()
bbox = self.get_spatial_extent_as_tuple()
z = abs(bbox[4] - bbox[5])
if z == 0:
z = 1.0
return area * z
[docs] def get_area(self):
"""Compute the area of the extent, extent in z direction is ignored"""
if self.get_projection() == "LL":
self.msgr.error(
_("Area computation is not supported " "for LL projections")
)
bbox = self.get_spatial_extent_as_tuple()
y = abs(bbox[0] - bbox[1])
x = abs(bbox[2] - bbox[3])
return x * y
[docs] def get_spatial_extent_as_tuple(self):
"""Return a tuple (north, south, east, west, top, bottom)
of the spatial extent"""
return (self.north, self.south, self.east, self.west, self.top, self.bottom)
[docs] def get_spatial_extent_as_tuple_2d(self):
"""Return a tuple (north, south, east, west,) of the 2d spatial extent"""
return (self.north, self.south, self.east, self.west)
[docs] def get_north(self):
"""Get the northern edge of the map
:return: None if not found"""
if "north" in self.D:
return self.D["north"]
else:
return None
[docs] def get_south(self):
"""Get the southern edge of the map
:return: None if not found"""
if "south" in self.D:
return self.D["south"]
else:
return None
[docs] def get_east(self):
"""Get the eastern edge of the map
:return: None if not found"""
if "east" in self.D:
return self.D["east"]
else:
return None
[docs] def get_west(self):
"""Get the western edge of the map
:return: None if not found"""
if "west" in self.D:
return self.D["west"]
else:
return None
[docs] def get_top(self):
"""Get the top edge of the map
:return: None if not found"""
if "top" in self.D:
return self.D["top"]
else:
return None
[docs] def get_bottom(self):
"""Get the bottom edge of the map
:return: None if not found"""
if "bottom" in self.D:
return self.D["bottom"]
else:
return None
id = property(fget=get_id, fset=set_id)
north = property(fget=get_north, fset=set_north)
south = property(fget=get_south, fset=set_south)
east = property(fget=get_east, fset=set_east)
west = property(fget=get_west, fset=set_west)
top = property(fget=get_top, fset=set_top)
bottom = property(fget=get_bottom, fset=set_bottom)
[docs] def print_info(self):
"""Print information about this class in human readable style"""
# 0123456789012345678901234567890
print(
" +-------------------- Spatial extent ----------------------------------------+" # noqa: E501
)
print(" | North:...................... " + str(self.get_north()))
print(" | South:...................... " + str(self.get_south()))
print(" | East:.. .................... " + str(self.get_east()))
print(" | West:....................... " + str(self.get_west()))
print(" | Top:........................ " + str(self.get_top()))
print(" | Bottom:..................... " + str(self.get_bottom()))
[docs] def print_shell_info(self):
"""Print information about this class in shell style"""
print("north=" + str(self.get_north()))
print("south=" + str(self.get_south()))
print("east=" + str(self.get_east()))
print("west=" + str(self.get_west()))
print("top=" + str(self.get_top()))
print("bottom=" + str(self.get_bottom()))
###############################################################################
[docs]class RasterSpatialExtent(SpatialExtent):
def __init__(
self,
ident=None,
north=None,
south=None,
east=None,
west=None,
top=None,
bottom=None,
):
SpatialExtent.__init__(
self, "raster_spatial_extent", ident, north, south, east, west, top, bottom
)
[docs]class Raster3DSpatialExtent(SpatialExtent):
def __init__(
self,
ident=None,
north=None,
south=None,
east=None,
west=None,
top=None,
bottom=None,
):
SpatialExtent.__init__(
self,
"raster3d_spatial_extent",
ident,
north,
south,
east,
west,
top,
bottom,
)
[docs]class VectorSpatialExtent(SpatialExtent):
def __init__(
self,
ident=None,
north=None,
south=None,
east=None,
west=None,
top=None,
bottom=None,
):
SpatialExtent.__init__(
self, "vector_spatial_extent", ident, north, south, east, west, top, bottom
)
[docs]class STRDSSpatialExtent(SpatialExtent):
def __init__(
self,
ident=None,
north=None,
south=None,
east=None,
west=None,
top=None,
bottom=None,
):
SpatialExtent.__init__(
self, "strds_spatial_extent", ident, north, south, east, west, top, bottom
)
[docs]class STR3DSSpatialExtent(SpatialExtent):
def __init__(
self,
ident=None,
north=None,
south=None,
east=None,
west=None,
top=None,
bottom=None,
):
SpatialExtent.__init__(
self, "str3ds_spatial_extent", ident, north, south, east, west, top, bottom
)
[docs]class STVDSSpatialExtent(SpatialExtent):
def __init__(
self,
ident=None,
north=None,
south=None,
east=None,
west=None,
top=None,
bottom=None,
):
SpatialExtent.__init__(
self, "stvds_spatial_extent", ident, north, south, east, west, top, bottom
)
###############################################################################
if __name__ == "__main__":
import doctest
doctest.testmod()
