Fix Segmentation for VBET NZ

[philipbaileynar]

VBET is failing to segment level path polygons. The end result, projects like this, simply produce VBET whole level path polygons. This started happening after Jordan introduced his new segmentation code. And I think it has something to do with transforming geometries and the linear EPSG that we use in New Zealand.

The problematic code is in this file, I think.

Can you try running/debugging/fixing VBET in New Zealand so that segmentation works and we get DGOs and IGOs please. YOu will know when it's working when you get DGOs that are 100-500 metres in length with corresponding IGOs. Currently the IGOs output is empty.

I strongly recommend using the new mini TEST project for New Zealand.

https://data.riverscapes.net/p/4d7ae7bf-9bf9-4eeb-bda0-6a8b0961ffa7/

[MattReimer]

Took me a while to get an environment going and I chased my tail a bit but I rebased the nz_lidar_dems branch on top of master and finally got an output for the mini TEST and it seems segmented to me:

I'm going to try queuing up a private run of the goofy project you linked above for huc 586561 and see if we're still getting a weird result.

UPDATE: I can't get VBET to even finish on cybercastor. Seems like it's the pitfill problem again which probably means raster crop issues. I'll investigate this tomorrow

pitfill.tif --dinfflowdir_ang $TAUDEM_DIR/intermediates/dinfflowdir_ang.tif --dinfflowdir_slp $TAUDEM_DIR/outputs/dinfflowdir_slp.tif --reach_codes 33400,46000,46003,46006,46007,55800 --meta "Runner=Cybercastor" --verbose --temp_folder $VBET_TEMP
0...10...20...30...40...50...60...70...80...90...<<VBET PROCESS ENDED WITH AN ERROR>>

CYBERCASTOR-76b7c882-f489-45b0-a916-4d4f9cbfd252_36748700-ecfd-4b8b-b204-327c0a94e5cb.log

[MattReimer]

I thought it might be a memory or resources issue so I tried it again with a bigger box but it still failed.

[info] [Clean Raster Regions] Generate regions
--
2
PROGRESS: 0 / 29,318  (0%)  Applying nodata mask     (Ellapsed: 0.0 seconds)
PROGRESS: 21,417 / 29,318  (73%)  Applying nodata mask     (Ellapsed: 20.0 seconds)
[info] [Applying nodata mask] Completed 29,318 operations.  Total Time: 27.0 seconds
[info] [threshold] Thresholding at 0.95
PROGRESS: 0 / 29,318  (0%)  Thresholding at 0.95     (Ellapsed: 0.0 seconds)
100 - done.
0...10...20...30...40...50...60...70...80...90...<<VBET PROCESS ENDED WITH AN ERROR>>

This really feels like a segfault error.

UPDATE: I increased the resources of the box to:

  "resources": {
    "cpu": 16384,
    "memory": 32768,
    "disk": 120
  }

and it seems to have pushed through this section. Fingers crossed....

UPDATE2: Ok, this is simply brutal

https://cybercastor.riverscapes.net/893deb5d-bcdf-4b83-9a50-0e5e00b4a410/4cf7355c-fbe1-4f60-b21d-59218a0745fa

The last time this ran it took a day and a half so I don't know what to do except wait. I can't see anything wrong with the sample HUC so we'll see if this job finishes.

[MattReimer]

Ok. this may be easier than I thought. My first clue was that the igo layer is not in the right place. Look at where it is next to the DGO and raster layers:

Then if you look at them side by side it's clearly transposing some coordinates:

Ok, new theory:

You're passing output_epsg into generate_igo_points() but then never using it. That seems weird. Probably a red herring though

The main difference I can see between all the north american stuff and the NZ stuff is that the NA stuff uses the "if" (below) and the NZ uses the "else"

        # if the line layer is not projected, we need to project it to a UTM zone
        # in order to get accurate lengths
        if line_lyr.spatial_ref.IsProjected() == 0:
            extent_poly = get_rectangle_as_geom(line_lyr.ogr_layer.GetExtent())
            extent_centroid = extent_poly.Centroid()
            utm_epsg = get_utm_zone_epsg(extent_centroid.GetX())
            transform_ref, transform = VectorBase.get_transform_from_epsg(
                line_lyr.spatial_ref, utm_epsg)
            transform_back = osr.CoordinateTransformation(
                transform_ref, line_lyr.spatial_ref)
        else:
            # if the line layer is projected, we can use the spatial reference directly
            transform_ref = line_lyr.spatial_ref
            transform = osr.CoordinateTransformation(line_lyr.spatial_ref, transform_ref)
            transform_back = transform

get_utm_zone_epsg includes a call to:

        # https://github.com/OSGeo/gdal/issues/1546
        # Set the axis mapping to be the same as the input. This might prove problematic in cases where the input is
        out_spatial_ref.SetAxisMappingStrategy(in_srs.GetAxisMappingStrategy())

but that never gets called in the "else" case.

which causes about 95% of all problems we've seen where coordinates are getting swapped. It could be that some dataset in the NZ data (but not the small sample one which is why we're not seeing the problem there) has its axis mapping strategy set improperly so that may be it.

[MattReimer]

For reference here is the small 0000001 I ran. For whatever reason this segmented just fine. As of now I can't explain why this is.

https://staging.data.riverscapes.net/rv/c6158dbb-4734-4114-b07f-a5001e388819?do=1&lo=1&vm=%7E1%7E8%7E9%7E11%7E33%7E45&bl=SAT&vp=172.122747%7E-43.112569%7E14.592

[MattReimer]

Final update from me: The cybercastor job finished

https://data.riverscapes.net/rv/f331aca6-60e4-446b-b74c-66ab0657d680?do=1&lo=0&vm=%7E23%7E17&bl=TOPO&vp=168.325854%7E-45.756798%7E12.575

Seems like it did segment but the rest of the rivershed is a bit of a mess so I don't really know. I've taken this as far as I can without any more input so I'm now passing this off.

FYI @philipbaileynar

[weblorin]

@jtgilbert Are you able to look at this?

[jtgilbert]

sorry, just getting back into the office, and I feel bad for @MattReimer chasing his tail so much when I'd been through a bunch of this haha.