Code covered by the BSD License  

Highlights from
KML-Matlab-Shapefile conversion

KML-Matlab-Shapefile conversion

by

 

21 Aug 2009 (Updated )

Matlab script to import KML files into Matlab as well as conversion to Shapefile.

kml_shapefile(filename,outname,varargin)
function [pts_struc] = kml_shapefile(filename,outname,varargin)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matlab script to read KML file into a mapstruct structure array and 
% offers the option of writing to a Shapefile. 
% KML files must be of a uniform type: Point, LineString, LinearRing or Polygon. 
% The output, called 'pts_struc', contains 5-6 fields, Geometry, X, Y, altitude,
% 'Id' and 'name'. The final field is created and populated only if there is 
% exactly one <name> tag for each feature. This mapstruct array can be easily 
% displayed in Matlab using the 'geoshow' command. This output format allows 
% exporting to a Shapefile using the Matlab command, 'shapewrite', which is 
% initiated by entering an output filename into the function. 
%
% INPUTS:
%   map_structure = kml_shapefile('myKML.kml');
%   map_structure = kml_shapefile('myKML.kml','mySHP.shp');
%
% OUTPUTS:
%   The code outputs a mapstruct structure array if only the filename is 
%   provided by the user. The mapstruct is put into memory. A shapefile is 
%   also created if an output shapefile filename is provided. This file is  
%   not brought into memory, however.
%
% by Michael Toomey, University of California Santa Barbara, 
% Department of Geography, Santa Barbara, CA 93106
% mtoomey@geog.ucsb.edu
% last modified: June 1, 2010
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%ARGUMENTS IN%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% determine number of arguments and whether they are the appropriate format
if nargin == 1
    fnm =ischar(filename);
    if fnm == 0
        disp('Filename must be a string')
    end    
elseif nargin == 2
    fnm=ischar(filename);
    onm=ischar(outname);
    if fnm == 0
        disp('Filename must be a string')
    end    
    if onm == 0 
        disp('Output Shapefile name must be a string');
    end
elseif nargin >= 3
    disp('Too many arguments. Arguments should include input KML filename and optionally, output Shapefile name')
end

%%%%%%%%%%%%%%%%%OPEN FILE AND EXTRACT CONTENTS AS CELL ARRAY%%%%%%%%%%%%%%
% file input name and open file, scan text, then close file
fid =fopen(filename);
if fid == -1
    disp('Cannot find file')
end
text = textscan(fid, '%s');
fclose(fid);
% determine size of scanned KML file - each line is one subelement in the
% cell array
sizes = numel(text{1});
body = zeros(sizes,1);

%%%%%%%%%%%%%%%%%%%FINDING TAGS FOR ALL REPORTED FEATURES%%%%%%%%%%%%%%%%%%
% start FOR loop to go through all text statements in the cell array and
% find instances of reported coordinates as well as description of object
% geometry (e.g. point, line, polygon)
for i=1:sizes
    % find indices of cell array, 'text', containing the string,
    % '<coordinates>', assign to 'body1st', and lines containing the string,
    % '</coordinates>', assign to 'body2nd'
    ans0 = strfind(char(text{1}(i)),'<coordinates>');
    body1st(i) = isempty(ans0);
    ans1 = strfind(char(text{1}(i)),'</coordinates>');
    body2nd(i) = isempty(ans1);
    % determine geometry of the objects by looking for one of the three
    % common Geometry tags
    ans2 = strfind(char(text{1}(i)),'<Point>');
    ispoint(i) = isempty(ans2);
    ans3 = strfind(char(text{1}(i)),'<Polygon>');
    ispolygon(i) = isempty(ans3);
    ans4 = strfind(char(text{1}(i)),'<LineString>');
    isline(i) = isempty(ans4);
    ans5 = strfind(char(text{1}(i)),'<LinearRing>');
    islinering(i) = isempty(ans5);
    % here we search for the <Placemark> and <name> tags, if present, to
    % create a name field for the output Shapefile.
    ans6 = strfind(char(text{1}(i)),'<Placemark>');
    isplacemk(i) = isempty(ans6);
    ans7 = strfind(char(text{1}(i)),'<name>');
    isname(i) = isempty(ans7);
    ans8 = strfind(char(text{1}(i)),'</name>');
    isendname(i) = isempty(ans8);
end

% determine which geometry types were reported and assign proper geometry
% to variable, 'geom' which will be later used to assign to the mapstruct
pt_inst = find(ispoint==0);
poly_inst = find(ispolygon==0);
line_inst = find(isline==0);
linering_inst = find(islinering==0);
if numel(pt_inst) >= 1
    geom = 'Point';
elseif numel(poly_inst) >= 1
    geom = 'Polygon';
elseif numel(line_inst) >= 1
    geom = 'Line';    
elseif numel(linering_inst) >= 1
    geom = 'Line';    
end
% determine locations within cell array of reported coordinates
reported_coords1st = find(body1st == 0);
reported_coords2nd = find(body2nd == 0);

% determine locations of 'Placemark' and 'name'
% tags. Instances where name comes right after Placemark, is an indication
% of a reported title/name for each feature. 'intersect' is used to
% determine whether and where the name tag immediately follows the 
% placemark tag, indicating it is an object name and not just a folder name
% or something similar
placemk_inst = find(isplacemk==0);
name_inst = find(isname==0);
endname_inst = find(isendname==0);
[intseca,intsecb,intsecc]=intersect(placemk_inst,name_inst-1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%CREATE MAPSTRUCT ARRAY%%%%%%%%%%%%%%%%%%%%%%%%%
% determine whether there are as many names available as there are features 
% to map out. if affirmative, 'names'=1 and we include 'name' field
% in the mapstruct structure array
if numel(reported_coords1st)==numel(intsecc)
    pts_struc = struct('Geometry',{},'X',{},'Y',{},'altitude',{},'Id',{},'name',{}); 
    names = 1;
else
    pts_struc = struct('Geometry',{},'X',{},'Y',{},'altitude',{},'Id',{}); 
    names = 0;
end

%%%%%%%%%%%%%%%%%EXTRACTING POINTS/NODES FROM KML TAGS%%%%%%%%%%%%%%%%%%%%%
% if all instances of '<coordinate>' and '</coordinate>' are on same line,
% is easier to deal with than if on separate lines. compare
% 'reported_coords1st' and 'reported_coords2nd' to see if equivalent  
equivtest=isequalwithequalnans(reported_coords2nd,reported_coords1st);
% if vectors are the same, assign these lines of text to untrim_coords
% because all information is found on the same line. if not, each point is 
% reported on a separate line and will be dealt with later 
if equivtest == 1
    untrim_coords = text{1}(reported_coords1st);
    % convert to character array
    untrim_coords = char(untrim_coords);
    % remove all instances of 'coordinates' from each string and then convert
    % to number for insertion into double array. 
    for i=1:numel(reported_coords1st)
        temp1 = strrep(untrim_coords(i,:),'<coordinates>','');
        temp2 = strrep(temp1,'</coordinates>','');
        % the first token extracted is the longitude - write into the
        % structure array
        [temp3,temp4] = strtok(temp2,','); 
        pts_struc(i).X = str2num(temp3);
        % the second token extracted is the latitude
        [temp5,temp6] = strtok(temp4,','); 
        pts_struc(i).Y = str2num(temp5);
        % the third token extracted is the altitude
        [temp7,temp8] = strtok(temp6,','); 
        pts_struc(i).altitude = str2num(temp7);
        % also insert the correct geometry into the structure
        pts_struc(i).Geometry = geom;
        % insert appropriate Id number
        pts_struc(i).Id = i;
        % if names=1 (i.e. there is a name tag for each object, then strip 
        % tags and insert appropriate text into 'name' field of pts_struc,
        % with some fancy character array flippage
        if names ==1
            nametxt=text{1}(name_inst(intsecc(i)):endname_inst(intsecc(i)));
            test=char(nametxt);ans=test';ans=rot90(ans(:));
            % strip tags
            ans=strrep(ans,'<name>','');
            pts_struc(i).name = strrep(ans,'</name>','');
        end
    end
% here we deal with the case if coordinates, with their tags are not all
% reported on the same line
elseif equivtest == 0
    % determine how many points in each coordinates cluster
    numpts = (reported_coords2nd-reported_coords1st)-1;
    for i=1:numel(numpts)
        temp1=zeros(numpts(i),3);
        temp2=char(text{1}(reported_coords1st(i)+1:reported_coords2nd(i)-1));
        for j=1:numpts(i)
            % the first token extracted is the longitude 
            [temp3,temp4] = strtok(temp2(j,:),','); 
            temp1(j,1) = str2num(temp3);
            % the second token extracted is the latitude
            [temp5,temp6] = strtok(temp4,','); 
            temp1(j,2)= str2num(temp5);
            % the third token extracted is the altitude
            [temp7,temp8] = strtok(temp6,','); 
            temp1(j,3) = str2num(temp7);
        end
        % load node coordinates for each object, stored in 'temp1' into the
        % corresponding place in the mapstruct
        pts_struc(i).X = temp1(:,1);
        pts_struc(i).Y = temp1(:,2);
        % you can uncomment this next line if you want to have an altitude
        % assigned to each node in a multipoint line or polygon feature,
        % but that does not aid the mapshow presentation and crashes the
        % shapewrite.m program
        %pts_struc(i).altitude = temp1(:,3);
        pts_struc(i).altitude = temp1(1,3);
        % also insert the correct geometry into the structure
        pts_struc(i).Geometry = geom;
        % insert appropriate Id number
        pts_struc(i).Id = i;
        % if names=1 (i.e. there is a name tag for each object, then strip 
        % tags and insert appropriate text into 'name' field of pts_struc,
        % with some fancy character array flippage
        if names ==1
            nametxt=text{1}(name_inst(intsecc(i)):endname_inst(intsecc(i)));
            test=char(nametxt);ans=test';ans=rot90(ans(:));
            % strip tags
            ans=strrep(ans,'<name>','');
            pts_struc(i).name = strrep(ans,'</name>','');
        end
        
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%WRITING SHAPEFILE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% write Shapefile to disk if output name has been provided
if nargin ==2
    shapewrite(pts_struc,outname);
    comment = cat(2,'Shapefile ',outname,' written to disk');
    disp(comment)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%END OF CODE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Contact us