[Octave-bug-tracker] [bug #53422] [octave forge] (mapping) Corrupt output when writing polygons to shapefile with shapewrite

[Philip Nienhuis]

Follow-up Comment #12, bug #53422 (project octave):

Ah I missed your new comment.

In the mean time I've pushed various shaperead and shapewrite fixes here
(essentially the same as the new file versions in my previous comment):
http://hg.code.sf.net/p/octave/mapping/rev/53cfa74ff769

and dbfwrite fixes here:
http://hg.code.sf.net/p/octave/io/rev/5308e5e42be1

I couldn't convince QGIS 3.0.1 to display the data, ILWIS 4alpha did a little
better. (That said, I find QGIS 2.18's ease-of-use a lot better than QGIS
3.0.1's.)
ArcMap did the best: only Line geometries still won't display, the other
geometries display fine and attributes of all geometries incl. Line are
editable.

Indeed in the .shx file the header wasn't updated. Otherwise AFAICS it just
conforms to the formal ESRI specs. But later on I'll check with the files you
uploaded.
For Point types the records lengths were computed incorrectly, for Multipoint
the XY coordinates were in the wrong order.
As to (poly)Line and Polygon geometries, within shapefiles these are exactly
identical save the geometry type bytes so if Polygons are properly written and
shown in ArcMap, Lines should as well.

Hopefully later this week I can have another look.

Anyway thank you for your assistance in fixing this; I've had very few chances
to test the shapewrite function with independent use case and user :-)


BTW if you want to check the low-level formats of .shp and .shx in Octave
yourself, don't use int8 but proceed along these lines (from the top of my
head):

fp = fopen ("line.shp");
## Read header
fread (fp, 7, "*int32", 0, "ieee-be") ## "*int32" maps double to int32
fread (fp, 2, "*int32", 0, "ieee-le") ## file version and -length in words (2
bytes)
fread (fp, 8, "double", 0, "ieee-le") ## XYMZ overall Bounding box
            ## <now we're at end of file header at 100 bytes from start>
## Read first feature record descriptor
fread (fp, 2, "*int32", 0, "ieee-be") ## feature number and -length in words
## Read first record byte
fread (fp, 1, "*int32", 0, "ieee-le") ## feature type
fread (fp, 4, "double", 0, "ieee-le") ## XY bounding box of this feature
fread (fp, 3, "*int32", 0, "ieee-le") ## - number of Line parts (here: 1),
                                      ## - total nr of points (here: 2),
                                      ## - relative pointer to first part
fread (fp, 4, "double", 0, "ieee-le") ## Read (in this case:) 2 XY pairs
ftell (fp)                            ## File position in bytes;
                                      ## subtract 100 and divide by 2 to
                                      ## get record offset in .shx of this
                                      ## feature

## Read second feature record descriptor
fread (fp, 2, "*int32", 0, "ieee-be") 
fread (fp, 1, "*int32", 0, "ieee-le") ## feature type
fread (fp, 4, "double", 0, "ieee-le") ## XY bounding box of 2nd feature
:
etc
:
fclose (fp)


Similar for the shx file.
See ESRI's shapefile doc (URL in an earlier comment here.)
If you err with nr of numbers or types to read, of big- or low-endian, just
do:
fseek (fp, 100, 'bof') 
to return to the end of the file header to start over :-)


    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?53422>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/