fix <espg::4326> geo2pixel2geo error

easyidp/geotiff.py

@@ -616,7 +616,7 @@ def geo2pixel(points_hv, header, return_index=False):
 
     Notes
     -----
-    Please note: gis coordinate, horizontal is x axis, vertical is y axis, origin at left upper.
+    Please note: gis UTM coordinate, horizontal is x axis, vertical is y axis, origin at left upper.
 
     To crop image ndarray:
 
@@ -652,17 +652,25 @@ def geo2pixel(points_hv, header, return_index=False):
                [16972, 26086]])
 
     """
+    if header['crs'].is_geographic:  # [lat, lon] -> vertical, horizontal -> y, x
+        # but the tie point and scale is [lon, lat] by tifffile.
+        gis_ph = points_hv[:, 1]
+        gis_pv = points_hv[:, 0]
+    else:   # is_projected or is_geocentric seems x,y,z order are correct, please refer metashape.convert_proj3d()
+        gis_ph = points_hv[:, 0]
+        gis_pv = points_hv[:, 1]
+
     gis_xmin = header['tie_point'][0]
     gis_ymax = header['tie_point'][1]
 
-    gis_ph = points_hv[:, 0]
-    gis_pv = points_hv[:, 1]
+    scale_x = header['scale'][0]
+    scale_y = header['scale'][1]
 
     # get float coordinate on pixels
     # - numpy_axis1 = x
-    np_ax_h = (gis_ph - gis_xmin) / header['scale'][0]
+    np_ax_h = (gis_ph - gis_xmin) / scale_x
     # - numpy_axis0 = y
-    np_ax_v = (gis_ymax - gis_pv) / header['scale'][1]
+    np_ax_v = (gis_ymax - gis_pv) / scale_y
 
     # get the pixel index (int)
     if return_index:  
@@ -692,31 +700,50 @@ def pixel2geo(points_hv, header):
     --------
     easyidp.geotiff.get_header
     """
+    if header['crs'].is_geographic:  # [lat, lon] -> vertical, horizontal -> y, x
+        # but the tie point and scale is [lon, lat] by tifffile.
+        gis_ph = points_hv[:, 1]
+        gis_pv = points_hv[:, 0]
+    else:   # is_projected or is_geocentric seems x,y,z order are correct, please refer metashape.convert_proj3d()
+        gis_ph = points_hv[:, 0]
+        gis_pv = points_hv[:, 1]
+
     gis_xmin = header['tie_point'][0]
     gis_ymax = header['tie_point'][1]
 
+    scale_x = header['scale'][0]
+    scale_y = header['scale'][1]
+
     # the px is numpy axis0 (vertical, h)
     # py is numpy axis1 (horizontal, w)
     if np.issubdtype(points_hv.dtype, np.integer):
         # all integer possible means the pixel index 
         #    rather than specific coordinates
         # +0.5 to get the pixel center rather than edge
         # but in QGIS, this will cause 0.5 pixel shift
-        pix_ph = points_hv[:, 0]  # + 0.5
-        pix_pv = points_hv[:, 1]  # + 0.5
+        # pix_ph = points_hv[:, 0]  # + 0.5
+        # pix_pv = points_hv[:, 1]  # + 0.5
+        pass
     elif np.issubdtype(points_hv.dtype, np.floating):
         # all floats possible means it is the pixel coordinates
         #    rather than pixel index
         # no need to +0.5 as image center
-        pix_ph = points_hv[:, 0]
-        pix_pv = points_hv[:, 1]
+        # pix_ph = points_hv[:, 0]
+        # pix_pv = points_hv[:, 1]
+        pass
     else:
         raise TypeError(f"The `points_hv` only accept numpy ndarray integer and float types")
 
-    gis_px = gis_xmin + pix_ph * header['scale'][0]
-    gis_py = gis_ymax - pix_pv * header['scale'][1]
+    pix_ph = points_hv[:, 0]
+    pix_pv = points_hv[:, 1]
 
-    gis_geo = np.vstack([gis_px, gis_py]).T
+    gis_px = gis_xmin + pix_ph * scale_x
+    gis_py = gis_ymax - pix_pv * scale_y
+
+    if header['crs'].is_geographic:  # [lat, lon] -> vertical, horizontal -> y, x
+        gis_geo = np.vstack([gis_py, gis_px]).T
+    else:
+        gis_geo = np.vstack([gis_px, gis_py]).T
 
     return gis_geo
 

easyidp/roi.py

@@ -268,7 +268,7 @@ def get_z_from_dsm(self, dsm, mode="face", kernel="mean", buffer=0, keep_crs=Fal
         else:
             global_z = None
 
-        # convert CRS is necessary
+        # convert CRS if necessary
         if self.crs.name != dsm.header["crs"].name and not keep_crs:
             self.change_crs(dsm.header["crs"])
             poly_dict = self.id_item.copy()

tests/test_geotiff.py

@@ -52,7 +52,7 @@ def test_def_get_imarray():
     assert lotus_part_np.shape == (lh["height"], lh["width"], lh["dim"])
 
 
-def test_def_geo2pixel2geo():
+def test_def_geo2pixel2geo_UTM():
     gis_coord = np.asarray([
         [ 484593.67474654, 3862259.42413431],
         [ 484593.41064743, 3862259.92582402],
@@ -69,7 +69,7 @@ def test_def_geo2pixel2geo():
     header = {'width': 19436, 'height': 31255, 'dim':4, 
               'scale': [0.001, 0.001], 'nodata': None,
               'tie_point': [484576.70205, 3862285.5109300003], 
-              'proj': pyproj.CRS.from_string("WGS 84 / UTM zone 53N")}
+              'crs': pyproj.CRS.from_string("WGS 84 / UTM zone 53N")}
 
     expected_pixel_idx = np.array([
         [16972, 26086],
@@ -107,6 +107,38 @@ def test_def_geo2pixel2geo():
     np.testing.assert_almost_equal(gis_revert_flt, gis_coord, decimal=3)
 
 
+def test_def_geo2pixel2geo_lonlat():
+    # using the source: https://github.com/UTokyo-FieldPhenomics-Lab/EasyIDP/discussions/44
+    gis_latlon_coord = np.array([
+        [ 25.78354364, -80.83957435],
+        [ 25.78354364, -80.83947435],
+        [ 25.78344364, -80.83947435],
+        [ 25.78344364, -80.83957435],
+        [ 25.78354364, -80.83957435]])
+
+    header = {
+        'height': 8748, 'width': 7941, 'dim': 1, 'nodata': -32767.0, 
+        'scale': [3.49222000000852e-07, 3.1617399999982425e-07], 
+        'tie_point': [-80.84039234705898, 25.784493471936425], 
+        'crs': pyproj.CRS.from_epsg(4326)
+    }
+
+    expected_pixel = np.array([
+        [2342.34114395, 3004.14308711],
+        [2628.6919466 , 3004.14308711],
+        [2628.6919466 , 3320.42462829],
+        [2342.34114395, 3320.42462829],
+        [2342.34114395, 3004.14308711]])
+
+    out = idp.geotiff.geo2pixel(gis_latlon_coord, header)
+
+    np.testing.assert_almost_equal(out, expected_pixel)
+
+    back = idp.geotiff.pixel2geo(out, header)
+
+    np.testing.assert_almost_equal(back, gis_latlon_coord, decimal=3)
+
+
 def test_def_tifffile_crop():
     maize_part_np = idp.geotiff.get_imarray(test_data.pix4d.maize_dom)
     # untiled tiff but 2 row as a patch