"""!@package grass.temporal @brief GRASS Python scripting module (temporal GIS functions) Temporal GIS related spatial extent functions to be used in Python scripts and tgis packages. Usage: @code import grass.temporal as tgis extent = tgis.raster_spatial_extent() ... @endcode (C) 2008-2011 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. @author Soeren Gebbert """ from base import * class spatial_extent(sql_database_interface): """!This is the spatial extent base class for all maps and space time datasets""" def __init__(self, table=None, ident=None, north=None, south=None, east=None, west=None, top=None, bottom=None, proj="XY"): sql_database_interface.__init__(self, table, ident) self.set_id(ident) self.set_spatial_extent(north, south, east, west, top, bottom) self.set_projection(proj) def overlapping_2d(self, extent): """!Return True if the two dimensional extents overlap. Code is lend from wind_overlap.c in lib/gis Overlapping includes the spatial relations: * contain * in * cover * covered * equivalent """ if self.get_projection() != extent.get_projection(): core.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 def overlapping(self, extent): """!Return True if the three dimensional extents overlap Overlapping includes the spatial relations: * contain * in * cover * covered * equivalent """ 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 def intersect_2d(self, extent): """!Return the two dimensional intersection as spatial_extent object or None in case no intersection was found. """ 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 = spatial_extent(north=nN, south=nS, east=nE, west=nW, top=0, bottom=0, proj=self.get_projection()) return new def intersect(self, extent): """!Return the three dimensional intersection as spatial_extent object or None in case no intersection was found. """ 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 def is_in_2d(self, extent): """Check two dimensional if the self is located in extent _____ |A _ | | |_| | |_____| """ if self.get_projection() != extent.get_projection(): core.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 def is_in(self, extent): """Check three dimensional if the self is located in extent """ 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 def contain_2d(self, extent): """Check two dimensional if self contains extent """ return extent.is_in_2d(self) def contain(self, extent): """Check three dimensional if self contains extent """ return extent.is_in(self) def equivalent_2d(self, extent): """Check two dimensional if self is equivalent to extent """ if self.get_projection() != extent.get_projection(): core.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 def equivalent(self, extent): """Check three dimensional if self is equivalent to extent """ 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 def cover_2d(self, extent): """Return True if two dimensional self covers extent _____ _____ _____ _____ |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 """ if self.get_projection() != extent.get_projection(): core.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 def cover(self, extent): """Return True if three dimensional self covers extent The following cases are excluded: * contain * in * equivalent """ if self.get_projection() != extent.get_projection(): core.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 def covered_2d(self, extent): """Check two dimensional if self is covered by extent """ return extent.cover_2d(self) def covered(self, extent): """Check three dimensional if self is covered by extent """ return extent.cover(self) def overlap_2d(self, extent): """Return True if the two dimensional extents overlap. Code is lend from wind_overlap.c in lib/gis _____ |A __|__ | | | B| |__|__| | |_____| The following cases are excluded: * contain * in * cover * covered * equivalent """ 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 def overlap(self, extent): """Return True if the three dimensional extents overlap The following cases are excluded: * contain * in * cover * covered * equivalent """ 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 def meet_2d(self,extent): """ Check if self and extent meet each other in two dimensions _____ _____ _____ _____ | A | B | | B | A | |_____| | | | | |_____| |_____|_____| ___ | A | | | |___| _____ | B | | B | | | | | |_____| |_____|_ | A | | | |_____| """ 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() # 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 def meet(self,extent): """ Check if self and extent meet each other in three dimensions""" 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 def disjoint_2d(self, extent): """Return True if the two dimensional extents are disjoint """ if self.overlapping_2d(extent) or self.meet_2d(extent): return False return True def disjoint(self, extent): """Return True if the three dimensional extents are disjoint """ if self.overlapping(extent) or self.meet(extent): return False return True def spatial_relation_2d(self, extent): """Returns the two dimensional spatial relation between self and extent Spatial relations are: * disjoint * meet * overlap * cover * covered * in * contain * equivalent """ 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" def spatial_relation(self, extent): """Returns the three dimensional spatial relation between self and extent Spatial relations are: * disjoint * meet * overlap * cover * covered * in * contain * equivalent """ 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" def set_spatial_extent(self, north, south, east, west, top, bottom): """Set the spatial extent""" self.set_north(north) self.set_south(south) self.set_east(east) self.set_west(west) self.set_top(top) self.set_bottom(bottom) 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 == None or (proj != "XY" and proj != "LL"): self.D["proj"] = "XY" else: self.D["proj"] = proj def set_spatial_extent_2d(self, north, south, east, west): self.set_id(ident) self.set_north(north) self.set_south(south) self.set_east(east) self.set_west(west) self.set_top(0) self.set_bottom(0) def set_id(self, ident): """Convenient method to set the unique identifier (primary key)""" self.ident = ident self.D["id"] = ident def set_north(self, north): """Set the northern edge of the map""" self.D["north"] = north def set_south(self, sourth): """Set the southern edge of the map""" self.D["south"] = sourth def set_west(self, west): """Set the western edge of the map""" self.D["west"] = west def set_east(self, east): """Set the eastern edge of the map""" self.D["east"] = east def set_top(self, top): """Set the top edge of the map""" self.D["top"] = top def set_bottom(self, bottom): """Set the bottom edge of the map""" self.D["bottom"] = bottom def get_id(self): """Convenient method to get the unique identifier (primary key) @return None if not found """ if self.D.has_key("id"): return self.D["id"] else: return None def get_projection(self): """Get the projection of the spatial extent""" return self.D["proj"] 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": core.error(_("Volume computation is not supported for LL projections")) area = self.get_area() bbox = self.get_spatial_extent() z = abs(bbox[4] - bbox[5]) if z == 0: z = 1.0 return area*z def get_area(self): """Compute the area of the extent, extent in z direction is ignored""" if self.get_projection() == "LL": core.error(_("Area computation is not supported for LL projections")) bbox = self.get_spatial_extent() y = abs(bbox[0] - bbox[1]) x = abs(bbox[2] - bbox[3]) return x*y def get_spatial_extent(self): """Return a tuple (north, south, east, west, top, bottom) of the spatial extent""" return (self.get_north(), self.get_south, self.get_east(), self.get_west(), \ self.get_top(), self.get_bottom()) def get_spatial_extent_2d(self): """Return a tuple (north, south, east, west,) of the 2d spatial extent""" return (self.get_north(), self.get_south, self.get_east(), self.get_west()) def get_north(self): """Get the northern edge of the map @return None if not found""" if self.D.has_key("north"): return self.D["north"] else: return None def get_south(self): """Get the southern edge of the map @return None if not found""" if self.D.has_key("south"): return self.D["south"] else: return None def get_east(self): """Get the eastern edge of the map @return None if not found""" if self.D.has_key("east"): return self.D["east"] else: return None def get_west(self): """Get the western edge of the map @return None if not found""" if self.D.has_key("west"): return self.D["west"] else: return None def get_top(self): """Get the top edge of the map @return None if not found""" if self.D.has_key("top"): return self.D["top"] else: return None def get_bottom(self): """Get the bottom edge of the map @return None if not found""" if self.D.has_key("bottom"): return self.D["bottom"] else: return None def print_info(self): """Print information about this class in human readable style""" # 0123456789012345678901234567890 print " +-------------------- Spatial extent ----------------------------------------+" 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()) 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()) ############################################################################### class raster_spatial_extent(spatial_extent): def __init__(self, ident=None, north=None, south=None, east=None, west=None, top=None, bottom=None): spatial_extent.__init__(self, "raster_spatial_extent", ident, north, south, east, west, top, bottom) class raster3d_spatial_extent(spatial_extent): def __init__(self, ident=None, north=None, south=None, east=None, west=None, top=None, bottom=None): spatial_extent.__init__(self, "raster3d_spatial_extent", ident, north, south, east, west, top, bottom) class vector_spatial_extent(spatial_extent): def __init__(self, ident=None, north=None, south=None, east=None, west=None, top=None, bottom=None): spatial_extent.__init__(self, "vector_spatial_extent", ident, north, south, east, west, top, bottom) class strds_spatial_extent(spatial_extent): def __init__(self, ident=None, north=None, south=None, east=None, west=None, top=None, bottom=None): spatial_extent.__init__(self, "strds_spatial_extent", ident, north, south, east, west, top, bottom) class str3ds_spatial_extent(spatial_extent): def __init__(self, ident=None, north=None, south=None, east=None, west=None, top=None, bottom=None): spatial_extent.__init__(self, "str3ds_spatial_extent", ident, north, south, east, west, top, bottom) class stvds_spatial_extent(spatial_extent): def __init__(self, ident=None, north=None, south=None, east=None, west=None, top=None, bottom=None): spatial_extent.__init__(self, "stvds_spatial_extent", ident, north, south, east, west, top, bottom)