My bad for unclear explanation. I think it's the latter. I want to get angles for 1, 2, and 3 in this drawing.

Originally posted by @Ken-Kuroki in #12 (comment)

Once init the dom = idp.GeoTiff(), still need idp.geotiff.pixel2geo(roi, dom.header) & idp.geotiff.geo2pixel(roi, dom.header) to calculate. It should be simplified to dom.geo2pixel(roi) and dom.pixel2geo(roi)

seems current code for shp reading is a little bit slow

The tqdm progressbar also support multiproceeing by imap

>>>img_coord
array([[7214.31561958, 3741.17729258],
       [6090.04392943, 3062.91735589],
       [6770.00648193, 1952.53149042],
       [7901.26625814, 2624.35137412],
       [7214.31561958, 3741.17729258]])
>>> imarray, offsets = idp.cvtools.imarray_crop(img_np, img_coord)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_17496/622690088.py in <module>
----> 1 imarray, offsets = idp.cvtools.imarray_crop(img_np, img_coord)

e:\github\easyidp\easyidp\cvtools.py in imarray_crop(imarray, polygon_hv, outside_value)
     79         # has 3 dimensions
     80         # e.g. DOM with RGB or RGBA band, other value outside changed alpha layer to 0
---> 81         roi_clipped = imarray[roi_top_left_offset[1]:roi_max[1], roi_top_left_offset[0]:roi_max[0], :]
     82 
     83         rh = roi_clipped.shape[0]

TypeError: slice indices must be integers or None or have an __index__ method

fix method:

>>> img_coord_int = img_coord.astype(int)
>>> imarray, offsets = idp.cvtools.imarray_crop(img_np, img_coord_int)

I wonder if you would consider eventually adding support for OpenDroneMap reconstructions. This could be useful for users aiming to build an entirely open-source workflow.

ODM internally uses OpenSfM for photogrammetry. OpenSfM saves reconstruction data (e.g., estimated camera locations and orientations) in a file reconstruction.json.

Thanks for considering, and thanks again for making and supporting this excellent software!

HI @HowcanoeWang I have managed to get the reverse calculation working with my Metashape imagery and the ROI from my shapefile. What would be the best way to export my reverse calculated images (.jpg) so that I can process them with Easy PCC? Would it be better to rotate the images, or add a mask based on the ROI? E.g here is an example plot which I've added a black mask to. Would this type of image be ok to process with Easy PCC? Thanks again for your help

Originally posted by @UQDannySmith in #50 (comment)

Current outputs:

| [0] id   | [1] line_id   |
|----------|---------------|
|          | Koshi_1       |
|          | Koshi_2       |
|          | Koshi_3       |
|          | Koshi_1       |
|          | Koshi_1       |
| ...      | ...           |
|          | Koshi_1       |
|          | Koshi_1       |
|          | Koshi_1       |
|          | Koshi_1       |
|          | Koshi_1       |

expected：

   | [0] id   | [1] line_id   |
   |----------|---------------|
1  |          | Koshi_1       |
2  |          | Koshi_2       |
3  |          | Koshi_3       |
4  |          | Koshi_1       |
5  |          | Koshi_1       |
...| ...      | ...           |
6  |          | Koshi_1       |
7  |          | Koshi_1       |
8  |          | Koshi_1       |
9  |          | Koshi_1       |
10 |          | Koshi_1       |

Also support using -1 and "index" to replace current None to using line number as id

The chunk_id is not the corrected order in metashape, chunk.label is the most correct.

Also roi.back2raw

>>> rotated
array([[[0.31690149, 0.38752052, 0.14294289],
        [0.3404107 , 0.41297721, 0.16013014],
        [0.38506028, 0.4331061 , 0.17238887],
        ...,
        [0.67001947, 0.70273325, 0.44726555],
        [0.71262498, 0.74353751, 0.49159185],
        [0.63883647, 0.66216885, 0.41865424]],
>>> rotated_geotiff_imarray, _ = idp.cvtools.imarray_crop(rotated, roi_geotiff_offseted.astype(int))
>>> plt.imshow(rotated_geotiff_imarray)

>>> rotated_geotiff_imarray[:,:,0]
array([[0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       ...,
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0]], dtype=uint8)

Can be fixed only if the input image type is uint8

>>> rotated = (rotated * 255).astype(np.uint8)

Reason:

Currently, call get_photo_position will cost some time to re-calculate.

Most of my photogrammetry (Metashape) projects involve processing photos that are nested within subfolders; for example:

images/
├── 100MEDIA/
│   ├── DJI_0001.jpg
│   ├── DJI_0002.jpg
│   ├── DJI_0003.jpg
│   └── ...
├── 101MEDIA/
│   ├── DJI_0001.jpg
│   ├── DJI_0002.jpg
│   ├── DJI_0003.jpg
│   └── ...
└── 102MEDIA/
    ├── DJI_0001.jpg
    ├── DJI_0002.jpg
    ├── DJI_0003.jpg
    └── ...

But when I try to load such a Metashape project in EasyIDP, I get: IndexError: Index [0] out of range (0, 0).

To reproduce it, you can download this example metashape project file (which is based on this image set) and run:

import easyidp as idp
from pathlib import Path
msproj = Path("path/to/project/example-multifolder.psx")
ms = idp.Metashape(msproj, chunk_id=0)

The full error returned is:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
/tmp/ipykernel_481076/3594889881.py in <module>
----> 1 ms = idp.Metashape(msproj, chunk_id=0)

~/Downloads/EasyIDP-2.0/easyidp/metashape.py in __init__(self, project_path, chunk_id)
    105         self.photos = self.photos
    106 
--> 107         self.open_project(project_path, chunk_id)
    108 
    109     def __repr__(self) -> str:

~/Downloads/EasyIDP-2.0/easyidp/metashape.py in open_project(self, project_path, chunk_id)
    172         if project_path is not None:
    173             self._open_whole_project(project_path)
--> 174             self.open_chunk(self.chunk_id)
    175         else:
    176             if chunk_id is not None:

~/Downloads/EasyIDP-2.0/easyidp/metashape.py in open_chunk(self, chunk_id, project_path)
    244 
    245         if chunk_id in self._project_chunks_dict.keys():
--> 246             chunk_content_dict = read_chunk_zip(
    247                 self.project_folder,
    248                 self.project_name,

~/Downloads/EasyIDP-2.0/easyidp/metashape.py in read_chunk_zip(project_folder, project_name, chunk_id, skip_disabled, return_label_only)
   1176         camera_meta, marker_meta = _decode_frame_xml(frame_xml_str)
   1177         for camera_idx, camera_path in camera_meta.items():
-> 1178             chunk_dict["photos"][camera_idx]._path = camera_path
   1179             # here need to resolve absolute path
   1180             # <photo path="../../../../source/220613_G_M600pro/DSC06035.JPG">

~/Downloads/EasyIDP-2.0/easyidp/__init__.py in __getitem__(self, key)
     72                 return self.id_item[key]
     73             else:
---> 74                 raise IndexError(f"Index [{key}] out of range (0, {len(self.item_label)})")
     75         elif isinstance(key, str):  # index by photo name
     76             if key in self.item_label.keys():

IndexError: Index [0] out of range (0, 0)

The error occurs even if the photo file names are not duplicated across the folders (which is true of the example I ran and linked to here). Ideally, EasyIDP would tolerate multiple subfolders and an arbitrary folder nesting depth (at least to 5 levels).

current support get_z_from_dsm using buffer, need adding buffer for ROI

The info can be found in the *.p4d file as xml format:

<images>
            <image path="Y:/.../maize_tanashi_3NA_20190729_Ins1Rgb_30m/DJI_0954.JPG" group="group1" enabled="true">
                <camera name="FC550_DJIMFT15mmF1.7ASPH_15.0_4608x3456" id="0"/>
                <time>2019:07:29 09:01:11</time>
                <gps alt="72.50900000000000034" lat="35.73757452777778099" lng="139.5406209722222286"/>
                <xyz x="0" y="0" z="0"/>
                <accuracyXY>5</accuracyXY>
                <accuracyZ>10</accuracyZ>
                <ori yaw="38.1" pitch="0.1" roll="0" perpCamera="true"/>
                <accuracyOri yaw="15" pitch="10" roll="10"/>
            </image>

            ...

for testing module, no need to manual specify the data path

The core reason caused this differences:

img_np = plt.imread(ms.photos["DJI_0479"].path)
VS
img_np = plt.imread(p4d.photos["DJI_0483.JPG"].path)

import easyidp as idp
print(idp.__version__)

2.0.0.dev3

ms = idp.Metashape("/Users/benweinstein/Dropbox/Weecology/everglades_species/easyidp/Hidden_Little_03_24_2022.psx", chunk_id=0)
roi = idp.ROI("/Users/benweinstein/Dropbox/Weecology/everglades_species/easyidp/example_bird.shp") 
img_dict = roi.back2raw(ms)

[shp][proj] Use projection [WGS 84] for loaded shapefile [example_bird.shp]

[shp] read shp [example_bird.shp]:   0%|          | 0/1 [00:00<?, ?it/s]
[shp] read shp [example_bird.shp]: 100%|##########| 1/1 [00:00<00:00, 212.31it/s]Traceback (most recent call last):
  Python Shell, prompt 22, line 5
  File "/Users/benweinstein/.conda/envs/DeepTreeAttention/lib/python3.8/site-packages/easyidp/roi.py", line 401, in back2raw
    out_dict = recons.back2raw(self, **kwargs)
  File "/Users/benweinstein/.conda/envs/DeepTreeAttention/lib/python3.8/site-packages/easyidp/metashape.py", line 259, in back2raw
    save_path = None
  File "/Users/benweinstein/.conda/envs/DeepTreeAttention/lib/python3.8/site-packages/easyidp/metashape.py", line 209, in back2raw_crs
    local_coord = self._world2local(self._crs2world(points_xyz))
  File "/Users/benweinstein/.conda/envs/DeepTreeAttention/lib/python3.8/site-packages/easyidp/metashape.py", line 112, in _crs2world
    return convert_proj3d(points_np, self.crs, self.world_crs)
  File "/Users/benweinstein/.conda/envs/DeepTreeAttention/lib/python3.8/site-packages/easyidp/metashape.py", line 1110, in convert_proj3d
    out = np.vstack([lat_m, lon_m, alt_m]).T
builtins.IndexError: index 2 is out of bounds for axis 1 with size 2

Looks like we are missing the z dimension.

Currently need to run one by one until all finished, hoping to give an empty list for one-time modificatiion

KeyError: "Meet with duplicated key [3058#1_5] for current shapefile, ...

idp.reconstruct.filter_closest_img()

    Also, when DSM did not cover all ROI, will return 

[x, y, nan]

also need to notice

Originally posted by @HowcanoeWang in #47 (comment)

Original point cloud import to ArcGIS:

Cropped point cloud import to ArcGIS:

The origina function built in caas_lite in EasyIDP v1.0, aimed to crop geotiff to several parts conveniently

Current single demo involves too many functions, splited them to small docs, focusing on simple function and its example

one figure contains all backward results preview, sorted by the distance to image center

Also, should support obtaining the imarray of backward ROI on raw images.

remove local files and redownload from Internet

To fix the poor quality DOM and also match with the raw2dom results

Usage:

top10 = roi[0:10]

This can decrease the workload for demo with thousands of roi that not used

Discussed in #71

^{Originally posted by youngdjn January 29, 2023}
Thank you for making this excellent package! It is incredibly useful. I have a use case where I need to use the mode="point" option for ROI.get_z_from_dsm(), but there appears to be a bug with this option. When I run the following:

import easyidp as idp
lotus = idp.data.Lotus()
roi = idp.ROI(lotus.shp, name_field = "plot_id")
roi.get_z_from_dsm(lotus.metashape.dsm, mode="point")

I get:

ValueError: all the input arrays must have same number of dimensions, but the array at index 0 has 2 dimension(s) and the array at index 1 has 3 dimension(s)

I need the mode="point" option because I am trying to extract images of the crowns of individual trees, which have a lot of height variation. When I use with mode="face", the projected ROI does not follow my feature quite right (though it's close, and generally still usable).

(I drew the "pacman" shape because I was hoping the central vertex would track the top of the tree.)

>>> dom = idp.GeoTiff("xxxx.tiff")
>>> dom
<matplotlib.pyplot.axes>

For large dom (> 2000x2000), use the subpixels, and only create a 200x200 preview (make grid point lists and point_query for pixel values of these points).

To load and batch process several reconstruction projects

Discussed in #39

roi = easyidp.ROI("roi.shp")
roi.get_z_from_dsm("dsm.tif")

Discussed in #75

import easyidp as idp
from pathlib import Path
msproj = Path("/path/to/project/example-multifolder.psx")
ms = idp.Metashape(msproj, chunk_id=0)
ms.photos[57].label

for i in range(len(ms.photos)):
    ms.photos[i].label = "newlabel" + str(i) 
ms.photos[57].label

img_dict_ms = roi.back2raw(ms)
img_dict_ms

{'1': {'100MEDIA-DJI_0165': array([[4741.6161679 ,  684.46603374],
         [4735.11907978,  712.14368673],
         [4718.00310275,  728.6865289 ],
         [4722.27333877,  728.17479155],
         ...

ms.show_roi_on_img(img_dict_ms, "1", "100MEDIA-DJI_0165")

ms.show_roi_on_img(img_dict_ms, "1", "newlabel57")

I'm sure i'll have questions in the next few weeks, but just making an issue saying thanks. This is exactly what we needed. As a python package developer for machine learning, https://deepforest.readthedocs.io/, I can see this all looks great. I'll be sure to share with others, i've heard many people wondering about this process.

Seems both cloudcompare and ArcGIS does not support this function very well

I would like to request an option for ROI.back2raw(..., save_folder="...") to save cropped images without masking them to the projected polygon; just crop the images to the bounding box of the projected polygon (ideally with a buffer option too). For example, for one image, the projected polygon looks like this:

but the image saved when using the save_folder option gets masked to this:

I would prefer for the saved image file to not be masked; just save the unmasked image for the full crop extent. It would be also nice if an optional bounding box buffer width (in pixels) could be applied before cropping the image. The buffer that is applied in the preview image (first image above) is nice.

Thanks for considering this enhancement.

I could replicate the error you mentioned here

>>> import easyidp as idp
>>> lotus = idp.data.Lotus()
>>> ms = idp.Metashape(lotus.metashape.project, chunk_id=0) 
>>> ms._world2crs(ms._local2world(ms.photos[0].transform[0:3, 3])) 
array([139.54053245,  35.73458169, 130.09433649])

But if you changed the ms.crs to geotiff and then do the same thing:

>>> import pyproj                                   
>>> ms.crs = pyproj.CRS.from_epsg(32654)              
>>> ms._world2crs(ms._local2world(ms.photos[0].transform[0:3, 3]))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "D:\OneDrive\Program\GitHub\EasyIDP\easyidp\metashape.py", line 99, in _world2crs
    return convert_proj3d(points_np, self.world_crs, self.crs)
  File "D:\OneDrive\Program\GitHub\EasyIDP\easyidp\metashape.py", line 1047, in convert_proj3d
    return out[0, :]
UnboundLocalError: local variable 'out' referenced before assignment

Temporary solution:

Set ms.crs = pyproj.from_epsg(4326) first and then get photo location, and transform to GIS crs by

pos = {p.label: ms._world2crs(ms._local2world(p.transform[0:3, 3]))   for p in ms.photos}
transformer = pyproj.Transformer.from_proj(pyproj.from_epsg(4326), pyproj.from_epsg(32654))

out = {}
for k,v in pos.items():
    # here the lon and lat need to be reversed for pyproj package
    transformed = transformer.transform(v[1], v[0], v[2])
    out[k] = np.asarray(transformed).T

Detailed explanation:
Ths crs in the Metashape is created from photos' meta_info, often epsg4326 with lon and lat, and the software will "display" in this coords.

local -> world -> crs(photo_meta)

When you export outputs from Metashape, you can specify another CRS for geotiff (in this case is epsg:32654)

local -> world -> crs(photo_meta) -> geotiff

Probably in the next release, the ms.crs -> ms.photo_meta_crs, and then ms.crs -> geotiff.crs

and you can get photo position at current geotiff crs by

>>> ms.get_photo_position(photos[0], crs=geotiff.crs)
array([[ 368030.75661044, 3955476.94877388, 129.8508284]])

Originally posted by @HowcanoeWang in #12 (reply in thread)

The current algorithm just moves points outside the image boundary along giving direction, but expect to crop by image boundary.

Current print style now showing the correct ROI key and looks ugly:

       [ 484585.84831165, 3862266.69649805]]), 116: array([[ 484586.95917975, 3862264.48670892],
       [ 484586.27618307, 3862265.79420158],
       [ 484586.74248551, 3862266.01187889],
       [ 484587.38897386, 3862264.72551779],
       [ 484586.95917975, 3862264.48670892]]), 117: array([[ 484588.07918062, 3862262.27468665],
       [ 484587.39618362, 3862263.58217923],
       [ 484587.86248614, 3862263.79985659],
       [ 484588.50897668, 3862262.51460452],
       [ 484588.07918062, 3862262.27468665]]), 118: array([[ 484589.18991156, 3862259.98283581],
       [ 484588.50691422, 3862261.2903283 ],
       [ 484588.97321497, 3862261.50689677],
       [ 484589.6197077 , 3862260.22275374],
       [ 484589.18991156, 3862259.98283581]]), 119: array([[ 484590.32812571, 3862257.73529679],
       [ 484589.64512804, 3862259.04278919],
       [ 484590.11143074, 3862259.26046666],
       [ 484590.75792192, 3862257.97521477],
       [ 484590.32812571, 3862257.73529679]])}

dom_parts = roi.crop(dom)

Should return {"roi_key": <GeoTiff object>}

And get the value by GeoTiff.values

Discussed in #33

import easyidp as idp
project_path = "/data/2017_tanashi_lotus/170531.Lotus.psx"
dsm_path = "/data/2017_tanashi_lotus/170531.Lotus.outputs/170531.Lotus_dsm.tif"
ms = idp.Metashape(project_path, chunk_id=0)
roi = idp.ROI("/data/2017_tanashi_lotus/plots.shp")
roi.get_z_from_dsm(dsm_path)
out = ms.back2raw(roi)

if RAM >> GeoTiff.filesize and len(roi) > 100:
    full load GeoTiff, and delete the variable after finish.

Have test fulling loading into memory, DOM.dtype(uint8) == DSM.dtype(float32), 20w x 20w dom requires 91GB RAM. Not applicable.

The optimal choice:

read one line tail for muitiple roi if applicable.

No need duplicately read geotiff for each roi.

> ROI.get_z_from_dsm():

1. calculate each roi and its bbox.
2. bbox to geotiff tail line & colume
3. revert previous index, to each geotiff tail, which columns are which roi and roi corresponding parts
4. For loop, for each geotiff tail, read to numpy, and crop & save to correspoing roi parts
    1. If one roi is completed loaded, return cropped results for calculation.
    2. delete the completed roi in memory.
    3. continue for the next tail.