How to read modern shapefile data

70 次查看（过去 30 天）
Chad Greene
Chad Greene 2022-11-8
评论： Tianren 2023-9-19
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?

请先登录，再回答此问题。

采纳的回答

Jacob Halbrooks
Jacob Halbrooks 2022-11-8
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);
  2 个评论
Tianren
Tianren 2023-9-19
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.

请先登录，再进行评论。

更多回答（1 个）

Matt Cooper
Matt Cooper 2022-12-8
编辑：Matt Cooper 2022-12-9
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

类别

在 Help Center 和 File Exchange 中查找有关 Other Formats 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by