Merge branch 'revision/revise_coord_trafo' into 'master'

Revision/revise coord trafo

See merge request !24

HISTORY.rst

@@ -2,6 +2,18 @@
 History
 =======
 
+0.16.2 (2020-11-18)
+-------------------
+
+* Fixed issue of remaining coverage artifacts after running 'make clean-test.
+* Revised coord_trafo.py. This fixes an issue that caused pixelToLatLon() to return Lon/Lat instead of Lat/Lon.
+* Fixed mapXY2imXY() and imXY2mapXY().
+* Added Test_mapXY2imXY(), Test_imXY2mapXY(), Test_pixelToLatLon(), Test_latLonToPixel().
+* Removed GDAL dataset input parameters from some functions.
+* Revised code style and some docstrings. Added some typehints.
+* Bugfix for always raising a RuntimeWarning in fill_holes_within_poly().
+
+
 0.16.1 (2020-11-03)
 -------------------
 

Makefile

@@ -43,10 +43,10 @@ clean-pyc: ## remove Python file artifacts
 	find . -name '__pycache__' -exec rm -fr {} +
 
 clean-test: ## remove test and coverage artifacts
-	## don't include 'coverage erase' lib here because clean-test is also executed during package setup and coverage is
-	## only a test requirement
+	## don't call coverage erase here because make install calls make clean which calls make clean-test
+	## -> since make install should run without the test requirements we can't use coverage erase here
 	rm -fr .tox/
-	rm -f .coverage
+	rm -fr .coverage.*
 	rm -fr htmlcov/
 	rm -fr nosetests.html
 	rm -fr nosetests.xml

py_tools_ds/geo/coord_calc.py

@@ -31,18 +31,15 @@ from shapely.geometry import Polygon
 __author__ = "Daniel Scheffler"
 
 
-def get_corner_coordinates(gdal_ds=None, gt=None, cols=None, rows=None):
+def get_corner_coordinates(gt, cols, rows):
     """Returns (ULxy, LLxy, LRxy, URxy) in the same coordinate units like the given geotransform."""
-    assert gdal_ds or (gt and cols and rows), \
-        "GEOP.get_corner_coordinates: Missing argument! Please provide either 'gdal_ds' or 'gt', 'cols' AND 'rows'."
-    gdal_ds_GT = gdal_ds.GetGeoTransform() if gdal_ds else gt
     ext = []
-    xarr = [0, gdal_ds.RasterXSize if gdal_ds else cols]
-    yarr = [0, gdal_ds.RasterYSize if gdal_ds else rows]
+    xarr = [0, cols]
+    yarr = [0, rows]
     for px in xarr:
         for py in yarr:
-            x = gdal_ds_GT[0] + (px * gdal_ds_GT[1]) + (py * gdal_ds_GT[2])
-            y = gdal_ds_GT[3] + (px * gdal_ds_GT[4]) + (py * gdal_ds_GT[5])
+            x = gt[0] + (px * gt[1]) + (py * gt[2])
+            y = gt[3] + (px * gt[4]) + (py * gt[5])
             ext.append([x, y])
         yarr.reverse()
     return ext

py_tools_ds/geo/coord_trafo.py

@@ -24,38 +24,47 @@
 import pyproj
 import numpy as np
 from shapely.ops import transform
-from typing import Union  # noqa F401  # flake8 issue
+from typing import Union, Sequence, List  # noqa F401  # flake8 issue
 from osgeo import gdal, osr
 
 __author__ = "Daniel Scheffler"
 
 
 def transform_utm_to_wgs84(easting, northing, zone, south=False):
+    # type: (float, float, int, bool) -> (float, float)
     """Transform an UTM coordinate to lon, lat, altitude."""
-    UTM = pyproj.Proj(proj='utm', zone=abs(zone), ellps='WGS84', south=(zone < 0 or south))
+    UTM = pyproj.Proj(proj='utm',
+                      zone=abs(zone),
+                      ellps='WGS84',
+                      south=(zone < 0 or south))
     return UTM(easting, northing, inverse=True)
 
 
 def get_utm_zone(longitude):
+    # type: (float) -> int
     """Get the UTM zone corresponding to the given longitude."""
     return int(1 + (longitude + 180.0) / 6.0)
 
 
 def transform_wgs84_to_utm(lon, lat, zone=None):
-    """ returns easting, northing, altitude. If zone is not set it is derived automatically.
-    :param lon:
-    :param lat:
-    :param zone:
+    # type: (float, float, int) -> (float, float)
+    """Return easting, northing, altitude for the given Lon/Lat coordinate.
+    :param lon:     longitude
+    :param lat:     latitude
+    :param zone:    UTM zone (if not set it is derived automatically)
     """
     if not (-180. <= lon <= 180. and -90. <= lat <= 90.):
         raise ValueError((lon, lat), 'Input coordinates for transform_wgs84_to_utm() out of range.')
 
-    UTM = pyproj.Proj(proj='utm', zone=get_utm_zone(lon) if zone is None else zone, ellps='WGS84')
+    UTM = pyproj.Proj(proj='utm',
+                      zone=get_utm_zone(lon) if zone is None else zone,
+                      ellps='WGS84')
     return UTM(lon, lat)
 
 
 def transform_any_prj(prj_src, prj_tgt, x, y):
-    """Transforms X/Y data from any source projection to any target projection.
+    # type: (Union[str, int], Union[str, int], float, float) -> (float, float)
+    """Transform X/Y data from any source projection to any target projection.
 
     :param prj_src:     GDAL projection as WKT string or EPSG code ('epsg:1234' or <EPSG_int>)
     :param prj_tgt:     GDAL projection as WKT string or EPSG code ('epsg:1234' or <EPSG_int>)
@@ -70,7 +79,9 @@ def transform_any_prj(prj_src, prj_tgt, x, y):
 
 
 def transform_coordArray(prj_src, prj_tgt, Xarr, Yarr, Zarr=None):
-    """Transforms geolocation arrays from one projection into another.
+    # type: (str, str, np.ndarray, np.ndarray, np.ndarray) -> Sequence[np.ndarray]
+    """Transform a geolocation array from one projection into another.
+
     HINT: This function is faster than transform_any_prj but works only for geolocation arrays.
 
     :param prj_src:     WKT string
@@ -113,166 +124,136 @@ def transform_coordArray(prj_src, prj_tgt, Xarr, Yarr, Zarr=None):
 
 
 def mapXY2imXY(mapXY, gt):
-    # type: (tuple, Union[list, tuple]) -> Union[tuple, np.ndarray]
-    """Translates given geo coordinates into pixel locations according to the given image geotransform.
+    # type: (Union[tuple, np.ndarray], Sequence) -> Union[tuple, np.ndarray]
+    """Translate the given geo coordinates into pixel locations according to the given image geotransform.
 
     :param mapXY:   <tuple, np.ndarray> The geo coordinates to be translated in the form (x,y) or as np.ndarray [Nx1].
     :param gt:      <list> GDAL geotransform
     :returns:       <tuple, np.ndarray> image coordinate tuple X/Y (column, row) or np.ndarray [Nx2]
     """
     if isinstance(mapXY, np.ndarray):
-        if mapXY.ndim == 1:
+        ndim = mapXY.ndim
+        if ndim == 1:
             mapXY = mapXY.reshape(1, 2)
         assert mapXY.shape[1] == 2, 'An array in shape [Nx2] is needed. Got shape %s.' % mapXY.shape
-        imXY = np.empty_like(mapXY)
+        imXY = np.empty_like(mapXY, dtype=np.float)
         imXY[:, 0] = (mapXY[:, 0] - gt[0]) / gt[1]
         imXY[:, 1] = (mapXY[:, 1] - gt[3]) / gt[5]
-        return imXY
+        return imXY if ndim > 1 else imXY.flatten()
     else:
-        return (mapXY[0] - gt[0]) / gt[1], (mapXY[1] - gt[3]) / gt[5]
+        return (mapXY[0] - gt[0]) / gt[1],\
+               (mapXY[1] - gt[3]) / gt[5]
 
 
 def imXY2mapXY(imXY, gt):
-    # type: (tuple, Union[list, tuple]) -> Union[tuple, np.ndarray]
-    """Translates given pixel locations into geo coordinates according to the given image geotransform.
+    # type: (Union[tuple, np.ndarray], Sequence) -> Union[tuple, np.ndarray]
+    """Translate the given pixel X/Y locations into geo coordinates according to the given image geotransform.
 
     :param imXY:    <tuple, np.ndarray> The image coordinates to be translated in the form (x,y) or as np.ndarray [Nx1].
     :param gt:      <list> GDAL geotransform
     :returns:       <tuple, np.ndarray> geo coordinate tuple X/Y (mapX, mapY) or np.ndarray [Nx2]
     """
     if isinstance(imXY, np.ndarray):
+        ndim = imXY.ndim
         if imXY.ndim == 1:
             imXY = imXY.reshape(1, 2)
         assert imXY.shape[1] == 2, 'An array in shape [Nx2] is needed. Got shape %s.' % imXY.shape
-        imXY = np.empty_like(imXY)
-        imXY[:, 0] = gt[0] + imXY[:, 0] * abs(gt[1])
-        imXY[:, 1] = gt[3] - imXY[:, 1] * abs(gt[5])
-        return imXY
+        mapXY = np.empty_like(imXY, dtype=np.float)
+        mapXY[:, 0] = gt[0] + imXY[:, 0] * abs(gt[1])
+        mapXY[:, 1] = gt[3] - imXY[:, 1] * abs(gt[5])
+        return mapXY if ndim > 1 else mapXY.flatten()
     else:
-        return (gt[0] + imXY[0] * abs(gt[1])), (gt[3] - imXY[1] * abs(gt[5]))
+        return (gt[0] + imXY[0] * abs(gt[1])),\
+               (gt[3] - imXY[1] * abs(gt[5]))
 
 
 def mapYX2imYX(mapYX, gt):
-    return (mapYX[0] - gt[3]) / gt[5], (mapYX[1] - gt[0]) / gt[1]
+    return (mapYX[0] - gt[3]) / gt[5],\
+           (mapYX[1] - gt[0]) / gt[1]
 
 
 def imYX2mapYX(imYX, gt):
-    return gt[3] - (imYX[0] * abs(gt[5])), gt[0] + (imYX[1] * gt[1])
+    return 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()
+def pixelToMapYX(pixelCoords, geotransform, projection):
+    # MapXYPairs = imXY2mapXY(np.array(pixelPairs), geotransform)
+
+    # Create a spatial reference object
     srs = osr.SpatialReference()
-    srs.ImportFromWkt(proj)
+    srs.ImportFromWkt(projection)
+
     # 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])
+    mapXYPairs = [ct.TransformPoint(pixX * geotransform[1] + geotransform[0],
+                                    pixY * geotransform[5] + geotransform[3])[:2]
+                  for pixX, pixY in pixelCoords]
 
-    return mapYmapXPairs
+    return [[mapY, mapX] for mapX, mapY in mapXYPairs]
 
 
-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
-    image resolution for pixel size to be insignificant
+def pixelToLatLon(pixelPairs, geotransform, projection):
+    # type: (List[Sequence[float, float]], Sequence, str) -> List[Sequence[float, float]]
+    """Translate the given pixel X/Y locations into latitude/longitude locations.
 
     :param pixelPairs:      Image coordinate pairs to be translated in the form [[x1,y1],[x2,y2]]
-    :param path_im:         The file location of the input image
     :param geotransform:    GDAL GeoTransform
-    :param projection:      GDAL Projection
+    :param projection:      WKT string
     :returns:               The lat/lon translation of the pixel pairings in the form [[lat1,lon1],[lat2,lon2]]
     """
-    assert path_im or geotransform and projection, \
-        "GEOP.pixelToLatLon: 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()
+    MapXYPairs = imXY2mapXY(np.array(pixelPairs), geotransform)
+    lons, lats = transform_any_prj(projection, 4326, MapXYPairs[:, 0], MapXYPairs[:, 1])
+
+    # validate output lat/lon values
+    if not (-180 <= np.min(lons) <= 180) or not (-180 <= np.max(lons) <= 180):
+        raise RuntimeError('Output longitude values out of bounds.')
+    if not (-90 <= np.min(lats) <= 90) or not (-90 <= np.max(lats) <= 90):
+        raise RuntimeError('Output latitudes values out of bounds.')
+
+    LatLonPairs = np.vstack([lats, lons]).T.tolist()
+
+    return LatLonPairs
 
-    # Create a spatial reference object for the dataset
-    srs = osr.SpatialReference()
-    srs.ImportFromWkt(proj)
-    # Set up the coordinate transformation object
-    srsLatLong = srs.CloneGeogCS()
-    ct = osr.CoordinateTransformation(srs, srsLatLong)
-    # Go through all the point pairs and translate them to pixel pairings
-    latLonPairs = []
-    for point in pixelPairs:
-        # Translate the pixel pairs into untranslated points
-        ulon = point[0] * gt[1] + gt[0]
-        ulat = point[1] * gt[5] + gt[3]
-        # Transform the points to the space
-        (lon, lat, holder) = ct.TransformPoint(ulon, ulat)
-        # Add the point to our return array
-        latLonPairs.append([lat, lon])
-
-    return latLonPairs
-
-
-def latLonToPixel(latLonPairs, path_im=None, geotransform=None, projection=None):
-    """The following method translates given latitude/longitude pairs into pixel locations on a given GEOTIF
-    This method does not take into account pixel size and assumes a high enough
-    image resolution for pixel size to be insignificant
+
+def latLonToPixel(latLonPairs, geotransform, projection):
+    # type: (List[Sequence[float, float]], Sequence, str) -> List[Sequence[float, float]]
+    """Translate the given latitude/longitude pairs into pixel X/Y locations.
 
     :param latLonPairs:     The decimal lat/lon pairings to be translated in the form [[lat1,lon1],[lat2,lon2]]
-    :param path_im:         The file location of the input image
     :param geotransform:    GDAL GeoTransform
-    :param projection:      GDAL Projection
+    :param projection:      WKT string
     :return:            The pixel translation of the lat/lon pairings in the form [[x1,y1],[x2,y2]]
     """
-    assert path_im or geotransform and projection, \
-        "GEOP.latLonToPixel: 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()
-    # Load the image dataset
-    # Create a spatial reference object for the dataset
-    srs = osr.SpatialReference()
-    srs.ImportFromWkt(proj)
-    # Set up the coordinate transformation object
-    srsLatLong = srs.CloneGeogCS()
-    ct = osr.CoordinateTransformation(srsLatLong, srs)
-    # Go through all the point pairs and translate them to latitude/longitude pairings
-    pixelPairs = []
-    for point in latLonPairs:
-        # Change the point locations into the GeoTransform space
-        (point[1], point[0], holder) = ct.TransformPoint(point[1], point[0])
-        # Translate the x and y coordinates into pixel values
-        x = (point[1] - gt[0]) / gt[1]
-        y = (point[0] - gt[3]) / gt[5]
-        # Add the point to our return array
-        pixelPairs.append([int(x), int(y)])
+    latLonArr = np.array(latLonPairs)
+
+    # validate input lat/lon values
+    if not (-180 <= np.min(latLonArr[:, 1]) <= 180) or not (-180 <= np.max(latLonArr[:, 1]) <= 180):
+        raise ValueError('Longitude value out of bounds.')
+    if not (-90 <= np.min(latLonArr[:, 0]) <= 90) or not (-90 <= np.max(latLonArr[:, 0]) <= 90):
+        raise ValueError('Latitude value out of bounds.')
+
+    mapXs, mapYs = transform_any_prj(4326, projection, latLonArr[:, 1], latLonArr[:, 0])
+    imXYarr = mapXY2imXY(np.vstack([mapXs, mapYs]).T, geotransform)
+
+    pixelPairs = imXYarr.tolist()
+
     return pixelPairs
 
 
 def lonlat_to_pixel(lon, lat, inverse_geo_transform):
-    """Translates the given lon, lat to the grid pixel coordinates in data array (zero start)"""
+    """Translate the given lon, lat to the grid pixel coordinates in data array (zero start)."""
     # transform to utm
-    utm_x, utm_y, utm_z = transform_wgs84_to_utm(lon, lat)
-    # apply inverse geo tranform
+    utm_x, utm_y = transform_wgs84_to_utm(lon, lat)
+
+    # apply inverse geo transform
     pixel_x, pixel_y = gdal.ApplyGeoTransform(inverse_geo_transform, utm_x, utm_y)
     pixel_x = int(pixel_x) - 1  # adjust to 0 start for array
     pixel_y = int(pixel_y) - 1  # adjust to 0 start for array
+
     return pixel_x, abs(pixel_y)  # y pixel is likly a negative value given geo_transform
 
 
@@ -291,7 +272,7 @@ def transform_GCPlist(gcpList, prj_src, prj_tgt):
 
 
 def reproject_shapelyGeometry(shapelyGeometry, prj_src, prj_tgt):
-    """Reprojects any shapely geometry from one projection to another.
+    """Reproject any shapely geometry from one projection to another.
 
     :param shapelyGeometry: any shapely geometry instance
     :param prj_src:         GDAL projection as WKT string or EPSG code ('epsg:1234' or <EPSG_int>)

py_tools_ds/geo/projection.py

@@ -156,11 +156,11 @@ def WKT2EPSG(wkt):
     return CRS.from_wkt(wkt.replace('\n', '').replace('\r', '').replace(' ', '')).to_epsg()
 
 
-def get_UTMzone(ds=None, prj=None):
-    assert ds or prj, "Specify at least one of the arguments 'ds' or 'prj'"
+def get_UTMzone(prj):
+    # type: (str) -> Union[int, None]
     if isProjectedOrGeographic(prj) == 'projected':
         srs = osr.SpatialReference()
-        srs.ImportFromWkt(ds.GetProjection() if ds else prj)
+        srs.ImportFromWkt(prj)
         return srs.GetUTMZone()
     else:
         return None

py_tools_ds/geo/raster/reproject.py

@@ -132,7 +132,6 @@ def warp_ndarray_rasterio(ndarray, in_gt, in_prj, out_prj, out_gt=None, outRowsC
 
             else:  # outRowsCols is None and outUL is None: use in_gt
                 left, bottom, right, top = gt2bounds(in_gt, rows, cols)
-                # ,im_xmax,im_ymin,im_ymax = corner_coord_to_minmax(get_corner_coordinates(self.ds_im2shift))
 
         elif out_extent:
             left, bottom, right, top = out_extent
@@ -158,8 +157,8 @@ def warp_ndarray_rasterio(ndarray, in_gt, in_prj, out_prj, out_gt=None, outRowsC
                 raise NotImplementedError('Different projections are currently not supported.')
 
             else:  # ('projected','geographic')
-                px_size_LatLon = np.array(pixelToLatLon([1, 1], geotransform=in_gt, projection=in_prj)) - \
-                                 np.array(pixelToLatLon([0, 0], geotransform=in_gt, projection=in_prj))
+                px_size_LatLon = np.array(pixelToLatLon([[1, 1]], geotransform=in_gt, projection=in_prj)) - \
+                                 np.array(pixelToLatLon([[0, 0]], geotransform=in_gt, projection=in_prj))
                 out_res = tuple(reversed(abs(px_size_LatLon)))
                 print('OUT_RES NOCHMAL CHECKEN: ', out_res)
 

py_tools_ds/geo/vector/conversion.py

@@ -56,8 +56,18 @@ def shapelyImPoly_to_shapelyMapPoly(shapelyBox, gt):
 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))
+
+    if coord_type == 'image':
+        UL_YX = ymin, xmin
+        UR_YX = ymin, xmax
+        LR_YX = ymax, xmax
+        LL_YX = ymax, xmin
+    else:
+        UL_YX = ymax, xmin
+        UR_YX = ymax, xmax
+        LR_YX = ymin, xmax
+        LL_YX = ymin, xmin
+
     return UL_YX, UR_YX, LR_YX, LL_YX
 
 

py_tools_ds/geo/vector/topology.py

@@ -198,7 +198,7 @@ def fill_holes_within_poly(poly):
             GDF_row.geometry.equals(largest_poly_filled),
             axis=1)
 
-        if False in gdf.within_or_equal:
+        if False in gdf.within_or_equal.values:
             n_disjunct_polys = int(np.sum(~gdf.within_or_equal))
             warnings.warn(RuntimeWarning('The given MultiPolygon contains %d disjunct polygone(s) outside of the '
                                          'largest polygone. fill_holes_within_poly() will only return the largest '

tests/test_geo/test_coord_trafo.py

@@ -34,13 +34,53 @@ from unittest import TestCase, main
 from shapely.geometry import Polygon
 import numpy as np
 
-from py_tools_ds.geo.coord_trafo import reproject_shapelyGeometry, transform_any_prj
-
-
-poly_local = Polygon([(5708.2, -3006), (5708, -3262), (5452, -3262), (5452, -3006), (5708, -3006)])
-
-ll_x, ll_y = (10.701359, 47.54536)
-utm_x, utm_y = (628023.1302739962, 5267173.503313034)
+from py_tools_ds.geo.projection import EPSG2WKT
+from py_tools_ds.geo.coord_trafo import (
+    imXY2mapXY, mapXY2imXY,
+    pixelToLatLon, latLonToPixel,
+    reproject_shapelyGeometry,
+    transform_any_prj,
+    get_utm_zone
+)
+
+poly_local = Polygon([(5708, -3006),
+                      (5708, -3262),
+                      (5452, -3262),
+                      (5452, -3006),
+                      (5708, -3006)])
+
+ll_x, ll_y = (10.701359,
+              47.54536)
+utm_x, utm_y = (628023.1302739962,
+                5267173.503313034)
+
+geotransform_utm = (331185.0, 30.0, -0.0, 5840115.0, -0.0, -30.0)
+wkt_utm = \
+    """
+    PROJCS["WGS 84 / UTM zone 33N",
+           GEOGCS["WGS 84",
+                  DATUM["WGS_1984",
+                        SPHEROID["WGS 84", 6378137, 298.257223563,
+                                 AUTHORITY["EPSG", "7030"]],
+                        AUTHORITY["EPSG", "6326"]],
+                  PRIMEM["Greenwich", 0,
+                         AUTHORITY["EPSG", "8901"]],
+                  UNIT["degree", 0.0174532925199433,
+                       AUTHORITY["EPSG", "9122"]],
+                  AUTHORITY["EPSG", "4326"]],
+           PROJECTION["Transverse_Mercator"],
+           PARAMETER["latitude_of_origin", 0],
+           PARAMETER["central_meridian", 15],
+           PARAMETER["scale_factor", 0.9996],
+           PARAMETER["false_easting", 500000],
+           PARAMETER["false_northing", 0],
+           UNIT["metre", 1,
+                AUTHORITY["EPSG", "9001"]],
+           AXIS["Easting", EAST],
+           AXIS["Northing", NORTH],
+           AUTHORITY["EPSG", "32633"]]
+    """
+wkt_utm = ' '.join(wkt_utm.split())
 
 
 class Test_reproject_shapelyGeometry(TestCase):
@@ -60,5 +100,150 @@ class Test_transform_any_prj(TestCase):
                                           np.array([utm_x, utm_y]))))
 
 
+def get_utm_gt_prj_from_lonlat(lon, lat):
+    utmepsg = int(f"32{6 if lat > 0 else 7}{get_utm_zone(lon)}")
+    utmx, utmy = transform_any_prj(4326, utmepsg, lon, lat)
+    gt = (utmx, 30, 0, utmy, 0, -30)
+    prj = EPSG2WKT(utmepsg)
+
+    return gt, prj
+
+
+class Test_mapXY2imXY(TestCase):
+
+    def test_singlecoord_tuple_at_origin(self):
+        out = mapXY2imXY((331185.0, 5840115.0), geotransform_utm)
+        self.assertIsInstance(out, tuple)
+        self.assertEqual(out, (0, 0))
+
+    def test_singlecoord_tuple_at_1_1(self):
+        out = mapXY2imXY((331215.0, 5840085.0), geotransform_utm)
+        self.assertIsInstance(out, tuple)
+        self.assertEqual(out, (1, 1))
+
+    def test_singlecoord_nparray_at_1_1(self):
+        out = mapXY2imXY(np.array([331215.0, 5840085.0]), geotransform_utm)
+        self.assertIsInstance(out, np.ndarray)
+        self.assertTrue(np.array_equal(out, np.array([1, 1])))
+
+    def test_multicoords_nparray(self):
+        out = mapXY2imXY(np.array([[331185.0, 5840115.0],
+                                   [331215.0, 5840085.0]]), geotransform_utm)
+        self.assertIsInstance(out, np.ndarray)
+        self.assertTrue(np.array_equal(out, np.array([[0, 0], [1, 1]])))
+
+    def test_tuple_ndarray_outputConsistency(self):
+        tupleout = mapXY2imXY((331215.0, 5840085.0), geotransform_utm)
+        npout = mapXY2imXY(np.array([331215.0, 5840085.0]), geotransform_utm)
+        self.assertTrue(np.array_equal(np.array(tupleout),
+                                       npout))
+
+
+class Test_imXY2mapXY(TestCase):
+
+    def test_singlecoord_tuple_at_origin(self):
+        out = imXY2mapXY((0, 0), geotransform_utm)
+        self.assertIsInstance(out, tuple)
+        self.assertEqual(out, (331185.0, 5840115.0))
+
+    def test_singlecoord_tuple_at_1_1(self):
+        out = imXY2mapXY((1, 1), geotransform_utm)
+        self.assertIsInstance(out, tuple)
+        self.assertEqual(out, (331215.0, 5840085.0))
+
+    def test_singlecoord_nparray_at_1_1(self):
+        out = imXY2mapXY(np.array([1, 1]), geotransform_utm)
+        self.assertIsInstance(out, np.ndarray)
+        self.assertTrue(np.array_equal(out, np.array([331215.0, 5840085.0])))
+
+    def test_multicoords_nparray(self):
+        out = imXY2mapXY(np.array([[0, 0], [1, 1]]), geotransform_utm)
+        self.assertIsInstance(out, np.ndarray)
+        self.assertTrue(np.array_equal(out, np.array([[331185.0, 5840115.0],
+                                                      [331215.0, 5840085.0]])))
+
+    def test_tuple_ndarray_outputConsistency(self):
+        tupleout = imXY2mapXY((1, 1), geotransform_utm)
+        npout = imXY2mapXY(np.array([1, 1]), geotransform_utm)
+        self.assertTrue(np.array_equal(np.array(tupleout),
+                                       npout))
+
+
+class Test_pixelToLatLon(TestCase):
+
+    def test_standard_case(self):
+        latlonPairs = pixelToLatLon([(0, 0)], geotransform=geotransform_utm, projection=wkt_utm)
+        self.assertIsInstance(latlonPairs, list)
+        self.assertEqual(len(latlonPairs), 1)
+
+    def test_negative_imcoords(self):
+        latlonPairs = pixelToLatLon([(-10, -10)], geotransform=geotransform_utm, projection=wkt_utm)
+        self.assertIsInstance(latlonPairs, list)
+        self.assertEqual(len(latlonPairs), 1)
+
+    def test_negative_lon(self):
+        for lon, lat in [(-1, 30), (-89, -30), (-90, -30), (-91, 30), (-130, -30), (-79.999999, -30)]:
+            latlonPairs = pixelToLatLon([(0, 0)], *get_utm_gt_prj_from_lonlat(lon=lon, lat=lat))
+            self.assertTrue(np.allclose(np.array(latlonPairs),
+                                        np.array([(lat, lon)])),
+                            msg=f"{latlonPairs} != {[(lat, lon)]}")
+
+    def test_negative_lat(self):
+        for lon, lat in [(-1, -1), (1, -30), (-1, -45), (1, -89.999999)]:
+            latlonPairs = pixelToLatLon([(0, 0)], *get_utm_gt_prj_from_lonlat(lon=lon, lat=lat))
+            self.assertTrue(np.allclose(np.array(latlonPairs),
+                                        np.array([(lat, lon)])),
+                            msg=f"{latlonPairs} != {[(lat, lon)]}")
+
+    def test_position_out_of_bounds(self):
+        with self.assertRaises(RuntimeError):
+            pixelToLatLon([(0, 0)], (0, 1, 0, -91, 0, -1), EPSG2WKT(4326))
+        with self.assertRaises(RuntimeError):
+            pixelToLatLon([(0, 0)], (361, 1, 0, 0, 0, -1), EPSG2WKT(4326))
+
+    def test_multiple_coords(self):
+        latlonPairs = pixelToLatLon([(-10, -10), (10, 10)], geotransform=geotransform_utm, projection=wkt_utm)
+        self.assertIsInstance(latlonPairs, list)