diff --git a/easyidp/geotiff.py b/easyidp/geotiff.py index ff13fc9..7c8f2cb 100644 --- a/easyidp/geotiff.py +++ b/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,9 +700,20 @@ 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): @@ -702,21 +721,29 @@ def pixel2geo(points_hv, header): # 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 diff --git a/easyidp/roi.py b/easyidp/roi.py index 5854a4a..f9af38e 100644 --- a/easyidp/roi.py +++ b/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() diff --git a/tests/test_geotiff.py b/tests/test_geotiff.py index cc23a06..0a878e0 100644 --- a/tests/test_geotiff.py +++ b/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