function updates for dome for Florida

* fix geo coords -> (lon, lat, alt) order
* Add progress bar for `get_z_from_dsm` #30 and roi `crop`
* provide function wrapper for distance filtering #37

docs/index.rst

@@ -47,6 +47,13 @@ EasyIDP (Easy Intermediate Data Processor) is a handy tool for dealing with regi
       * `日本語(翻訳募集) <https://easyidp.readthedocs.io/ja/latest/>`_
 
 
+
+Coordinates Order
+=================
+
+In the EasyIDP, we use the (horizontal, vertical, dim) order as the coords order. When it applies to the GIS coordintes, is (longitude, latitude, altitude)
+
+
 Examples
 ========
 

easyidp/__init__.py

@@ -1,4 +1,4 @@
-__version__ = "2.0.0.dev3"
+__version__ = "2.0.0.dev4"
 
 import os
 import warnings

easyidp/geotiff.py

@@ -3,6 +3,7 @@
 import numpy as np
 import tifffile as tf
 import warnings
+from tqdm import tqdm
 from pathlib import Path
 from pyproj.exceptions import CRSError
 
@@ -149,8 +150,9 @@ def crop(self, roi, is_geo=True, save_folder=None):
         if not isinstance(roi, (dict, idp.ROI)):
             raise TypeError(f"Only <dict> and <easyidp.ROI> are accepted, not {type(roi)}")
 
+        pbar = tqdm(roi.items(), desc=f"Crop roi from geotiff [{os.path.basename(self.file_path)}]")
         out_dict = {}
-        for k, polygon_hv in roi.items():
+        for k, polygon_hv in pbar:
             if save_folder is not None and Path(save_folder).exists():
                 save_path = Path(save_folder) / (k + ".tif")
             else:
@@ -548,21 +550,28 @@ def get_header(tif_path):
         # -> 'ModelTiepoint': [0.0, 0.0, 0.0, 419509.89816000004, 3987344.8286, 0.0]
         header["tie_point"] = page.geotiff_tags["ModelTiepoint"][3:5]
 
-        # pix4d:
+        # pix4d UTM CRS:
         #    page.geotiff_tags
         #    -> 'GTCitationGeoKey': 'WGS 84 / UTM zone 54N'
-        # metashape:
-        #    page.geotiff_tags
-        #    -> 'PCSCitationGeoKey': 'WGS 84 / UTM zone 54N'
-        # metashape: wgs 84
-        #    page.geotiff_tags
-        #    -> 'GeogCitationGeoKey': 'WGS 84'
-        for k, v in page.geotiff_tags.items():
-            if "CitationGeoKey" in k:
-                proj_str = v
-                break
+        if "GTCitationGeoKey" in page.geotiff_tags.keys():
+            proj_str = page.geotiff_tags["GTCitationGeoKey"]
+        # metashape UTM CRS:
+        #     page.geotiff_tags
+        #     -> 'PCSCitationGeoKey': 'WGS 84 / UTM zone 54N'
+        #     -> 'GeogCitationGeoKey': 'WGS 84'
+        elif "PCSCitationGeoKey" in page.geotiff_tags.keys():
+            proj_str = page.geotiff_tags["PCSCitationGeoKey"]
         else:
-            raise KeyError(f"Can not find key '**CitationGeoKey' in Geotiff tages {page.geotiff_tags}")
+            # metashape: wgs 84
+            #    page.geotiff_tags
+            #    -> 'GeogCitationGeoKey': 'WGS 84'
+            for k, v in page.geotiff_tags.items():
+                if "CitationGeoKey" in k:
+                    proj_str = v
+                    print(f"Could not find prefered [GTCitationGeoKey, PCSCItationGeoKey], but find [{k}]={v} instead for Coordinate Reference System (CRS)")
+                    break
+            else:
+                raise KeyError(f"Can not find Coordinate Reference System (CRS) keys '**CitationGeoKey' in Geotiff tages {page.geotiff_tags}")
 
         try:
             crs = pyproj.CRS.from_string(proj_str)
@@ -652,13 +661,8 @@ 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_ph = points_hv[:, 0]
+    gis_pv = points_hv[:, 1]
 
     gis_xmin = header['tie_point'][0]
     gis_ymax = header['tie_point'][1]
@@ -700,13 +704,8 @@ 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_ph = points_hv[:, 0]
+    gis_pv = points_hv[:, 1]
 
     gis_xmin = header['tie_point'][0]
     gis_ymax = header['tie_point'][1]
@@ -740,10 +739,7 @@ def pixel2geo(points_hv, header):
     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
+    gis_geo = np.vstack([gis_px, gis_py]).T
 
     return gis_geo
 

easyidp/metashape.py

@@ -219,7 +219,7 @@ def back2raw_crs(self, points_xyz, distortion_correct=True, save_folder=None, ig
 
             # find out those correct images
             if log:
-                print(f'[Calculator][Judge]{photo.label}w:{photo.w}h:{photo.h}->', end='')
+                print(f'[Calculator][Judge]{photo.label}w:{photo.sensor.width}h:{photo.sensor.height}->', end='')
 
             coords = photo.sensor.in_img_boundary(projected_coord, ignore=ignore, log=log)
             if coords is not None:
@@ -295,6 +295,29 @@ def get_photo_position(self, to_crs=None):
 
         return out
 
+
+    def sort_img_by_distance(self, img_dict_all, roi, distance_thresh=None, num=None):
+        """Advanced wrapper of sorting back2raw img_dict results by distance from photo to roi
+
+        Parameters
+        ----------
+        img_dict_all : dict
+            All output dict of roi.back2raw(...)
+            e.g. img_dict = roi.back2raw(...) -> img_dict
+        roi: idp.ROI
+            Your roi variable
+        num : None or int
+            Keep the closest {x} images
+        distance_thresh : None or float
+            Keep the images closer than this distance to ROI.
+
+        Returns
+        -------
+        dict
+            the same structure as output of roi.back2raw(...)
+        """
+        return idp.reconstruct.sort_img_by_distance(self, img_dict_all, roi, distance_thresh, num)
+
 ###############
 # zip/xml I/O #
 ###############
@@ -1102,9 +1125,9 @@ def convert_proj3d(points_np, crs_origin, crs_target, is_xyz=True):
             x, y, z = ts.transform(*points_np.T)
             out =  np.vstack([x, y, z]).T
         elif crs_target.is_geographic:
-            lat, lon, alt = ts.transform(*points_np.T)
+            lon, lat, alt = ts.transform(*points_np.T)
             # the pyproj output order is reversed
-            out = np.vstack([lon, lat, alt]).T
+            out = np.vstack([lat, lon, alt]).T
         elif crs_target.is_projected:
             lat_m, lon_m, alt_m = ts.transform(*points_np.T)
             out = np.vstack([lat_m, lon_m, alt_m]).T
@@ -1121,9 +1144,10 @@ def convert_proj3d(points_np, crs_origin, crs_target, is_xyz=True):
             out = np.vstack([lon, lat, alt]).T
         elif crs_target.is_projected and crs_target.is_derived:
             lat_m, lon_m, alt_m = ts.transform(lat, lon, alt)
-            out = np.vstack([lat_m, lon_m, alt_m]).T
+            out = np.vstack([lon_m, lat_m, alt_m]).T
         else:
             raise TypeError(f"Given crs is neither `crs.is_geocentric=True` nor `crs.is_geographic` nor `crs.is_projected`")
+
     if is_single:
         return out[0, :]
     else:

easyidp/pix4d.py

@@ -468,6 +468,28 @@ def get_photo_position(self, to_crs=None):
 
         return out
 
+    def sort_img_by_distance(self, img_dict_all, roi, distance_thresh=None, num=None):
+        """Advanced wrapper of sorting back2raw img_dict results by distance from photo to roi
+
+        Parameters
+        ----------
+        img_dict_all : dict
+            All output dict of roi.back2raw(...)
+            e.g. img_dict = roi.back2raw(...) -> img_dict
+        roi: idp.ROI
+            Your roi variable
+        num : None or int
+            Keep the closest {x} images
+        distance_thresh : None or float
+            Keep the images closer than this distance to ROI.
+
+        Returns
+        -------
+        dict
+            the same structure as output of roi.back2raw(...)
+        """
+        return idp.reconstruct.sort_img_by_distance(self, img_dict_all, roi, distance_thresh, num)
+
 ####################
 # code for file IO #
 ####################

easyidp/pointcloud.py

@@ -3,6 +3,7 @@
 from datetime import datetime
 from tabulate import tabulate
 from pathlib import Path
+from tqdm import tqdm
 
 import numpy as np
 import numpy.lib.recfunctions as rfn
@@ -121,7 +122,7 @@ def offset(self, o):
         if self._points is not None:
             self._update_btf_print()
 
-    def set_offset(self, o):
+    def set_offset(self, o, show_warn=True):
         # the point values not change
         # --------------------------------
         # points =  _point + offset
@@ -132,7 +133,8 @@ def set_offset(self, o):
         if self._points is not None:
             self._points = self._points + self._offset - o
             self._offset = o
-            warnings.warn("This will not change the value of point xyz values, if you want to just change offset value, please operate `pcd._offset = offset; pcd._update_btf_print()` directly")
+            if show_warn:
+                warnings.warn("This will not change the value of point xyz values, if you want to just change offset value, please operate `pcd._offset = offset; pcd._update_btf_print()` directly")
             self._update_btf_print()
         else:
             self._offset = o
@@ -280,7 +282,8 @@ def crop(self, roi, save_folder=None):
             raise TypeError(f"Only <dict> and <easyidp.ROI> are accepted, not {type(roi)}")
 
         out_dict = {}
-        for k, polygon_hv in roi.items():
+        pbar = tqdm(roi.items(), desc=f"Crop roi from point cloud [{os.path.basename(self.file_path)}]")
+        for k, polygon_hv in pbar:
             if save_folder is not None and Path(save_folder).exists():
                 save_path = Path(save_folder) / (k + self.file_ext)
             else:
@@ -347,7 +350,7 @@ def crop_point_cloud(self, polygon_xy):
             # create new Point Cloud object
             crop_pcd = PointCloud()
             crop_pcd.points = self.points[pick_idx, :]
-            crop_pcd.set_offset(self.offset)
+            crop_pcd.set_offset(self.offset, show_warn=False)
             if self.has_colors():
                 crop_pcd.colors = self.colors[pick_idx, :]
             if self.has_normals():

easyidp/reconstruct.py

@@ -2,6 +2,7 @@
 import pyproj
 import numpy as np
 from pathlib import Path
+from tqdm import tqdm
 
 import easyidp as idp
 
@@ -416,41 +417,45 @@ def __init__(self):
         self.matrix_inv = None
 
 
-def filter_closest_img(p4d, img_dict, plot_geo, dist_thresh=None, num=None):
-    """[summary]
+def _sort_img_by_distance_one_roi(recons, img_dict, plot_geo, cam_pos, distance_thresh=None, num=None):
+    """Sort the back2raw img_dict results by distance from photo to roi
 
     Parameters
     ----------
+    recons: idp.Metashape or idp.Pix4D
+        The reconsturction project class
     img_dict : dict
-        The outputs dict of geo2raw.get_img_coords_dict()
+        One ROI output dict of roi.back2raw()
+        e.g. img_dict = roi.back2raw(ms) -> img_dict["N1W1"]
     plot_geo : nx3 ndarray
         The plot boundary polygon vertex coordinates
     num : None or int
         Keep the closest {x} images
-    dist_thresh : None or float
+    distance_thresh : None or float
         If given, filter the images smaller than this distance first
 
     Returns
     -------
     dict
-        the same structure as output of geo2raw.get_img_coords_dict()
+        the same structure as output of roi.back2raw()
     """
     dist_geo = []
     dist_name = []
-    
+
     img_dict_sort = {}
-    for img_name, img_coord in img_dict.items():
 
+    for img_name in img_dict.keys():
         xmin_geo, ymin_geo = plot_geo[:,0:2].min(axis=0)
         xmax_geo, ymax_geo = plot_geo[:,0:2].max(axis=0)
+
         xctr_geo = (xmax_geo + xmin_geo) / 2
         yctr_geo = (ymax_geo + ymin_geo) / 2
 
-        ximg_geo, yimg_geo, _ = p4d.img[img_name].cam_pos
+        ximg_geo, yimg_geo, _ = cam_pos[img_name]
 
         image_plot_dist = np.sqrt((ximg_geo-xctr_geo) ** 2 + (yimg_geo - yctr_geo) ** 2)
 
-        if dist_thresh is not None and image_plot_dist > dist_thresh:
+        if distance_thresh is not None and image_plot_dist > distance_thresh:
             # skip those image-plot geo distance greater than threshold
             continue
         else:
@@ -466,5 +471,41 @@ def filter_closest_img(p4d, img_dict, plot_geo, dist_thresh=None, num=None):
 
     dist_geo_idx = np.asarray(dist_geo).argsort()[:num]
     img_dict_sort = {dist_name[idx]:img_dict[dist_name[idx]] for idx in dist_geo_idx}
-
-    return img_dict_sort
+
+    return img_dict_sort
+
+
+def sort_img_by_distance(recons, img_dict_all, roi, distance_thresh=None, num=None):
+    """Advanced wrapper of sorting back2raw img_dict results by distance from photo to roi
+
+    Parameters
+    ----------
+    recons: idp.Metashape or idp.Pix4D
+        The reconsturction project class
+    img_dict_all : dict
+        All output dict of roi.back2raw(...)
+        e.g. img_dict = roi.back2raw(...) -> img_dict
+    roi: idp.ROI
+        Your roi variable
+    num : None or int
+        Keep the closest {x} images
+    distance_thresh : None or float
+        Keep the images closer than this distance to ROI.
+
+    Returns
+    -------
+    dict
+        the same structure as output of roi.back2raw()
+    """
+    cam_pos = recons.get_photo_position(to_crs=roi.crs)
+
+    img_dict_sort_all = {}
+    pbar = tqdm(roi.keys(), desc=f"Filter by distance to ROI")
+    for roi_name in pbar:
+        sort_dict = _sort_img_by_distance_one_roi(
+            recons, img_dict_all[roi_name], roi[roi_name], 
+            cam_pos, distance_thresh, num
+        )
+        img_dict_sort_all[roi_name] = sort_dict
+
+    return img_dict_sort_all

easyidp/roi.py

@@ -2,6 +2,7 @@
 import pyproj
 import warnings
 import numpy as np
+from tqdm import tqdm
 from copy import copy as ccopy
 from shapely.geometry import Point, Polygon
 from pathlib import Path
@@ -269,13 +270,16 @@ def get_z_from_dsm(self, dsm, mode="face", kernel="mean", buffer=0, keep_crs=Fal
             global_z = None
 
         # convert CRS if necessary
-        if self.crs.name != dsm.header["crs"].name and not keep_crs:
+        if self.crs.name == dsm.header["crs"].name:
+            poly_dict = self.id_item.copy()
+        elif self.crs.name != dsm.header["crs"].name and not keep_crs:
             self.change_crs(dsm.header["crs"])
             poly_dict = self.id_item.copy()
         else:
             poly_dict = idp.shp.convert_proj(self.id_item, self.crs, dsm.header["crs"])
 
-        for key, poly in poly_dict.items():
+        pbar = tqdm(poly_dict.items(), desc=f"Read z values of roi from DSM [{dsm.file_path.name}]")
+        for key, poly in pbar:
             # only get the x and y of coords
             poly = poly[:, 0:2]