ncssar / sartopo_python

Python calls for the caltopo / sartopo API
GNU General Public License v3.0
14 stars 2 forks source link

crop: incorrect results #41

Closed caver456 closed 2 years ago

caver456 commented 2 years ago

have only observed this a few times. Both times it involved multiple enter-and-exits, resulting in multiple crop result lines:

image

see exported json (of the pre-cropped map, on the left) incorrect_crop_20220218.json

caver456 commented 2 years ago

this is a bug in the home-brew 'intersection2' function which is used by crop() instead of shapely.intersection, since shapely.intersection does not preserve complex lines (i.e. with internal crossovers).

The 'known issue' at the bottom of these comments describes the exact problem. So, now is the time to fix it:

# intersection2(targetGeom,boundaryGeom)
# we want a function that can take the place of shapely.ops.intersection
#  when the target is a LineString and the boundary is a Polygon,
#  which will preserve complex (is_simple=False) lines i.e. with internal crossovers

# walk thru the points (except fot the last point) in the target shape(line):
#    A = current point, B = next point
#  A and B both inside boundary?  --> append A to output coord list
#  A inside, B outside --> append A; append point at instersection of this segment with boundary; don't append B
#  A outside, B inside --> append B; append point at intersection of this segment with boundary; don't append A
#  A outside, B outside --> don't append either

# known issue: straight segments that are 'clipped' by the boundary corner,
#  i.e. A and B are both outside, but a portion of the AB segment is inside,
#  will be omitted from the result, since only the drawn vertices are checked.
caver456 commented 2 years ago

Added code to handle the fourth case (A outside, B outside):

            else: # neither endpoint is inside the boundary
                abl=LineString([ap,bp])
                mp=abl.intersection(boundaryGeom.exterior)
                # the result will be a single disjoint line segment inside the boundary,
                #  with both vertices touching the boundary;
                # the result is a multipoint, which has no .coords attribute
                #  see https://stackoverflow.com/a/51060918
                mpcoords=[(p.x,p.y) for p in mp]
                nextInsidePointStartsNewLine=True
                outLines.append(mpcoords)

image

Result looks good. Fixed - still in the plans_console repo during development.