bugGNU Octave - Bugs: bug #53422, [octave forge] (mapping) Corrupt...

 
 

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

Submitted by:  None
Submitted on:  Thu 22 Mar 2018 06:00:24 PM UTC  
 
Category:  Octave Forge Package Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  Philip Nienhuis <philipnienhuis>
Originator Name:  Matt Parkan Originator Email:  -email is unavailable-
Open/Closed:  Closed Release:  other
Operating System:  Any

Add a New Comment (Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

( Jump to the original submission)

Thu 05 Apr 2018 01:07:29 PM UTC, comment #15: 

I've just tried reading and writing SHP files with the latest shaperead, shapewrite and dbfwrite functions. Everything seems to work correctly (the geometry and attribute table are displayed correctly in QGIS).

I'll let you know if I find any remaining issues.

Many thanks for the work.

For reference, I've attached the test files in a ZIP and the Octave test code is below:


%% read / write points
s1 = shaperead('point_test_qgis.shp');
shapewrite(s1, 'point_test_octave.shp');
s2 = shaperead('point_test_octave.shp');

%% read / write lines
s1 = shaperead('line_test_qgis.shp');
shapewrite(s1, 'line_test_octave.shp');
s2 = shaperead('line_test_octave.shp');

%% read / write polygons
s1 = shaperead('polygon_test_qgis.shp');
shapewrite(s1, 'polygon_test_octave.shp');
s2 = shaperead('polygon_test_octave.shp');

As a side note, the code I used previously (... F1 = fread(fid1, inf, 'uint8'); ...) which read the SHP and SHX files into an uint8 array was just used to compare the files, not to extract meaningful attribute values (which obviously requires parsing according to the file specification).

(file #43809)

Anonymous
Wed 04 Apr 2018 08:50:21 PM UTC, comment #14: 

I retried on my desktop with QGIS 301, and it also worked there for all geometries.

Closing report as fixed; if you disagree it can always be reopened.
Thanks Matt for your assistance.

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Wed 04 Apr 2018 11:49:08 AM UTC, comment #13: 

With the fixed shapewrite & shaperead (& dbfwrite) I can produce shapefiles for Point, Multipoint, Line and Polygon that display perfectly in ArcMap and QGIS 2.18.
My comment about Line turned out to be due to coordinates being out of the map range - "zoom to layer" helped.

Please confirm if it works for you now; then I can close this bug report.

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Tue 03 Apr 2018 08:50:15 PM UTC, comment #12: 

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 :-)

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Tue 03 Apr 2018 05:51:44 PM UTC, comment #11: 

Thanks for the updates Philip.

I've made a few additional tests in Octave and QGIS for the Point type with the new functions you provided.

I still can't display the geometry in QGIS. However, I think (at least for the Point type) that the problem is coming from the index file (SHX). 

To test this, I prepared a shapefile (point_test_qgis.shp) with QGIS 2.14.4, read it in Octave with shaperead.m and wrote the resulting structure to a new shapefile (point_test_octave.shp) with shapewrite.m.

I then checked if the content of the SHP and SHX files were the same for the Octave and QGIS versions. The SHP files were the same, but the SHX are different. To confirm this, if I substitute the Octave SHX file with the QGIS SHX file, QGIS is able to display the layer correctly.

I've included the test files as an attachement. 


% test SHP equality
fid1 = fopen('point_test_qgis.shp');
F1 = fread(fid1, inf, 'uint8');
fclose(fid1);

fid2 = fopen('point_test_octave.shp');
F2 = fread(fid2, inf, 'uint8');
fclose(fid2);

isequal(F1, F2) % YES

% test SHX equality
fid1 = fopen('point_test_qgis.shx');
I1 = fread(fid1, inf, 'uint8');
fclose(fid1);

fid2 = fopen('point_test_octave.shx');
I2 = fread(fid2, inf, 'uint8');
fclose(fid2);

isequal(I1, I2) % NO

(file #43784, file #43785, file #43786, file #43787)

Anonymous
Tue 03 Apr 2018 12:55:52 PM UTC, comment #10: 

Please try with the attached shapewrite.m and shaperead.m.
Also swap attached dbfwrite.m in place (in the io package).

FYI, ArcMap (ArcGIS) happily swallows and displays all the resulting shapefiles: Point, Multipoint, and Polygon.
Line still seems problematic although initially I / ArcMap had no problems; maybe next Wednesday or Thursday I can look further.

My QGIS desktop 2.18 gave no read errors either but couldn't display them (may be due to some lacking QGIS proficiency between chair and keyboard here) but it does allow editing the attributes.

(file #43780, file #43781, file #43782)

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Sun 01 Apr 2018 10:02:20 PM UTC, comment #9: 

A progress report:

Upon close inspection, shapewrite's code for "Point" geometries; AFAICS conforms exactly to the specs in http://www.esri.com/library/whitepapers/pdfs/shapefile.pdf.
Indeed, I found that the utility program shapechk.exe [*] rather screws up the .shx file.

However, along the way I fixed a bug in dbfwrite.m (io package), it prepended a space to text attributes.
I also adapted shapewrite.m to always write a .dbf file.

As to your remark in comment #4 that shaperead.m produces a different struct than you wrote with shapewrite.m:
Note that the Point geometry has no BoundingBox field - it isn't in the specs, see URL above.

With these changes I'll try again with ArcGIS.

[*] https://allaroundgis.wordpress.com/2013/12/10/shapefile-shp-shape-checker-v3-2/

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Wed 28 Mar 2018 08:54:47 PM UTC, comment #8: 

Í just finished some polygon clipping stuff for the geometry (and thus mapping) package so I'll take up fixing the Point / Multipoint geometries now.
Processing a few more contributions, fixing shaperead's clipping, and then a new mapping package release can be made. Hopefully I have time this coming Easter weekend.

BTW to me ArcView is a bit more authoritative than QGIS -  after all it is ESRI who invented the shapefile format. But don't worry, I'd like to "support" QGIS as well.

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Wed 28 Mar 2018 03:30:41 PM UTC, comment #7: 

FWIW, for the code in the original report (Polygon), ArcMap give me the attached "OP" pic (ArcMap complains about projection but displays the data).

As to comment #5:

With the result of the Multipoint code I see only one point => fails.

Point is failing, too, see attached pic, while the program says it should be legal ArcMap doesn't swallow it.

In case of Line I see two specs that after zooming in on either of them give me 2 complete straight lines whose XY coordinates match both your j=1 and j=2 data.
==> So I think Line works.
(Line pic attached)

For Polygon it seems to work as well.

BTW there is a program called shapechk.exe which sometimes fixes up shapefiles. I used it here; it seems a .dbf file is required in all cases, even if no extra attributes are present; shapewrite doesn't write a .dbf file if no extra attributes are present so that is one thing to improve.


Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Wed 28 Mar 2018 02:36:52 PM UTC, comment #6: 

Thanks very much for the investigation.

I'll pick up later, lots of other things to do.

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Wed 28 Mar 2018 08:09:35 AM UTC, comment #5: 

Part of the code missing for the Point geometry type:


%% Point

s = struct;
j = 1;
s(j).Geometry = 'Point';
s(j).X = 597600;
s(j).Y = 172250;
s(j).BoundingBox = [min(s(j).X), min(s(j).Y); max(s(j).X), max(s(j).Y)];
s(j).Custom1 = 'abc';
s(j).Custom2 = 42;

j = 2;
s(j).Geometry = 'Point';
s(j).X = 598600;
s(j).Y = 170250;
s(j).BoundingBox = [min(s(j).X), min(s(j).Y); max(s(j).X), max(s(j).Y)];
s(j).Custom1 = 'cdf';
s(j).Custom2 = 43;

% test if input and output have the same field values
shapewrite(s, 'point.shp');
s2 = shaperead('point.shp');
test = all(cellfun(@(k) isequal(s.(k), s2.(k)), fieldnames(s2)))

Anonymous
Wed 28 Mar 2018 08:04:55 AM UTC, comment #4: 

I've tried shapewrite (latest snapshot on https://sourceforge.net/p/octave/mapping/ci/default/tree/inst/) with the Point, MultiPoint, Line and Polygon geometry types. I observe the same problem as previously (no geometry appears on screen) when tying to open the files in QGIS 2.14.4.

The range of coordinates was set to fit in the coordinate reference system EPSG:21781 and the projection settings were set accordingly in QGIS, so I don't think this is just a visualization issue.

Moreover, If I try reading the files with shaperead(), I do not get the same structure as the input structure provided to shapewrite().

The code below was used for testing:


%% Multipoint

s = struct;
j = 1;
s(j).Geometry = 'MultiPoint';
s(j).X = [597600 597601 597602 597603 597604];
s(j).Y = [170250 170251 170252 170253 170254];
s(j).BoundingBox = [min(s(j).X), min(s(j).Y); max(s(j).X), max(s(j).Y)];
s(j).Custom1 = 'abc';
s(j).Custom2 = 42;

j = 2;
s(j).Geometry = 'MultiPoint';
s(j).X = [590600 590601 590602 590603 590604 590605];
s(j).Y = [178250 178251 178252 178253 178254 178255];
s(j).BoundingBox = [min(s(j).X), min(s(j).Y); max(s(j).X), max(s(j).Y)];
s(j).Custom1 = 'cdf';
s(j).Custom2 = 43;

% test if input and output have the same field values
shapewrite(s, 'multipoint.shp');
s2 = shaperead('multipoint.shp');
test = all(cellfun(@(k) isequal(s.(k), s2.(k)), fieldnames(s2)))


%% Polygon (no holes)

s = struct;
j = 1;
s(j).Geometry = 'Polygon';
s(j).X = [590610 590620 590620 590610 590610];
s(j).Y = [178220 178220 178210 178210 178220];
s(j).BoundingBox = [min(s(j).X), min(s(j).Y); max(s(j).X), max(s(j).Y)];
s(j).Custom1 = 'abc';
s(j).Custom2 = 42;

j = 2;
s(j).Geometry = 'Polygon';
s(j).X = [597620 597630 597630 597620 597620];
s(j).Y = [170230 170230 170220 170220 170230];
s(j).BoundingBox = [min(s(j).X), min(s(j).Y); max(s(j).X), max(s(j).Y)];
s(j).Custom1 = 'cdf';
s(j).Custom2 = 43;

% test if input and output have the same field values
shapewrite(s, 'polygon.shp');
s2 = shaperead('polygon.shp');
test = all(cellfun(@(k) isequal(s.(k), s2.(k)), fieldnames(s2)))


%% Line

s = struct;
j = 1;
s(j).Geometry = 'Line';
s(j).X = [597600 597601 597602 597603 597604];
s(j).Y = [170250 170251 170252 170253 170254];
s(j).BoundingBox = [min(s(j).X), min(s(j).Y); max(s(j).X), max(s(j).Y)];
s(j).Custom1 = 'abc';
s(j).Custom2 = 42;

j = 2;
s(j).Geometry = 'Line';
s(j).X = [590600 590601 590602 590603 590604 590605];
s(j).Y = [178250 178251 178252 178253 178254 178255];
s(j).BoundingBox = [min(s(j).X), min(s(j).Y); max(s(j).X), max(s(j).Y)];
s(j).Custom1 = 'cdf';
s(j).Custom2 = 43;

% test if input and output have the same field values
shapewrite(s, 'line.shp');
s2 = shaperead('line.shp');
test = all(cellfun(@(k) isequal(s.(k), s2.(k)), fieldnames(s2)))

Anonymous
Fri 23 Mar 2018 04:26:59 PM UTC, comment #3: 

OK I'll try myself.
Busy times so it may take a while.
In the mean time you can try with Point, Multipoint and Line shape.

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Fri 23 Mar 2018 08:12:18 AM UTC, comment #2: 

Thanks for the info Philip.

I've just tested the same minimal example with the new shapewrite function you suggested on Octave Forge. This time the three files are written (.shp, .shx, .dbf) and can be opened with Quantum GIS (2.14.4). However, no polygons are displayed. Perhaps the part of shapewrite which writes the geometry file (.shp) still has a problem with polygons. The attributes table is ok.

The polygon coordinates are listed in clockwise order, they are closed and not self-intersecting, so I don't think there is a problem with the input structure s.

Anonymous
Thu 22 Mar 2018 07:01:18 PM UTC, comment #1: 

Can you please try the shapewrite.m version in the repo?

https://sourceforge.net/p/octave/mapping/ci/default/tree/inst/shapewrite.m

mapping-1.2.1. is from June 2016, since then many fixes to various files, a.o., shapewrite.m, have been pushed.

FYI, as mapping pkg maintainer I run a much newer mapping dev version and that gives me 3 files with your code. I haven't tied reading it in QGIS or so but I think they're OK.

Admittedly a new mapping pkg version is long overdue.
I'm busy fixing shaperead.m to do proper polygon /polyline clipping when reading with a Bounding Box, once that's fixed and a new geometry pkg is released (a mapping pkg dependency, and its maintainer is waiting for a last contribution by me) I'll make a new mapping pkg release.

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Thu 22 Mar 2018 06:00:24 PM UTC, original submission:  

Hi all,

When writing polygons to a shapefile with shapewrite, only a .shp file is produced (no .shx or .dbf files). The .shp file seems to be corrupt and cannot be opened with Quantum GIS (2.14.4) or read with ogrinfo. However, it can be opened with shaperead.

Example:


s = struct;
j = 1;
s(j).Geometry = 'Polygon';
s(j).X = [1 2 2 1 1];
s(j).Y = [2 2 1 1 2];
s(j).BoundingBox = [1 1; 2 2];
s(j).Name = 'att1';

j = 2;
s(j).Geometry = 'Polygon';
s(j).X = [2 3 3 2 2];
s(j).Y = [3 3 2 2 3];
s(j).BoundingBox = [2 3; 2 3];
s(j).Name = 'att2';

shapewrite(s, 'test.shp');
I = shaperead('test.shp');

Anonymous

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #43809:  octave_mapping_test.zip added by None (4KiB - application/x-zip-compressed)
file #43784:  point_test_octave.shp added by None (184B - application/octet-stream)
file #43785:  point_test_octave.shx added by None (124B - application/octet-stream)
file #43786:  point_test_qgis.shp added by None (184B - application/octet-stream)
file #43787:  point_test_qgis.shx added by None (124B - application/octet-stream)
file #43780:  shapewrite.m added by philipnienhuis (23KiB - text/plain)
file #43781:  dbfwrite.m added by philipnienhuis (6KiB - text/plain)
file #43782:  shaperead.m added by philipnienhuis (38KiB - text/plain)
file #43707:  bug53422_OP.PNG added by philipnienhuis (55KiB - image/png)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by jwe (Updated the item)
  • -email is unavailable- added by None (Submitted the item)
  •  

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

     

     

    Follow 23 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-04-05 None Attached File- => Added octave_mapping_test.zip, #43809
    2018-04-04 philipnienhuis StatusReady For Test => Fixed
        Open/ClosedOpen => Closed
    2018-04-03 None Attached File- => Added point_test_octave.shp, #43784
        Attached File- => Added point_test_octave.shx, #43785
        Attached File- => Added point_test_qgis.shp, #43786
        Attached File- => Added point_test_qgis.shx, #43787
    2018-04-03 philipnienhuis Attached File- => Added shapewrite.m, #43780
        Attached File- => Added dbfwrite.m, #43781
        Attached File- => Added shaperead.m, #43782
        StatusIn Progress => Ready For Test
    2018-04-01 philipnienhuis StatusConfirmed => In Progress
    2018-03-28 philipnienhuis StatusNeed Info => Confirmed
    2018-03-28 philipnienhuis Attached File- => Added bug53422_OP.PNG, #43707
        Attached File- => Added bug53422_Line.PNG, #43708
        Attached File- => Added bug53422_Point2.PNG, #43709
    2018-03-23 philipnienhuis Assigned toNone => philipnienhuis
    2018-03-23 jwe Summary[octave-forge] (mapping) Corrupt output when writing polygons to shapefile with shapewrite => [octave forge] (mapping) Corrupt output when writing polygons to shapefile with shapewrite
    2018-03-22 philipnienhuis Summary[octave-forge] (mapping) Corrupt output when writing polygons to shapefile with shapewrite (Mapping package) => [octave-forge] (mapping) Corrupt output when writing polygons to shapefile with shapewrite
    2018-03-22 philipnienhuis SummaryCorrupt output when writing polygons to shapefile with shapewrite (Mapping package) => [octave-forge] (mapping) Corrupt output when writing polygons to shapefile with shapewrite (Mapping package)
    2018-03-22 philipnienhuis StatusNone => Need Info
        Release4.2.1 => other
        Operating SystemMicrosoft Windows => Any

    Back to the top


    Powered by Savane 3.5