Skip to content

Commit

Permalink
introduced new point to tile functionality;
Browse files Browse the repository at this point in the history
  • Loading branch information
cnavacch committed Aug 21, 2022
1 parent 8dbebbc commit a07b3a7
Show file tree
Hide file tree
Showing 2 changed files with 114 additions and 52 deletions.
84 changes: 48 additions & 36 deletions src/geospade/raster.py
Original file line number Diff line number Diff line change
Expand Up @@ -880,7 +880,7 @@ def xy2rc(self, x, y, sref=None, px_origin=None) -> Tuple[int, int]:
"""

if sref is not None:
x, y = transform_coords(x, y, sref, self.sref.osr_sref)
x, y = transform_coords(x, y, sref, self.sref)
px_origin = self.px_origin if px_origin is None else px_origin
c, r = xy2ij(x, y, self.geotrans, origin=px_origin)
return r, c
Expand Down Expand Up @@ -1841,18 +1841,59 @@ def xy2tile(self, x, y, sref=None) -> "Tile":
osr_sref = self.sref.osr_sref if sref is None else sref.osr_sref
point.AssignSpatialReference(osr_sref)

return self.poi2tile(point, sref=sref)

@_align_geom(align=True)
def poi2tile(self, poi, sref=None) -> Tile:
"""
Returns the tile intersecting with the given point of interest in world system coordinates. If the coordinates
are outside the mosaic boundary, no tile is returned.
Parameters
----------
poi : ogr.wkbPoint
Point of interest.
sref : SpatialRef, optional
Spatial reference system of the world system coordinates. If None, the spatial reference system of the
coordinates and the spatial reference system of the mosaic are assumed to be the same (default).
Returns
-------
geospade.raster.Tile :
Tile intersecting/matching with the given world system coordinates.
"""

tile_oi = None
if point.Intersects(self.boundary): # point is inside grid definition
if self._poi_intersects(poi): # point is inside grid definition
for tile in self.all_tiles:
if tile.intersects(point): # point is inside tile
if tile.intersects(poi): # point is inside tile
tile_oi = self._mask_tile(tile)
break
else:
wrn_msg = f"The given coordinate tuple ({x},{y}) is not within the mosaic boundary."
warnings.warn(wrn_msg)

return tile_oi

def _poi_intersects(self, poi) -> bool:
"""
Checks if the given point intersects with the mosaic. If not, a warning is raised.
Parameters
----------
poi : ogr.wkbPoint
Point of interest.
Returns
-------
poi_intscts : bool
True if the point intersects with the mosaic, false if not.
"""
poi_intscts = poi.Intersects(self.boundary)
if not poi_intscts:
wrn_msg = f"The given point ({poi.GetX()},{poi.GetY()}) is not within the mosaic boundary."
warnings.warn(wrn_msg)
return poi_intscts

def name2tile(self, tile_name, active_only=True, apply_mask=True) -> "Tile":
"""
Converts a tile name to a tile.
Expand Down Expand Up @@ -2586,35 +2627,6 @@ def shape(self) -> Tuple[int, int]:
"""
return self._adjacency_matrix.shape

def xy2tile(self, x, y, sref=None) -> Tile:
"""
Returns the tile intersecting with the given world system coordinates. If the coordinates are outside the
mosaic boundary, no tile is returned.
Parameters
----------
x : number
World system coordinate in X direction.
y : number
World system coordinate in Y direction.
sref : SpatialRef, optional
Spatial reference system of the world system coordinates. If None, the spatial reference system of the
coordinates and the spatial reference system of the mosaic are assumed to be the same (default).
Returns
-------
geospade.raster.Tile :
Tile intersecting/matching with the given world system coordinates.
"""
# create OGR point
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(x, y)
osr_sref = self.sref.osr_sref if sref is None else sref.osr_sref
point.AssignSpatialReference(osr_sref)

return self.poi2tile(point, sref=sref)

@_align_geom(align=True)
def poi2tile(self, poi, sref=None) -> Tile:
"""
Expand All @@ -2637,7 +2649,7 @@ def poi2tile(self, poi, sref=None) -> Tile:
"""

tile_oi = None
if poi.Intersects(self.boundary): # point is inside grid definition
if self._poi_intersects(poi): # point is inside grid definition
mosaic_row, mosaic_col = self._raster_geom.xy2rc(poi.GetX(), poi.GetY(), sref=sref)
tile_oi = self.__tile_from_rc(mosaic_row, mosaic_col)

Expand Down
82 changes: 66 additions & 16 deletions src/geospade/tools.py
Original file line number Diff line number Diff line change
Expand Up @@ -327,46 +327,96 @@ def any_geom2ogr_geom(geom, sref=None) -> ogr.Geometry:
return ogr_geom


def _round_geom_coords(geom, decimals) -> ogr.Geometry:
def _round_polygon_coords(polygon, decimals) -> ogr.wkbPolygon:
"""
'Cleans' the coordinates, so that it has rounded coordinates.
'Cleans' the coordinates of an OGR polygon, so that it has rounded coordinates.
Parameters
----------
geom : ogr.Geometry
An OGR geometry.
polygon : ogr.wkbPolygon
An OGR polygon.
decimals : int
Number of significant digits to round to.
Returns
-------
geometry_out : ogr.Geometry
An OGR geometry with 'cleaned' and rounded coordinates.
geometry_out : ogr.wkbPolygon
An OGR polygon with rounded coordinates.
"""
if not geom.ExportToWkt().startswith('POLYGON'):
err_msg = "Only OGR Polygons are supported at the moment."
raise NotImplementedError(err_msg)

ring = geom.GetGeometryRef(0)
ring = polygon.GetGeometryRef(0)

rounded_ring = ogr.Geometry(ogr.wkbLinearRing)

n_points = ring.GetPointCount()

for p in range(n_points):
lon, lat, z = ring.GetPoint(p)
rlon, rlat, rz = [np.round(lon, decimals=decimals),
np.round(lat, decimals=decimals),
np.round(z, decimals=decimals)]
rounded_ring.AddPoint(rlon, rlat, rz)
x, y, z = ring.GetPoint(p)
rx, ry, rz = [np.round(x, decimals=decimals),
np.round(y, decimals=decimals),
np.round(z, decimals=decimals)]
rounded_ring.AddPoint(rx, ry, rz)

geometry_out = ogr.Geometry(ogr.wkbPolygon)
geometry_out.AddGeometry(rounded_ring)

return geometry_out


def _round_point_coords(point, decimals) -> ogr.wkbPoint:
"""
'Cleans' the coordinates of an OGR polygon, so that it has rounded coordinates.
Parameters
----------
point : ogr.wkbPoint
An OGR point.
decimals : int
Number of significant digits to round to.
Returns
-------
geometry_out : ogr.wkbPoint
An OGR point with rounded coordinates.
"""
rx, ry = np.round(point.GetX(), decimals=decimals), np.round(point.GetY(), decimals=decimals)
rpoint = ogr.Geometry(ogr.wkbPoint)
rpoint.AddPoint(rx, ry)

return rpoint


def _round_geom_coords(geom, decimals) -> ogr.Geometry:
"""
'Cleans' the coordinates of an OGR geometry, so that it has rounded coordinates.
Parameters
----------
geom : ogr.Geometry
An OGR geometry.
decimals : int
Number of significant digits to round to.
Returns
-------
geometry_out : ogr.Geometry
An OGR geometry with 'cleaned' and rounded coordinates.
Notes
-----
Only OGR polygons and points are supported.
"""
if geom.ExportToWkt().startswith('POLYGON'):
return _round_polygon_coords(geom, decimals)
elif geom.ExportToWkt().startswith('POINT'):
return _round_point_coords(geom, decimals)
else:
err_msg = "Only OGR Polygons and Points are supported at the moment."
raise NotImplementedError(err_msg)


def rasterise_polygon(geom, x_pixel_size, y_pixel_size, extent=None) -> np.ndarray:
"""
Rasterises a Shapely polygon defined by a clockwise list of points.
Expand Down

0 comments on commit a07b3a7

Please sign in to comment.