added a lot of functions, checked importability, refactored map module to 'geo'

ptds/__init__.py

+#import py_tools_ds.ptds.io
+#import py_tools_ds.ptds.geo
+#import py_tools_ds.ptds.numeric
+#import py_tools_ds.ptds.processing
 from . import io
-from . import map
+from . import geo
 from . import numeric
 from . import processing
 
+from .io.raster.GeoArray import GeoArray
+#from py_tools_ds.ptds.io.raster.GeoArray import GeoArray
+
 __all__=['io',
          'map',
          'numeric',
-         'processing']
+         'processing',
+         'GeoArray']
 
 __version__ = '20160927.01'
 __author__='Daniel Scheffler'
\ No newline at end of file

ptds/map/coord_calc.py → ptds/geo/coord_calc.py

@@ -7,8 +7,6 @@ import numpy as np
 # custom
 from shapely.geometry import Polygon
 
-from .vector import topology as TOPO
-
 
 
 def get_corner_coordinates(gdal_ds=None, gt=None, cols=None, rows=None):
@@ -106,10 +104,11 @@ def calc_FullDataset_corner_positions(mask_1bit, assert_four_corners=True, algor
                           (first_dataCol, StartStopRows_in_first_dataCol[0]))
                 line_O = ((last_dataCol, StartStopRows_in_last_dataCol[1]),
                           (StartStopCols_in_last_dataRow[1], last_dataRow))
-                corners = [list(reversed(TOPO.find_line_intersection_point(line_N, line_W))),
-                           list(reversed(TOPO.find_line_intersection_point(line_N, line_O))),
-                           list(reversed(TOPO.find_line_intersection_point(line_S, line_W))),
-                           list(reversed(TOPO.find_line_intersection_point(line_S, line_O)))]
+                from .vector.topology import find_line_intersection_point # import here avoids circular dependencies
+                corners = [list(reversed(find_line_intersection_point(line_N, line_W))),
+                           list(reversed(find_line_intersection_point(line_N, line_O))),
+                           list(reversed(find_line_intersection_point(line_S, line_W))),
+                           list(reversed(find_line_intersection_point(line_S, line_O)))]
             else:
                 dataRow_in_first_dataCol = np.mean(StartStopRows_in_first_dataCol)
                 dataRow_in_last_dataCol  = np.mean(StartStopRows_in_last_dataCol)
@@ -124,4 +123,16 @@ def calc_FullDataset_corner_positions(mask_1bit, assert_four_corners=True, algor
             # if UL[0] == 0 and UR[0] == 0:
             #     envi.save_image(os.path.abspath('./testing/out/__faulty_mask_1bit.hdr'), self.mask_1bit, dtype= 'uint8', interleave='bsq', ext='.bsq',force=True)
             # print(UL,UR,LL,LR)
-    return out
\ No newline at end of file
+    return out
+
+
+def corner_coord_to_minmax(corner_coords):
+    """Returns the bounding coordinates for a given set of XY coordinates.
+
+    :param corner_coords:   # four XY tuples of corner coordinates. Their order does not matter.
+    :return: xmin,xmax,ymin,ymax
+    """
+    x_vals = [i[0] for i in corner_coords]
+    y_vals = [i[1] for i in corner_coords]
+    xmin, xmax, ymin, ymax = min(x_vals), max(x_vals), min(y_vals), max(y_vals)
+    return xmin, xmax, ymin, ymax
\ No newline at end of file

ptds/map/coord_grid.py → ptds/geo/coord_grid.py

@@ -3,8 +3,10 @@ __author__ = "Daniel Scheffler"
 
 
 import numpy as np
+from shapely.geometry import box
 
 from ..numeric.vector import find_nearest
+from .coord_calc      import get_corner_coordinates
 
 
 def snap_bounds_to_pixGrid(bounds, gt, roundAlg='auto'):
@@ -40,3 +42,27 @@ def is_point_on_grid(pointXY,xgrid,ygrid):
         return False
 
 
+def find_nearest_grid_coord(valXY,gt,rows,cols,direction='NW',extrapolate=True):
+    UL, LL, LR, UR = get_corner_coordinates(gt=gt, rows=rows, cols=cols)  # (x,y) tuples
+    round_x = {'NW': 'off', 'NO': 'on', 'SW': 'off', 'SO': 'on'}[direction]
+    round_y = {'NW': 'on', 'NO': 'on', 'SW': 'off', 'SO': 'off'}[direction]
+    tgt_xgrid = np.arange(UL[0], UR[0] + gt[1], gt[1])
+    tgt_ygrid = np.arange(LL[1], UL[1] + abs(gt[5]), abs(gt[5]))
+    tgt_x = find_nearest(tgt_xgrid, valXY[0], roundAlg=round_x, extrapolate=extrapolate)
+    tgt_y = find_nearest(tgt_ygrid, valXY[1], roundAlg=round_y, extrapolate=extrapolate)
+    return tgt_x,tgt_y
+
+
+def move_shapelyPoly_to_image_grid(shapelyPoly, gt, rows, cols, moving_dir='NW'):
+    polyULxy = (min(shapelyPoly.exterior.coords.xy[0]), max(shapelyPoly.exterior.coords.xy[1]))
+    polyLRxy = (max(shapelyPoly.exterior.coords.xy[0]), min(shapelyPoly.exterior.coords.xy[1]))
+    UL, LL, LR, UR = get_corner_coordinates(gt=gt, rows=rows, cols=cols)  # (x,y) tuples
+    round_x = {'NW': 'off', 'NO': 'on', 'SW': 'off', 'SE': 'on'}[moving_dir]
+    round_y = {'NW': 'on', 'NO': 'on', 'SW': 'off', 'SE': 'off'}[moving_dir]
+    tgt_xgrid = np.arange(UL[0], UR[0] + gt[1], gt[1])
+    tgt_ygrid = np.arange(LL[1], UL[1] + abs(gt[5]), abs(gt[5]))
+    tgt_xmin = find_nearest(tgt_xgrid, polyULxy[0], roundAlg=round_x)
+    tgt_xmax = find_nearest(tgt_xgrid, polyLRxy[0], roundAlg=round_x)
+    tgt_ymin = find_nearest(tgt_ygrid, polyLRxy[1], roundAlg=round_y)
+    tgt_ymax = find_nearest(tgt_ygrid, polyULxy[1], roundAlg=round_y)
+    return box(tgt_xmin, tgt_ymin, tgt_xmax, tgt_ymax)
\ No newline at end of file

ptds/map/coord_trafo.py → ptds/geo/coord_trafo.py

@@ -79,9 +79,9 @@ def transform_any_prj(prj_src,prj_tgt,x,y):
 
 def mapXY2imXY(mapXY, gt):
     # type: (tuple,list) -> tuple
-    """Translates given map coordinates into pixel locations according to the given image geotransform.
+    """Translates given geo coordinates into pixel locations according to the given image geotransform.
 
-    :param mapXY:   <tuple> The map coordinates to be translated in the form (x,y).
+    :param mapXY:   <tuple> The geo coordinates to be translated in the form (x,y).
     :param gt:      <list> GDAL geotransform
     :returns:       <tuple> image coordinate tuple X/Y (column, row)
     """
@@ -90,15 +90,47 @@ def mapXY2imXY(mapXY, gt):
 
 def imXY2mapXY(imXY, gt):
     # type: (tuple,list) -> tuple
-    """Translates given pixel locations into map coordinates according to the given image geotransform.
+    """Translates given pixel locations into geo coordinates according to the given image geotransform.
 
     :param imXY:    <tuple> The image coordinates to be translated in the form (x,y).
     :param gt:      <list> GDAL geotransform
-    :returns:       <tuple> map coordinate tuple X/Y (mapX, mapY)
+    :returns:       <tuple> geo coordinate tuple X/Y (mapX, mapY)
     """
     return (gt[0] + imXY[0] * abs(gt[1])), (gt[3] - imXY[1] * abs(gt[5]))
 
 
+mapYX2imYX = lambda mapYX,gt: ((mapYX[0]-gt[3])/gt[5],     (mapYX[1]-gt[0])/gt[1])
+imYX2mapYX = lambda imYX ,gt: (gt[3]-(imYX[0]*abs(gt[5])), gt[0]+(imYX[1]*gt[1]))
+
+
+def pixelToMapYX(pixelCoords, path_im=None, geotransform=None, projection=None):
+    assert path_im or geotransform and projection, \
+        "pixelToMapYX: Missing argument! Please provide either 'path_im' or 'geotransform' AND 'projection'."
+    if geotransform and projection:
+        gt,proj = geotransform, projection
+    else:
+        ds       = gdal.Open(path_im)
+        gt, proj = ds.GetGeoTransform(), ds.GetProjection()
+    srs = osr.SpatialReference()
+    srs.ImportFromWkt(proj)
+    # Set up the coordinate transformation object
+    ct  = osr.CoordinateTransformation(srs, srs)
+
+    mapYmapXPairs = []
+    pixelCoords   = [pixelCoords] if not type(pixelCoords[0]) in [list,tuple] else pixelCoords
+
+    for point in pixelCoords:
+        # Translate the pixel pairs into untranslated points
+        u_mapX = point[0] * gt[1] + gt[0]
+        u_mapY = point[1] * gt[5] + gt[3]
+        # Transform the points to the space
+        (mapX, mapY, holder) = ct.TransformPoint(u_mapX, u_mapY)
+        # Add the point to our return array
+        mapYmapXPairs.append([mapY, mapX])
+
+    return mapYmapXPairs
+
+
 def pixelToLatLon(pixelPairs, path_im=None, geotransform=None, projection=None):
     """The following method translates given pixel locations into latitude/longitude locations on a given GEOTIF
     This method does not take into account pixel size and assumes a high enough

ptds/map/map_info.py → ptds/geo/map_info.py

@@ -11,10 +11,10 @@ from .coord_trafo import transform_utm_to_wgs84
 
 
 def geotransform2mapinfo(gt,prj):
-    """Builds an ENVI map info from given GDAL GeoTransform and Projection (compatible with UTM and LonLat projections).
+    """Builds an ENVI geo info from given GDAL GeoTransform and Projection (compatible with UTM and LonLat projections).
     :param gt:  GDAL GeoTransform, e.g. (249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0)
     :param prj: GDAL Projection - WKT Format
-    :returns:   ENVI map info, e.g. [ UTM , 1 , 1 , 256785.0 , 4572015.0 , 30.0 , 30.0 , 43 , North , WGS-84 ]
+    :returns:   ENVI geo info, e.g. [ UTM , 1 , 1 , 256785.0 , 4572015.0 , 30.0 , 30.0 , 43 , North , WGS-84 ]
     :rtype:     list
     """
 
@@ -37,8 +37,8 @@ def geotransform2mapinfo(gt,prj):
 
 def mapinfo2geotransform(map_info):
     # type: (list) -> list
-    """Builds GDAL GeoTransform tuple from an ENVI map info.
-    :param map_info: ENVI map info (list), e.g. ['UTM', 1, 1, 192585.0, 5379315.0, 30.0, 30.0, 41, 'North', 'WGS-84']
+    """Builds GDAL GeoTransform tuple from an ENVI geo info.
+    :param map_info: ENVI geo info (list), e.g. ['UTM', 1, 1, 192585.0, 5379315.0, 30.0, 30.0, 41, 'North', 'WGS-84']
     :returns:        GDAL GeoTransform, e.g. [249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0]
     :rtype:          list
     """
@@ -49,7 +49,7 @@ def mapinfo2geotransform(map_info):
         ULmapX = float(map_info[3])-(float(map_info[1])*float(map_info[5])-float(map_info[5]))
         ULmapY = float(map_info[4])+(float(map_info[2])*float(map_info[6])-float(map_info[6]))
     assert isinstance(map_info,list) and len(map_info)==10,\
-        "The map info argument has to be a list of length 10. Got %s." %map_info
+        "The geo info argument has to be a list of length 10. Got %s." %map_info
     return [ULmapX,float(map_info[5]),0.,ULmapY,0.,-float(map_info[6])]
 
 def get_corner_coordinates(gdal_ds=None, gt=None, cols=None, rows=None):

ptds/map/projection.py → ptds/geo/projection.py

@@ -51,6 +51,7 @@ def prj_equal(prj1, prj2):
 
 def isProjectedOrGeographic(prj):
     """:param prj: accepts EPSG, Proj4 and WKT projections"""
+    if prj is None: return None
     srs = osr.SpatialReference()
     if prj.startswith('EPSG:'):
         srs.ImportFromEPSG(int(prj.split(':')[1]))
@@ -75,21 +76,24 @@ def EPSG2WKT(EPSG_code):
     return srs.ExportToWkt()
 
 
-def WKT2EPSG(wkt, epsg=os.environ['GDAL_DATA'].replace('/gdal','/proj/epsg')):
+def WKT2EPSG(wkt, epsg=os.environ['GDAL_DATA'].replace('/gdal', '/proj/epsg')):
     """ Transform a WKT string to an EPSG code
-    :param wkt:  WKT definition, e.g. the result of gdal_ds.GetProjection()
+    :param wkt:  WKT definition
     :param epsg: the proj.4 epsg file (defaults to '/usr/local/share/proj/epsg')
     :returns:    EPSG code
     http://gis.stackexchange.com/questions/20298/is-it-possible-to-get-the-epsg-value-from-an-osr-spatialreference-class-using-th
     """
+    # FIXME this function returns None if datum=NAD27 but works with datum=WGS84, e.g.:
+    # FIXME {PROJCS["UTM_Zone_33N",GEOGCS["GCS_North_American_1927",DATUM["D_North_American_1927",SPHEROID["Clarke_1866",6378206.4,294.9786982]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",15.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]}
     code = None #default
     p_in = osr.SpatialReference()
-    s = p_in.ImportFromWkt(wkt)
+    s    = p_in.ImportFromWkt(wkt)
     if s == 5:  # invalid WKT
         raise Exception('Invalid WKT.')
     if p_in.IsLocal():
         raise Exception('The given WKT is a local coordinate system.')
     cstype = 'GEOGCS' if p_in.IsGeographic() else 'PROJCS'
+    p_in.AutoIdentifyEPSG()
     an = p_in.GetAuthorityName(cstype)
     assert an in [None,'EPSG'], "No EPSG code found. Found %s instead." %an
     ac = p_in.GetAuthorityCode(cstype)
@@ -106,7 +110,7 @@ def WKT2EPSG(wkt, epsg=os.environ['GDAL_DATA'].replace('/gdal','/proj/epsg')):
                 code = code if code else None # match or no match
     else:
         code = ac
-    return int(code)
+    return code
 
 
 def get_UTMzone(ds=None,prj=None):

ptds/geo/vector/__init__.py

+# -*- coding: utf-8 -*-
+__author__ = "Daniel Scheffler"
+
+from . import geometry
+from . import topology
+from . import conversion
+
+__all__=['geometry',
+         'topology',
+         'conversion']
\ No newline at end of file

ptds/geo/vector/conversion.py

+# -*- coding: utf-8 -*-
+__author__ = "Daniel Scheffler"
+
+import numpy as np
+
+# custom
+from shapely.geometry import shape, mapping, box
+
+from ..coord_trafo import imYX2mapYX, mapYX2imYX, pixelToMapYX
+
+
+
+def shapelyImPoly_to_shapelyMapPoly_withPRJ(shapelyImPoly, gt,prj):
+    #ACTUALLY PRJ IS NOT NEEDED BUT THIS FUNCTION RETURNS OTHER VALUES THAN shapelyImPoly_to_shapelyMapPoly
+    geojson   = mapping(shapelyImPoly)
+    coords    = list(geojson['coordinates'][0])
+    coordsYX = pixelToMapYX(coords,geotransform=gt,projection=prj)
+    coordsXY = tuple([(i[1],i[0]) for i in coordsYX])
+    geojson['coordinates'] = (coordsXY,)
+    return shape(geojson)
+
+
+def shapelyImPoly_to_shapelyMapPoly(shapelyBox, gt):
+    xmin, ymin, xmax, ymax = shapelyBox.bounds
+    ymax, xmin = imYX2mapYX((ymax, xmin), gt)
+    ymin, xmax = imYX2mapYX((ymin, xmax), gt)
+    return box(xmin, ymin, xmax, ymax)
+
+
+def shapelyBox2BoxYX(shapelyBox,coord_type='image'):
+    xmin, ymin, xmax, ymax = shapelyBox.bounds
+    assert coord_type in ['image','map']
+    UL_YX, UR_YX, LR_YX, LL_YX = ((ymin,xmin),(ymin,xmax),(ymax,xmax),(ymax,xmin)) if coord_type=='image' else \
+                                 ((ymax,xmin),(ymax,xmax),(ymin,xmax),(ymin,xmin))
+    return UL_YX, UR_YX, LR_YX, LL_YX
+
+
+def get_boxImXY_from_shapelyPoly(shapelyPoly, im_gt):
+    # type: (shapely.Polygon,tuple) -> np.ndarray
+    """Converts each vertex coordinate of a shapely polygon into image coordinates corresponding to the given
+    geotransform without respect to invalid image coordinates. Those must be filtered later.
+    :param shapelyPoly:     <shapely.Polygon>
+    :param im_gt:           <list> the GDAL geotransform of the target image
+    """
+    get_coordsArr = lambda shpPoly: np.swapaxes(np.array(shpPoly.exterior.coords.xy), 0, 1)
+    coordsArr     = get_coordsArr(shapelyPoly)
+    boxImXY       = [mapYX2imYX((Y, X), im_gt) for X, Y in coordsArr.tolist()] # FIXME incompatible to GMS version
+    boxImXY       = [(i[1],i[0]) for i in boxImXY]
+    return boxImXY
+
+
+def round_shapelyPoly_coords(shapelyPoly, precision=10,out_dtype=None):
+    geojson = mapping(shapelyPoly)
+    geojson['coordinates'] = np.round(np.array(geojson['coordinates']), precision)
+    if out_dtype:
+        geojson['coordinates'] = geojson['coordinates'].astype(out_dtype)
+    return shape(geojson)
\ No newline at end of file

ptds/geo/vector/geometry.py

+# -*- coding: utf-8 -*-
+__author__ = "Daniel Scheffler"
+
+
+from shapely.geometry import Polygon, box
+
+
+from .conversion import shapelyImPoly_to_shapelyMapPoly, get_boxImXY_from_shapelyPoly, \
+                        shapelyBox2BoxYX, round_shapelyPoly_coords
+from ..coord_calc  import corner_coord_to_minmax
+from ..coord_trafo import imYX2mapYX
+
+
+class boxObj(object):
+    def __init__(self, **kwargs):
+        """Create a dynamic/self-updating box object that represents a rectangular or quadratic coordinate box
+        according to the given keyword arguments.
+        Note: Either mapPoly+gt or imPoly+gt or wp+ws or boxMapYX+gt or boxImYX+gt must be passed.
+
+        :Keyword Arguments:
+            - gt (tuple):                   GDAL geotransform (default: (0, 1, 0, 0, 0, -1))
+            - prj (str):                    projection as WKT string
+            - mapPoly (shapely.Polygon):    Polygon with map unit vertices
+            - imPoly (shapely.Polygon):     Polygon with image unit vertices
+            - wp (tuple):                   window position in map units (x,y)
+            - ws (tuple):                   window size in map units (x,y)
+            - boxMapYX (list):              box map coordinates like [(ULy,ULx), (URy,URx), (LRy,LRx), (LLy,LLx)]
+            - boxImYX (list):               box image coordinates like [(ULy,ULx), (URy,URx), (LRy,LRx), (LLy,LLx)]
+        """
+        # FIXME self.prj is not used
+        self.gt        = kwargs.get('gt',       (0, 1, 0, 0, 0, -1))
+        self.prj       = kwargs.get('prj',      '')
+        self._mapPoly  = kwargs.get('mapPoly',  None)
+        self._imPoly   = kwargs.get('imPoly',   None)
+        self.wp        = kwargs.get('wp',       None)
+        self._ws       = kwargs.get('ws',       None)
+        self._boxMapYX = kwargs.get('boxMapYX', None)
+        self._boxImYX  = kwargs.get('boxImYX',  None)
+
+        if self._mapPoly:
+            assert self.gt, "A geotransform must be passed if mapPoly is given."
+        else:
+            # populate self._mapPoly
+            if self._imPoly:
+                self._mapPoly = shapelyImPoly_to_shapelyMapPoly(self._imPoly,self.gt)
+            elif self._boxMapYX:
+                self.boxMapYX = self._boxMapYX
+            elif self._boxImYX:
+                self.boxImYX = self._boxImYX
+            elif self.wp or self._ws: # asserts wp/ws in map units
+                assert self.wp and self._ws,\
+                    "'wp' and 'ws' must be passed together. Got wp=%s and ws=%s." %(self.wp,self._ws)
+                (wpx,wpy),(wsx,wsy) = self.wp,self._ws
+                self._mapPoly = box(wpx-wsx/2, wpy-wsy/2, wpx+wsx/2, wpy+wsy/2)
+            else:
+                raise Exception("No proper set of arguments received.")
+
+    # all getters and setters synchronize using self._mapPoly
+    @property
+    def mapPoly(self):
+        imPoly = Polygon(get_boxImXY_from_shapelyPoly(self._mapPoly, self.gt))
+        return shapelyImPoly_to_shapelyMapPoly(imPoly, self.gt)
+
+    @mapPoly.setter
+    def mapPoly(self, shapelyPoly):
+        self._mapPoly = shapelyPoly
+
+    @property
+    def imPoly(self):
+        return Polygon(get_boxImXY_from_shapelyPoly(self.mapPoly, self.gt))
+
+    @imPoly.setter
+    def imPoly(self, shapelyImPoly):
+        self.mapPoly = shapelyImPoly_to_shapelyMapPoly(shapelyImPoly, self.gt)
+
+    @property
+    def boxMapYX(self):
+        return shapelyBox2BoxYX(self.mapPoly, coord_type='map')
+
+    @boxMapYX.setter
+    def boxMapYX(self, mapBoxYX):
+        mapBoxXY = [(i[1], i[0]) for i in mapBoxYX]
+        xmin, xmax, ymin, ymax = corner_coord_to_minmax(mapBoxXY)
+        self.mapPoly = box(xmin, ymin, xmax, ymax)
+
+    @property
+    def boxImYX(self):
+        temp_imPoly  = round_shapelyPoly_coords(self.imPoly, precision=0, out_dtype=int)
+        floatImBoxYX = shapelyBox2BoxYX(temp_imPoly, coord_type='image')
+        return [[int(i[0]),int(i[1])] for i in floatImBoxYX]
+
+    @boxImYX.setter
+    def boxImYX(self, imBoxYX):
+        imBoxXY = [(i[1], i[0]) for i in imBoxYX]
+        xmin, xmax, ymin, ymax = corner_coord_to_minmax(imBoxXY)
+        self.imPoly  = box(xmin, ymin, xmax, ymax)
+
+    @property
+    def boundsMap(self):
+        """Returns xmin,xmax,ymin,ymax in map coordinates."""
+        boxMapYX = shapelyBox2BoxYX(self.mapPoly, coord_type='image')
+        boxMapXY = [(i[1], i[0]) for i in boxMapYX]
+        return corner_coord_to_minmax(boxMapXY)
+
+    @property
+    def boundsIm(self):
+        """Returns xmin,xmax,ymin,ymax in image coordinates."""
+        boxImXY = get_boxImXY_from_shapelyPoly(self.mapPoly, self.gt)
+        return corner_coord_to_minmax(boxImXY) # xmin,xmax,ymin,ymax
+
+    @property
+    def imDimsYX(self):
+        xmin, xmax, ymin, ymax = self.boundsIm
+        return (ymax-ymin),(xmax-xmin)
+
+    @property
+    def mapDimsYX(self):
+        xmin, xmax, ymin, ymax = self.boundsMap
+        return (ymax-ymin),(xmax-xmin)
+
+    def buffer_imXY(self, buffImX=0, buffImY=0):
+        # type: (float, float)
+        """Buffer the box in X- and/or Y-direction.
+        :param buffImX:     <float> buffer value in x-direction as IMAGE UNITS (pixels)
+        :param buffImY:     <float> buffer value in y-direction as IMAGE UNITS (pixels)
+        """
+        xmin,xmax,ymin,ymax = self.boundsIm
+        xmin,xmax,ymin,ymax = xmin-buffImX, xmax+buffImX, ymin-buffImY, ymax+buffImY
+        self.imPoly         = box(xmin, ymin, xmax, ymax)
+
+    def is_larger_DimXY(self,boundsIm2test):
+        # type: (tuple) -> bool,bool
+        """Checks if the boxObj is larger than a given set of bounding image coordinates (in X- and/or Y-direction).
+        :param boundsIm2test:   <tuple> (xmin,xmax,ymin,ymax) as image coordinates
+        """
+        b2t_xmin, b2t_xmax, b2t_ymin, b2t_ymax = boundsIm2test
+        xmin, xmax, ymin, ymax                 = self.boundsIm
+        x_is_larger = xmin<b2t_xmin or xmax>b2t_xmax
+        y_is_larger = ymin<b2t_ymin or ymax>b2t_ymax
+        return x_is_larger,y_is_larger
+
+
+def get_winPoly(wp_imYX, ws, gt, match_grid=0):
+    # type: (tuple, tuple, list, bool) -> shapely.Polygon, tuple, tuple
+    """Creates a shapely polygon from a given set of image cordinates, window size and geotransform.
+    :param wp_imYX:
+    :param ws:          <tuple> X/Y window size
+    :param gt:
+    :param match_grid:  <bool>
+    """
+    ws = (ws, ws) if isinstance(ws, int) else ws
+    xmin, xmax, ymin, ymax = (
+    wp_imYX[1] - ws[0] / 2, wp_imYX[1] + ws[0] / 2, wp_imYX[0] + ws[1] / 2, wp_imYX[0] - ws[1] / 2)
+    if match_grid:
+        xmin, xmax, ymin, ymax = [int(i) for i in [xmin, xmax, ymin, ymax]]
+    box_YX = (ymax, xmin), (ymax, xmax), (ymin, xmax), (ymin, xmin)  # UL,UR,LR,LL
+    UL, UR, LR, LL = [imYX2mapYX(imYX, gt) for imYX in box_YX]
+    box_mapYX = UL, UR, LR, LL
+    return Polygon([(UL[1], UL[0]), (UR[1], UR[0]), (LR[1], LR[0]), (LL[1], LR[0])]), box_mapYX, box_YX
\ No newline at end of file

ptds/map/vector/topology.py → ptds/geo/vector/topology.py

@@ -2,14 +2,19 @@
 __author__ = "Daniel Scheffler"
 
 
-from shapely.geometry import shape
+import math
 
+from shapely.geometry import shape, Polygon, box
+from ..coord_trafo    import mapYX2imYX
+from ..coord_grid     import find_nearest_grid_coord
 
-def get_overlap_polygon(poly1, poly2):
+
+def get_overlap_polygon(poly1, poly2, v=True):
     """ Returns a dict with the overlap of two shapely.Polygon() objects, the overlap percentage and the overlap area.
 
     :param poly1:   first shapely.Polygon() object
     :param poly2:   second shapely.Polygon() object
+    :param v:       verbose mode
     :return:        overlap polygon as shapely.Polygon() object
     :return:        overlap percentage as float value [%]
     :return:        area of overlap polygon
@@ -17,21 +22,52 @@ def get_overlap_polygon(poly1, poly2):
     overlap_poly = poly1.intersection(poly2)
     if not overlap_poly.is_empty:
         overlap_percentage = 100 * shape(overlap_poly).area / shape(poly2).area
+        if v: print('%.2f percent of the image to be shifted is covered by the reference image.' % overlap_percentage)
         return {'overlap poly':overlap_poly, 'overlap percentage':overlap_percentage, 'overlap area':overlap_poly.area}
     else:
         return {'overlap poly':None, 'overlap percentage':0, 'overlap area':0}
 
 
+def get_footprint_polygon(CornerLonLat):
+    """ Converts a list of coordinates into a shapely polygon object.
+    :param CornerLonLat:    a list of coordinate tuples like [[lon,lat], [lon. lat], ..]
+                            in clockwise or counter clockwise order
+    :return:                a shapyly.Polygon() object
+    """
+    outpoly = Polygon(CornerLonLat)
+    assert outpoly. is_valid, 'The given coordinates result in an invalid polygon. Check coordinate order.'
+    return outpoly
+
+
+def get_smallest_boxImYX_that_contains_boxMapYX(box_mapYX, gt_im):
+    xmin,ymin,xmax,ymax = Polygon([(i[1],i[0]) for i in box_mapYX]).bounds # map-bounds box_mapYX
+    (ymin,xmin),(ymax,xmax) = mapYX2imYX([ymin,xmin],gt_im), mapYX2imYX([ymax,xmax],gt_im) # image coord bounds
+    xmin,ymin,xmax,ymax = int(xmin),math.ceil(ymin), math.ceil(xmax), int(ymax) # round min coords off and max coords on
+    return (ymax,xmin),(ymax,xmax),(ymin,xmax),(ymin,xmin) # UL_YX,UR_YX,LR_YX,LL_YX
+
+
+def get_largest_onGridPoly_within_poly(outerPoly,gt,rows,cols):
+    oP_xmin,oP_ymin,oP_xmax,oP_ymax = outerPoly.bounds
+    xmin,ymax = find_nearest_grid_coord((oP_xmin,oP_ymax),gt,rows,cols,direction='SE')
+    xmax,ymin = find_nearest_grid_coord((oP_xmax,oP_ymin),gt,rows,cols,direction='NW')
+    return box(xmin, ymin, xmax, ymax)
+
+
+def get_smallest_shapelyImPolyOnGrid_that_contains_shapelyImPoly(shapelyPoly):
+    """Returns the smallest box that matches the coordinate grid of the given geotransform.
+    The returned shapely polygon contains image coordinates."""
+    xmin,ymin,xmax,ymax = shapelyPoly.bounds # image_coords-bounds
+
+    return box(int(xmin),int(ymin), math.ceil(xmax), math.ceil(ymax)) # round min coords off and max coords on
+
 def find_line_intersection_point(line1, line2):
     xdiff = (line1[0][0] - line1[1][0], line2[0][0] - line2[1][0])
     ydiff = (line1[0][1] - line1[1][1], line2[0][1] - line2[1][1])
 
-    def det(a, b):
-        return a[0] * b[1] - a[1] * b[0]
-
+    det = lambda a, b: a[0] * b[1] - a[1] * b[0]
     div = det(xdiff, ydiff)
     if div == 0:
-        raise Exception('lines do not intersect')
+        return None, None
     d = (det(*line1), det(*line2))
     x = det(d, xdiff) / div
     y = det(d, ydiff) / div

ptds/io/raster/GeoArray.py

@@ -18,13 +18,12 @@ except ImportError:
     import gdalnumeric
 
 
-
-from ...map        import coord_calc  as CC
-from ...map        import coord_grid  as CG
-from ...map        import coord_trafo as CT
-from ...map        import projection  as PRJ
-from ...map.vector import topology    as TOPO
-from ...map.raster.reproject import warp_ndarray
+from ...geo.coord_calc       import get_corner_coordinates, calc_FullDataset_corner_positions
+from ...geo.coord_grid       import snap_bounds_to_pixGrid
+from ...geo.coord_trafo      import mapXY2imXY, imXY2mapXY, transform_any_prj
+from ...geo.projection       import prj_equal
+from ...geo.vector.topology  import get_overlap_polygon
+from ...geo.raster.reproject import warp_ndarray
 
 
 
@@ -331,7 +330,7 @@ class GeoArray(object):
         B     = shape[2] if len(shape) > 2 else 1
 
         #if not bandslist or (bandslist and sorted(bandslist)==list(range(B))): # FIXME activate if bug is fixed
-        if PRJ.prj_equal(arr_prj,mapBounds_prj):   ############# WORKAROUND FOR corrupt return values in case of 3D array input that has to be warped
+        if prj_equal(arr_prj,mapBounds_prj):   ############# WORKAROUND FOR corrupt return values in case of 3D array input that has to be warped
             """process whole array at once"""
             sub_arr, sub_gt, sub_prj = get_array_at_mapPos(self, arr_gt, arr_prj, mapBounds, mapBounds_prj,
                                                            fillVal=fillVal)
@@ -394,11 +393,11 @@ def _clip_array_at_mapPos(arr, mapBounds, arr_gt, band2clip=None, fillVal=0):
     rows, cols             = arr.shape[:2]
     bands                  = arr.shape[2] if len(arr.shape) == 3 else 1
     arr_dtype              = arr.dtype
-    ULxy, LLxy, LRxy, URxy = CC.get_corner_coordinates(gt