Clear Filters
Clear Filters

How to read modern shapefile data

80 views (last 30 days)
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?

Accepted Answer

Jacob Halbrooks
Jacob Halbrooks on 8 Nov 2022
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);
Tianren on 19 Sep 2023
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.

Sign in to comment.

More Answers (1)

Matt Cooper
Matt Cooper on 8 Dec 2022
Edited: Matt Cooper on 9 Dec 2022
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;
[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.')
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,:)]);
S = m_shaperead(strrep(fname,'.shp',''));
ok = true;
catch ME2
rethrow(ME2) % throw the error (revisit)
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;
A = struct2table(S.dbf); % attribute table
Aok = true;
% note: make sure the lat,lon values that go into elements of S are row vectors
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});
% back to geostruct
T = table2struct(movevars(T,'BoundingBox','After','Geometry'));
S = T; % send back the geostruct (overwrite m_shapefile S)
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,
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) );
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);




Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!