kvos / CoastSat

Global shoreline mapping tool from satellite imagery
http://coastsat.space
GNU General Public License v3.0
696 stars 252 forks source link

Question: Why is ` n_intersect[i] = len(xy_rot[0, :])` used instead of `n_intersect[i] = np.sum(~np.isnan(xy_rot[0, :]))` in `compute_intersection_QC` #500

Closed 2320sharon closed 4 months ago

2320sharon commented 5 months ago

Hi Killian,

I was going through compute_intersection_QC and I noticed you used n_intersect[i] = len(xy_rot[0, :]) on line 344 in SDS_transects.py instead of n_intersect[i] = np.sum(~np.isnan(xy_rot[0, :])). The reason I ask because n_intersect[i] = len(xy_rot[0, :]) includes NaN points in the count while np.sum(~np.isnan(xy_rot[0, :])) does not count the NaN points. In the line xy_rot[0, xy_rot[0,:] < settings['min_chainage']] = np.nan the points that were below the min_chainage are set to NaN, so shouldn't these points not be considered an intersection point?

Here is the relevant portion of the code

            # in case there are no shoreline points close to the transect 
            if len(idx_close) == 0:
                std_intersect[i] = np.nan
                med_intersect[i] = np.nan
                max_intersect[i] = np.nan
                min_intersect[i] = np.nan
                n_intersect[i] = np.nan
            else:
                # change of base to shore-normal coordinate system
                xy_close = np.array([sl[idx_close,0],sl[idx_close,1]]) - np.tile(np.array([[X0],
                                   [Y0]]), (1,len(sl[idx_close])))
                xy_rot = np.matmul(Mrot, xy_close)
                # remove points that are too far landwards relative to the transect origin (i.e., negative chainage)
                xy_rot[0, xy_rot[0,:] < settings['min_chainage']] = np.nan

                # compute std, median, max, min of the intersections
                std_intersect[i] = np.nanstd(xy_rot[0,:])
                med_intersect[i] = np.nanmedian(xy_rot[0,:])
                max_intersect[i] = np.nanmax(xy_rot[0,:])
                min_intersect[i] = np.nanmin(xy_rot[0,:])
                n_intersect[i] = len(xy_rot[0,:])
kvos commented 5 months ago

hi @2320sharon , good that you're having a look at this level as the devil is in the details with these things. I think you are right, good catch. The goal of having multiple points to calculate the intersection is to get a more robust estimate when calculating the median intersection, as when done with only 1 point is not a median intersection anymore. So by default I have set it to 3 points. If 1 of those 3 points was below the min_chainage, then it will still pass the n_intersect condition, although the median was calculated with only 2 points. I am happy with this change, if you have time could you pls PR? thanks

2320sharon commented 5 months ago

Thank you for clearing up my confusion @kvos and for explaining the logic. I'll submit a PR later today with the new code added in. It makes sense that you won't compute a median intersection with less than 3 points.

As always thanks for your all your work maintaining coastsat @kvos !

kvos commented 4 months ago

thanks for this update @2320sharon