jblindsay / whitebox-tools

An advanced geospatial data analysis platform
https://www.whiteboxgeo.com/
MIT License
967 stars 161 forks source link

FillBurn using streams with no fields #434

Open Aethelred-Von-Stick opened 2 weeks ago

Aethelred-Von-Stick commented 2 weeks ago

Hi, I've been using FillBurn with a streams shapefile without any fields. You can create a shapefile with no fields by creating a dataset with a temporary field, and then deleting it. I've mocked up a script that will do this. I don't know if it's a great idea to be creating these kinds of shapefiles, but the error message that was produced wasn't helpful in figuring out what was wrong.

OS: Ubuntu 22.04 WhiteboxTools Open Core: v2.4.0 GDAL: v3.8.2

Create one vector with fields and one without using the script below.

./create_vectors.py

Create a DEM.

gdal_create -outsize 10 10 -burn 1 -a_srs EPSG:4326 -a_ullr -1 1 1 0 dem.tif

Running FillBurn with the streams with fields works fine.

./whitebox_tools -r=FillBurn --dem=dem.tif --streams=with_fields.shp --output=burned.tif
****************************
* Welcome to FillBurn      *
* Powered by WhiteboxTools *
* www.whiteboxgeo.com      *
****************************
Reading streams data...
Reading DEM data...
Rasterizing Streams: 18446744073709551615%
Line Thinning 1: 0%
...
Line Thinning 1: 100%
Initializing output: 0%
...
Initializing output: 100%
Progress: 2%
...
Progress: 36%
Filling DEM: 37%
...
Filling DEM: 137%
Saving data...
Output file written
Elapsed Time (excluding I/O): 0.0s

Running FillBurn on streams without fields produces this unhelpful back trace.

RUST_BACKTRACE=full whitebox_tools -r=FillBurn --dem=test.tif --streams=no_fields.shp --output=burned_2.tif
****************************
* Welcome to FillBurn      *
* Powered by WhiteboxTools *
* www.whiteboxgeo.com      *
****************************
Reading streams data...
thread 'main' panicked at /home/runner/work/whitebox-tools/whitebox-tools/whitebox-common/src/utils/byte_order_reader.rs:73:44:
called `Result::unwrap()` on an `Err` value: Error { kind: UnexpectedEof, message: "failed to fill whole buffer" }
stack backtrace:
   0:     0x7f33d4311332 - <unknown>
   1:     0x7f33d43581fc - <unknown>
   2:     0x7f33d430e60f - <unknown>
   3:     0x7f33d4311104 - <unknown>
   4:     0x7f33d431271b - <unknown>
   5:     0x7f33d4312473 - <unknown>
   6:     0x7f33d4312bbd - <unknown>
   7:     0x7f33d4312a92 - <unknown>
   8:     0x7f33d4311806 - <unknown>
   9:     0x7f33d43127c4 - <unknown>
  10:     0x7f33d356e465 - <unknown>
  11:     0x7f33d356e953 - <unknown>
  12:     0x7f33d429590d - <unknown>
  13:     0x7f33d428a036 - <unknown>
  14:     0x7f33d4286c4a - <unknown>
  15:     0x7f33d387f60a - <unknown>
  16:     0x7f33d3d9f17e - <unknown>
  17:     0x7f33d3f3606d - <unknown>
  18:     0x7f33d3f331fa - <unknown>
  19:     0x7f33d35c2313 - <unknown>
  20:     0x7f33d3eda1b9 - <unknown>
  21:     0x7f33d4309155 - <unknown>
  22:     0x7f33d3f370b5 - <unknown>

create_vectors.py

#! /usr/bin/env python3

"""create_vectors.py"""

from osgeo import gdal, ogr, osr

srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
drv = gdal.GetDriverByName("ESRI Shapefile")
geom = ogr.CreateGeometryFromWkt("LINESTRING (0 0,0 1)")

# Create a dataset without any fields.
ds_1 = drv.Create("no_fields.shp", 0, 0, 0, gdal.GDT_Unknown)
lyr_1 = ds_1.CreateLayer("no_fields", srs, ogr.wkbLineString)
field_defn = ogr.FieldDefn("tmp_field", ogr.OFTInteger)

# If no field is created, an FID field will be created instead.
lyr_1.CreateField(field_defn)
feat_1 = ogr.Feature(lyr_1.GetLayerDefn())

feat_1.SetGeometry(geom)
feat_1.SetField("tmp_field", 1)
lyr_1.CreateFeature(feat_1)

# Delete the only field.
lyr_1.DeleteField(0)

# Create a dataset with a field.
ds_2 = drv.Create("with_fields.shp", 0, 0, 0, gdal.GDT_Unknown)
lyr_2 = ds_2.CreateLayer("with_fields", srs, ogr.wkbLineString)

field_defn = ogr.FieldDefn("ID", ogr.OFTInteger)

lyr_2.CreateField(field_defn)
feat_2 = ogr.Feature(lyr_2.GetLayerDefn())
feat_2.SetField("ID", 1)
feat_2.SetGeometry(geom)
lyr_2.CreateFeature(feat_2)