How to read modern shapefile data

It seems increasingly common these days that data providers are saving shapefile data as PolygonZ or PolyLineZ data, rather than the old elevationless Polygn or PolyLine format. This is a big problem, because Matlab cannot read PolygonZ or PolyLineZ format:

S = shaperead'myfile.shp');
Error using openShapeFiles>readHeaderTypeCode
Unsupported shape type PolyLineZ (type code = 13). Use readgeotable instead.
Error in openShapeFiles (line 24)
headerTypeCode = readHeaderTypeCode(shpFileId,callingFcn);
Error in shaperead (line 210)
= openShapeFiles(filename,'shaperead'); 

The error message suggests using readgeotable, but that's no good because there's no way to get the actual coordinate data from a geopolyshape. Until now, I've been dealing with this by opening shapefiles in QGIS, saving them without Z data, and then reading the data into Matlab via shaperead. However, today I need to open a dataset containing many hundreds of shapefiles in PolyLineZ, and manual conversion is not feasible with so many files.

Is there any way to read PolygonZ or PolyLineZ data into Matlab?

Answer 1

[Jacob Halbrooks]

You can access latitude-longitude coordinates from a geopolyshape using the geotable2table function. We recognize that there is a need to make coordinate data access easier, but in the meantime you can use this approach as a workaround. For example, assuming your table contains geopolyshape data:

GT = readgeotable("myfile.shp");
T = geotable2table(GT,["Lat","Lon"]); 

You can also extract the data into NaN-delimited arrays using polyjoin:

[lat,lon] = polyjoin(T.Lat,T.Lon);

[Chad Greene]

That's helpful; thanks @Jacob Halbrooks!

[Tianren]

This is still inconvenient for PolyLineM type.

The part information will be lost.

Why matlab can't even support what has already be realized long ago in PyShp?

PyShp can easily deal with all kinds of shapefiles.

Answer 2

[Matt Cooper]

For users without readgeotable and/or geotable2table (pre-2021b), here is a function that will try shaperead and if a PolyLineZ error is thrown, it will use the m_map toolbox function m_shaperead instead, which reads -z type shapefiles. It converts the output of m_shaperead into a geostruct that mimics the output of shaperead. The second input namedargs is a varargin-like name-value cell-array because the function copied below is nested inside another function that accepts all name-value pairs accepted by shaperead, so just delete it or replace it with 'varargin' if you don't want that. It only supports PolyLineZ right now but it should be clear how to add support for additional -z type shapefiles.

function [S,A] = tryshaperead(fname,namedargs);

ok = false;
try
    [S,A] = shaperead(fname,namedargs{:});
catch ME
    
    if strcmp(ME.message,'Unsupported shape type PolyLineZ (type code = 13).')
        
        % try m_map/m_shaperead. note: m_shaperead argument UBR (User Bounding
        % Rectangle) assumes format [minX  minY  maxX maxY] whereas shaperead
        % BoundingBox is [minX  minY ; maxX maxY]. do the conversion if needed.
        warning('Requested shapefile is type PolyLineZ. Using m_shaperead.')
        try
            if ismember({'BoundingBox'},namedargs(1:2:end)) % elements 1:2:end are parameters
                B = namedargs{find(ismember('BoundingBox',namedargs(1:2:end)))+1};
                S = m_shaperead(strrep(fname,'.shp',''),[B(1,:),B(2,:)]);
            else
                S = m_shaperead(strrep(fname,'.shp',''));
            end
            ok = true;
        catch ME2
            rethrow(ME2) % throw the error (revisit)
        end
    end
end

if ok
    % the Data produced by m_shaperead may need to be wrangled to a ueseable
    % format. return to this later. for now, convert the attributes to A and the
    % lat,lon to S
    Aok = false;
    try
        A = struct2table(S.dbf); % attribute table
        Aok = true;
    catch
    end
    
    % note: make sure the lat,lon values that go into elements of S are row vectors
    try
        lat = cellfun(@transpose,cellfun(@(x)vertcat(x(:,2)),S.ncst,'Uni',0),'Uni',0);
        lon = cellfun(@transpose,cellfun(@(x)vertcat(x(:,1)),S.ncst,'Uni',0),'Uni',0);
        
        % above returns latlon in cell array format, suitable for geostruct-type
        % objects and plotting with geoshow. below converts to nan-separated
        % lists, suitable for axesm-based functions like plotm.
        
        % activate this and add it to output for axesm-based compatibility.
        % [lat,lon] = polyjoin(lat,lon);
        
        switch S.ctype
            case 'polylineZ'
                T = geostructinit('Line',numel(lat),'fieldnames',S.fieldnames);
                [T(1:numel(lon)).Lon] = lon{:};
                [T(1:numel(lat)).Lat] = lat{:};
                T = updategeostruct(T); % get the bounding box of each element
                T = struct2table(T); % for movevars
                
                if Aok % we have an attribute table, join it with S
                    fields = S.fieldnames;
                    for n = 1:numel(fields)
                        T.(fields{n}) = A.(fields{n});
                    end
                end
                
                % back to geostruct
                T = table2struct(movevars(T,'BoundingBox','After','Geometry'));
        end
        S = T; % send back the geostruct (overwrite m_shapefile S)
    catch
    end
end


function S = geostructinit(geometry,numfeatures,varargin);
%geostructinit initializes a geostructure S of size 1xnumfeatures
% 
%  Syntax
% 
%  S = geostructinit(geometry,numfeatures); 
%  Creates a geostruct S of size(1,numfeatures) and specified geometry
%  'line','point', or 'polygon'.
% 
%  S = geostructinit(geometry,numfeatures,'fieldnames',fieldnames);
%  Creates a geostruct S of size(1,numfeatures) and specified geometry
%  'line','point', or 'polygon' and fields specified by 'fieldnames'
%  name-value input.
% 
%  Author: Matt Cooper, Sep-23-2022, https://github.com/mgcooper
%--------------------------------------------------------------------------
p                 = inputParser;
p.FunctionName    = 'geostructinit';
p.addRequired(   'geometry',              @(x)ischar(x)           );
p.addRequired(   'numfeatures',           @(x)isnumeric(x)        );
p.addParameter(  'fieldnames',   '',      @(x)ischar(x)|iscell(x) );

p.parse(geometry,numfeatures,varargin{:});

geometry    = p.Results.geometry;
numfeatures = p.Results.numfeatures;
fieldnames  = p.Results.fieldnames;
%--------------------------------------------------------------------------

% make sure the geometry is capitalized
geometry(1) = upper(geometry(1));

% init a geo struct
[S(1:numfeatures).Geometry]  = deal(geometry);
[S(1:numfeatures).Lon]       = deal(nan);
[S(1:numfeatures).Lat]       = deal(nan);

for n = 1:numel(fieldnames)
   [S(1:numfeatures).(fieldnames{n})] = deal(nan);
end