Shoreline extraction issues

[conlin-matt]

@kvos,

I am (sometimes) getting a strange error message when I call "SDS_shoreline.extract_shorelines(metadata,settings)", where the error message is "ValueError: LineStrings must have at least 2 coordinate tuples" (screenshot attached). Sometimes the function works without error, and sometimes I get this message. I cannot seem to detect a pattern of when it errors (i.e. only for certain satellites or dates, etc.), other than that it only seems to happen when I apply the function to a relatively long interval (e.g. 15 yr). For example, I can call the function, for the same satellites and end date as established in "metadata", for 3 months and it will work fine, but will error if I try it for 10 yr. As the screenshot shows, the function gets all the way to 100% for each satellite before it errors out. For reference, I'm running MacOS 10.14 with Python 3.7. Any ideas?

[kvos]

I have corrected this bug in the latest update to the toolbox...if you pull the latest changes it shouldn't happen any more.

The problem was that if you accept an image that has no mapped shoreline it would then crash when trying to write a .geojson file as there are no coordinates for that shoreline.

[conlin-matt]

Yep, this issue seems to be fixed. Thanks!

[kvos]

ok perfect. Did you also figure out why it could not map shorelines (the message: "Could not map shoreline for this image" that appears in your terminal)?

[conlin-matt]

No, that issue has stumped me so far. The warning only seems to appear for S2 images. When the "find_wl_contours2" function (i.e. line 734-735 of SDS_shoreline.py) is run on an image which results in the error, I get the ValueError shown in the screenshot. So it looks like it has something to do with an inability to detect ranges for different classes in the thesholding maybe? I should also note that the older, less robust thresholding method (implemented in find_wl_contours1) works on the image that generates the error.

[kvos]

could you post the image @conlin-matt ? I think the error is that Otsu's threholding function doesn't work if there are NaNs in the histogram, but I don't see how you could get NaNs in there... what you can do also is that you can put pdb.set_trace() before it calls threshold_otsu and it will stop the code there and then you can see what's in int_all and check if there are any NaNs or inf values...

[conlin-matt]

Github doesn't support uploads of TIFFs, so I haven't uploaded the image. After some debugging this morning, here's what I found:

You are correct that this issue appears when "int_all" has NaNs in it. This happens when the SWIR1 and green bands of the image both contain a value of 0 in the same pixel (i.e. im_ms[:,:,4] and im_ms[:,:,1] in your code both have 0s for the same pixels). Here is what happens when this occurs:

1. SDS_shorelines.extract_shorelines calls SDS_shorelines.find_wl_contours2 (line 734-735 of SDS_shorelines)
2. SDS_shorelines.find_wl_contours2 calls SDS_tools.nd_index (line 248 of SDS_shorelines)
3. In SDS_tools.nd_index you have a calculation on lines 181-182 where you np.divide the two bands. But when you run np.divide(0,0) the result is a nan. So the variable im_mwi which results from the call to this function has NaNs in it, which propagate into int_all, causing the error.

Perhaps a fix to this could be to replace nans with zeros in im_mwi?

[kvos]

oh yeah, when trying to divide 0 by 0 it returns a NaN, good catch @conlin-matt I didn't think an image could have 0s at the same pixel for different bands. could you post a JPEG of that image (generated with SDS_preprocess.save_jpg(metadata, settings) ? the image must have some very dark areas...

I think replacing nans with zeros should be fine... another option could be to add those pixels to the cloud_mask so they are ignored

[conlin-matt]

Yep, looks like there are definitely some black spots in the image...

[kvos]

the black spots located on the clouds are part of the cloud mask (I know the cloud mask is terrible for Sentinel-2 images!)... what I am more worried about is what's happening on the top-left corner of the image, it seems that the pixels are a bit funny (blueish colours) at the edge of the no data region. My guess would be that the 0 division is happening there, if you don't mind sending me the coordinates of your polygon I can also investigate that 2019-01-01 S2 image.

[conlin-matt]

Gotcha. Polygon coordinates are below:

polygon = [[[-80.680891,28.699524], [-80.587325,28.526645], [-80.510662,28.537539], [-80.620483,28.698629], [-80.680891,28.699524]]]

[kvos]

the image below shows in red the pixels that have an intensity of 0 for ALL BANDS.

I added to the cloud mask the pixels that have 0s in the G, NIR and SWIR bands in S2 images. I just pushed this fix into the development branch. Let me know if it works.

Also @conlin-matt , the area that you are mapping is quite large, if you see that processing the shorelines is too slow you can split the polygon in 2/3 smaller polygons to speed up the processing (more images to manually check though).

[conlin-matt]

Yep, that fix seems to be working for me. Glad we were able to figure these things out!

[kvos]

cool! did you run into this issue also with Landsat images or it was only S2?

[conlin-matt]

It only ever popped up with S2 images for me.

[kvos]

issue solved