Trouble with the split_raster command

[jeffgillan]

Hi DeepForest. Thanks for making this python library and helping to simplify deep learning and CNNs for those of us not well versed in pytorch. My name is Jeff Gillan and I'm a researcher/data scientist at University of Arizona. I am also associated with the NSF funded Open Forest Observatory project. Over the last few months I've been playing with DeepForest and trying to identify trees in some of my drone imagery (~ 1 cm gsd). The prebuilt model kind of works but it seems some fine-tuning is necessary to really make it work well.

I have a shapefile of tree crown labels and have successfully converted them to the annotation data frame using the shapefile_to_annotations. The conversion from geographic coordinates to pixel coordinates seems to work. The problem comes when I use:

output_crops = split_raster(path_to_raster=train_image_path,
                            annotations_file=annotation_path,
                            base_dir=crop_dir,
                            patch_size=1100,  
                            patch_overlap=0.25, 
                            allow_empty=False)

At first it appears that everything was done correctly. The command takes the large geotif and splits it into many equally sized smaller images along with the associated annotation. But upon further inspection, the pixel coordinates of the cropped image annotation do not make any sense. The bounding boxes do not appear to be where the should be, (ie, on top of the tree crowns).

I was oblivious to this for a long time and assumed I was making a mistake during model.trainer.fit(model). But I think I tracked it down to the split_raster step.

Here is an cropped image output from model.evaluate that hopefully helps illustrate my problem.

I assume the orange boxes are the manual labels that I provided and the blue boxes are the prediction. I assure you, those orange boxes are not where they should be.

Any ideas what could be the problem here? I'm totally stumped. If you want to see my entire code, it is here Thanks for your help.

[bw4sz]

Thanks for this discussion. At the start, I want to underline that we never think of DeepForest as a zero-shot tool for tree detection globally. Retraining will always be the right choice.

Can you confirm you are using the latest DeepForest version?

There must be something specific about your use case, as many users use split_raster all the time in the way you described. I happened to run it myself like this yesterday. Let's try to unravel what's happening with your data. Can you share a single image and annotation file? Then we can visualize each step and figure out what's happening.

The general steps are to use deepforest.visualize.plot_predictions to confirm that the initial annotations look correct.

Then let's try patch_size = 0, to see if its the overlap.

I suspect that this is a flip_y_axis issue, in what program are you looking at the image above? Numpy/matplotlib/qgis all use very slightly different axis conventions. But that's just a guess.

[jeffgillan]

Thank you for the quick response.

I'm using deepforest 1.3.3 from pypi and pre-built model 1.0.0 from your github repo
I'm running this code in a jupyter notebook hosted on Cyverse

I started by creating tree crown label in QGIS (following your method). I digitized rectangular bounding boxes (n=150) on top of a drone orthomosaic geotif (RGB; EPSG:4326 - WGS 84) and saved them as a shapefile (EPSG:4326 - WGS 84).

Next, I converted the shapefile to annotation and used the same orthomosaic geotif as the rgb argument

df = shapefile_to_annotations(shapefile=shp_path, rgb=train_image_path, geometry_type='bbox', savedir=savedir)

This produced a dataframe with 6 columns (image_path | xmin | ymin | xmax | ymax | label) and 150 rows
How can I prove that these image coordinates are correct spatially with respect to my orthomosaic?

Next, I did the split_raster command. I tried the setting overlap=0. This did not improve my problem.

Attached is a single image (.png) that was created from the split_raster command and the annotation file (csv) that was output from the split_raster

1a_g2_ortho.csv

Sorry, but I don't understand how to use deepforest.visualize.plot_predictions. Perhaps you could provide some python code for me to use it correctly.

The image I provided in my first message was a png what was output from your deepforest command.
results = model.evaluate(training_file, os.path.dirname(training_file), iou_threshold = 0.4, savedir= save_dir)

I did not use any python software to look at the this first image. It was a .png that I simply clicked on and opened on my computer.

Something seems to be going wrong in the either the split_raster or maybe even the shapefile_to_annotations command. The annotations don't seem to line up on the correct place on the orthomosaic. No idea why.

[jeffgillan]

We could go back in forth forever on this topic. It may be more useful to do a Zoom call together. My email is [email protected]

[jeffgillan]

Do y-axis image coordinates start from the top or start from the bottom?

[jeffgillan]

Because I am struggling to execute annotation of geospatial data (qgis, shapefiles), perhaps it makes sense to do the tree crown labeling on my many non-georeferenced jpeg drone images. In this way I could skip any conversion from geographic to image coordinates. Would you recommend labeling of non-georeferenced drone images? Do you have a suggested workflow to do this? I probably just need to get the annotation data into your format (image_path; xmin, ymin, xmax, ymax, label)? Have you used Label Studio to do this? Thank you.

[bw4sz]

Hey @jeffgillan, let's get this done, I sense there isn't anything wrong here, just the docs need to be better.

your command looks fine here

"This produced a dataframe with 6 columns (image_path | xmin | ymin | xmax | ymax | label) and 150 rows
How can I prove that these image coordinates are correct spatially with respect to my orthomosaic?"

you can confirm they are correct by using https://deepforest.readthedocs.io/en/latest/source/deepforest.html#deepforest.visualize.plot_predictions.

Can you slow down, provide full code, annotations and image data and we can go through this. Let see the plots and make this reproducible.

[jeffgillan]

@bw4sz Can you please suggest some python snippet to use your plot_predictions tools? I'm weak at python and don't know the exact code needed to use your visualization tool. I have spent many hours trying with no success.

My entire code is here https://github.com/ua-datalab/Geospatial_Workshops/blob/main/notebooks/deepforest_2.ipynb

Thank you

[jeffgillan]

@bw4sz I've managed to cobble some python together and have made some progress. I now believe that the root of my problem may be occurring in the shapefile_to_annotations step.

My evidence:
When I bring the orthomosaic and the training labels shapefile into my jupyter notebook, everything seems to display correctly. My labels are directly on top of the trees.

After I do the shapefile_to_annotations step, I used your plot_predictions tool to display the annotations on top of the orthomosaic. It does not look correct. It appears that all of the boxes have shifted northward in the y axis.

Any ideas?

[jeffgillan]

@bw4sz I have done a bit of testing to figure out the behavior of shapefile_to_annotation. It seems to do the 'shift' when I use an orthomosaic and shapefile that do not have projections. I am using geographic coordinate system WGS84 for the data. When I use data (orthomosaic and shapefile) with projections (e.g., utm) the tool seems to work as it should. Does seem plausible? thank you.

[bw4sz]

I really cannot help you unless you include code and data. That is really the only way you can get assistance with something specific. Everything else is just guesses. Are we talking about the same data and code as above? Show me the code that produces the result with WGS84, include the annotations and ortho in this thread. That way we can tell if this is a genuine bug, or something specific to your dataset.

[jeffgillan]

Here is the code, image data, and shapefile to reproduce the 'shift' error that I seem to be seeing when using data without map projections.

code
Geotiff orthomosaic
shapefile

[bw4sz]

Okay, I can confirm the error with supplied data. When using degrees, the transform does not work as expected. Current workaround is to use a utm or other meters based transform.

from deepforest import main
from deepforest.utilities import shapefile_to_annotations
from deepforest.visualize import plot_predictions
import rasterio
from matplotlib import pyplot as plt

m = main.deepforest()
m.use_release()

shp = shapefile_to_annotations(rgb="/Users/benweinstein/Downloads/1a_g2_ortho.tif",
                                shapefile="/Users/benweinstein/Downloads/srer_train.shp")

image = rasterio.open("/Users/benweinstein/Downloads/1a_g2_ortho.tif").read()
plotted_image = plot_predictions(image=image, df=shp, thickness=3, color=(255, 0, 0))
plt.imshow(plotted_image)
plt.show()

Tiny note: https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.13472 is the correct citation, from the notebook, I can see a different paper.

[jeffgillan]

Thank you Ben for the confirmation, work around, and updated reference.

On another note, I have successfully fine-tuned the original model using my own tree crown labels. Just a few labels significantly improved prediction so I'm hoping many more labels will really dial in the accuracy. Thanks again for DeepForest!

[bw4sz]

What is strange is that I cannot recreate your error on a new dataset, trying to decipher why.

def test_shapefile_to_annotations_latlon(tmpdir):
    image_path = get_data("OSBS_029.tif")
    annotation_path = get_data("OSBS_029.csv")
    df = pd.read_csv(annotation_path)
    gdf = utilities.boxes_to_shapefile(df, root_dir=os.path.dirname(image_path), projected=True)
    gdf.to_crs("EPSG:4326", inplace=True)
    gdf.to_file("{}/annotations.shp".format(tmpdir))
    
    with rasterio.open(image_path) as src:
        # Reproject image to EPSG:4326
        transform, width, height = rasterio.warp.calculate_default_transform(src.crs, "EPSG:4326", src.width, src.height, *src.bounds)
        kwargs = src.meta.copy()
        kwargs.update({
            'crs': "EPSG:4326",
            'transform': transform,
            'width': width,
            'height': height
        })
        
        # Create a new reprojected image
        reprojected_image_path = "{}/reprojected_image.tif".format(tmpdir)
        with rasterio.open(reprojected_image_path, 'w', **kwargs) as dst:
            for i in range(1, src.count + 1):
                rasterio.warp.reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=transform,
                    dst_crs="EPSG:4326",
                    resampling=rasterio.warp.Resampling.nearest
                )
    
    shp = utilities.shapefile_to_annotations(shapefile="{}/annotations.shp".format(tmpdir), rgb=reprojected_image_path, savedir=tmpdir)
    
    # plot the image
    image = rio.open(reprojected_image_path).read()
    plotted_image = plot_predictions(image=image, df=shp, thickness=3, color=(255, 0, 0))
    plt.imshow(plotted_image)
    plt.show()
    assert shp.shape[0] == 2

[bw4sz]

I suspect your orthophoto is not properly referenced. That it was assigned but not transformed in some prior mode. I cannot say for certain, but when you look at in QGIS, I don't like the askew angle despite the fact that the coordinate systems are correct.

This leads me to believe there was a forced reference system change without transformation. I'm going to leave this issue open, but the fact that i cannot recreate outside the data you provide suggests there is not a bug here. I'm continuing to investigate with more independent data and will post.

[jeffgillan]

Thanks for your further inquiry. I have done some testing myself. The orthomosaics are outputs of Agisoft Metashape. Yesterday, I output 2 orthomosaics, one with EPSG:4326 and one with UTM coordinates. The orthomosaics were brought into QGIS and both were labeled with polygon shapefiles. The CRS of the shapefiles matched the CRS of the orthomosaics. The orthomosaic and shapefile with EPGS:4326 again showed the same 'shift' when using shapefile_to_annotations. The orthomosaic and shapefile with UTM worked correctly. I'm not really sure what's going on, but at least I know a fix.

[bw4sz]

Okay, i'm starting to get a handle on this.

[bw4sz]

I think we can close this in favor of a dedication issue and bug fix. #695