1## Copyright (C) 2018-2020 Philip Nienhuis 2## 3## This program is free software: you can redistribute it and/or modify it 4## under the terms of the GNU General Public License as published by 5## the Free Software Foundation, either version 3 of the License, or 6## (at your option) any later version. 7## 8## This program is distributed in the hope that it will be useful, but 9## WITHOUT ANY WARRANTY; without even the implied warranty of 10## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11## GNU General Public License for more details. 12## 13## You should have received a copy of the GNU General Public License 14## along with this program. If not, see 15## <https://www.gnu.org/licenses/>. 16 17## -*- texinfo -*- 18## @deftypefn {} {@var{kml} =} kmlread (@var{fname}) 19## (EXPERIMENTAL) Read a Google kml file and return its contents in a struct. 20## 21## @var{name} is the name of a .kml file. 22## Currently kmlread can read Point, LineString, Polygon, Track and Multitrack 23## entries. 24## 25## @var{kml} is a struct with fields Type, Lat, Lon, Ele, Time and Name. If 26## the .kml file contained LineStrings and/or Polygons also field BoundingBox 27## is present containing suitable values for the relevant items. 28## 29## @seealso{kmzread} 30## @end deftypefn 31 32## Author: Philip Nienhuis <prnienhuis@users.sf.net> 33## Created: 2018-05-03 34 35function outp = kmlread (fname) 36 37 fid = fopen (fname, "r"); 38 if (fid > 0) 39 xml = fread (fid, Inf, "*char")'; 40 fclose (fid); 41 else 42 error ("kmlread: couldn't read %s", fname); 43 endif 44 45 outp = repmat (struct ("Type", "-", ... 46 "Lat", [], ... 47 "Lon", [], ... 48 "Ele", [], ... 49 "Time", [], ... 50 "Name", "-"), 0, 1); 51 52 ## Points 53 is = 1; 54 ipt = 0; 55 while (is > 0) 56 [pnt, il, is] = getxmlnode (xml, "Point", is, 1); 57 if (is > 0) 58 outp(end+1).Type = "Point"; 59 ++ipt; 60 pnam = regexp (xml(max (1, il - 1000):il+7), ... 61 '<Placemark.*?>.*?<name>(.+?)</name>.*?<Point', "tokens"); 62 if (! isempty (pnam)) 63 pname = cell2mat (pnam){end}; 64 outp(end).Name = pnam; 65 else 66 outp(end).Name = ["Point #" num2str(ipt)]; 67 endif 68 coords = strrep (getxmlnode (pnt, "coordinates", 1, 1), ",", " "); 69 xyz = sscanf (coords, "%f", Inf); 70 outp(end).Lon = xyz(1); 71 outp(end).Lat = xyz(2); 72 if (numel (xyz) == 3) 73 outp(end).Ele = xyz(3); 74 endif 75 endif 76 endwhile 77 78 ## Linestrings / Polylines 79 is = 1; 80 iln = 0; 81 while (is > 0) 82 [lin, il, is] = getxmlnode (xml, "LineString", is, 1); 83 if (is > 0) 84 outp(end+1).Type = "Line"; 85 ++iln; 86 lnam = regexp (xml(max (1, il - 1000):il+12), ... 87 '<Placemark.*?>.*?<name>(.+?)</name>.*?<LineString', "tokens"); 88 if (! isempty (lnam)) 89 lnam = cell2mat (lnam){end}; 90 outp(end).Name = lnam; 91 else 92 outp(end).Name = ["Line #" num2str(iln)]; 93 endif 94 coords = getxmlnode (lin, "coordinates", 1, 1); 95 ## Check if we have Ele (Z). Coordinate tuples are separated by spaces 96 lines = strsplit (strtrim (coords), {" ", "\n"}); 97 nd = numel (strfind (lines{1}, ",")) + 1; 98 xyz = reshape (sscanf (strrep (coords, ",", " "), "%f", Inf)', nd, [])'; 99 outp(end).Lon = xyz(:, 1); 100 outp(end).Lat = xyz(:, 2); 101 outp(end).BoundingBox = [min(outp(end).Lon), max(outp(end).Lon); ... 102 min(outp(end).Lat), max(outp(end).Lat)]; 103 if (nd > 2) 104 outp(end).Ele = xyz(:, 3); 105 outp(end).BoundingBox = [outp(end).BoundingBox; 106 min(outp(end).Ele), max(outp(end).Ele)]; 107 else 108 outp(end).BoundingBox = [outp(end).BoundingBox; 109 NaN, NaN]; 110 endif 111 endif 112 endwhile 113 114 ## Polygons 115 is = 1; 116 ip = 0; 117 while (is > 0) 118 [pol, il, is] = getxmlnode (xml, "Polygon", is, 1); 119 if (is > 0) 120 outp(end+1).Type = "Polygon"; 121 ++ip; 122 pnam = regexp (xml(max (1, il - 1000):il+12), ... 123 '<Placemark.*?>.*?<name>(.+?)</name>.*?<Polygon', "tokens"); 124 if (! isempty (pnam)) 125 pnam = cell2mat (pnam){end}; 126 outp(end).Name = pnam; 127 else 128 outp(end).Name = ["Polygon #" num2str (ip)]; 129 endif 130 ir = 1; 131 ## First get outer ring 132 [ring, ~, ir] = getxmlnode (pol, "outerBoundaryIs", ir); 133 coords = strrep (getxmlnode (ring, "coordinates", 1, 1), ",", " "); 134 xyz = reshape (sscanf (coords, "%f", Inf)', 3, [])'; 135 ## FIXME check on CCW 136 outp(end).Lon = xyz(:, 1); 137 outp(end).Lat = xyz(:, 2); 138 outp(end).Ele = xyz(:, 3); 139 ## Next, any inner rings 140 while (ir > 0) 141 [ring, ~, ir] = getxmlnode (pol, "innerBoundaryIs", ir); 142 if (ir > 0) 143 coords = strrep (getxmlnode (ring, "coordinates", 1, 1), ",", " "); 144 xyz = reshape (sscanf (coords, "%f", Inf)', 3, [])'; 145 ## FIXME check on CW 146 outp(end).Lon = [outp(end).Lon; NaN; xyz(:, 1)]; 147 outp(end).Lat = [outp(end).Lat; NaN; xyz(:, 2)]; 148 outp(end).Ele = [outp(end).Ele; NaN; xyz(:, 3)]; 149 endif 150 endwhile 151 outp(end).BoundingBox = [min(outp(end).Lon), max(outp(end).Lon); ... 152 min(outp(end).Lat), max(outp(end).Lat); ... 153 min(outp(end).Ele), max(outp(end).Ele)]; 154 endif 155 endwhile 156 157 ## Tracks and MultiTracks 158 ptrnT = '<when>(.*?)</when>'; 159 is = em = 1; 160 it = 0; 161 while (is > 0) 162 if (em <= is) 163 ## Try to get extent of multiTrack node 164 [mtrk, ~, em] = getxmlnode (xml, "gx:MultiTrack", is, 1); 165 if (em > 0) 166 ## Get specific attributes for this Multitrack 167 mtid = getxmlattv (mtrk, "id"); 168 intpm = str2double (getxmlnode (mtrk, "gx:interpolate", 10, 1)); 169 ## MultiTrack; start new struct item for next tracks 170 outp(end+1).Type = "Track"; 171 ++it; 172 if (isempty (mtid)) 173 outp(end).Name = ["Track #" num2str(it)]; 174 else 175 outp(end).Name = mtid; 176 endif 177 endif 178 endif 179 180 [trk, ~, is] = getxmlnode (xml, "gx:Track", is, 1); 181 if (is > 0) 182 if (! em) 183 ## Separate track; start new struct item 184 outp(end+1).Type = "Track"; 185 ++it; 186 tnam = getxmlattv (trk, "id"); 187 if (isempty (tnam)) 188 outp(end).Name = ["Track #" num2str(it)]; 189 else 190 outp(end) = tnam; 191 endif 192 outp(end).BoundingBox = [min(outp(end).Lon), max(outp(end).Lon); ... 193 min(outp(end).Lat), max(outp(end).Lat); ... 194 min(outp(end).Ele), max(outp(end).Ele)]; 195 endif 196 times = cell2mat (regexp (trk, ptrnT, "tokens")); 197 ltime = cellfun ("numel", times); 198 if (any (ltime < 20)) 199 ## Year resolution 200 times(ltime == 4) = ... 201 cellfun (@(x) [x "-00-00T00:00:00Z"], times(ltime == 4), "uni", 0); 202 ## Month resolution 203 times(ltime == 7) = ... 204 cellfun (@(x) [x "-00T00:00:00Z"], times(ltime == 7), "uni", 0); 205 ## Day resolution 206 times(ltime == 10) = ... 207 cellfun (@(x) [x "T00:00:00Z"], times(ltime == 10), "uni", 0); 208 endif 209 times = datenum (times, "yyyy-mm-ddTHH:MM:SSZ"); 210 ptrnP = '<gx:coord>(.+?) (.+?) (.+?)</gx:coord>'; 211 xyz = reshape (str2double (cell2mat (regexp (trk, ptrnP, "tokens"))), 3, [])'; 212 if (isempty (outp(end).Lat)) 213 outp(end).Lon = xyz(:, 1); 214 outp(end).Lat = xyz(:, 2); 215 outp(end).Ele = xyz(:, 3); 216 if (! isempty (times)) 217 outp(end).Time = times; 218 endif 219 elseif (intpm) 220 outp(end).Lon = [outp(end).Lon; xyz(:, 1)]; 221 outp(end).Lat = [outp(end).Lat; xyz(:, 2)]; 222 outp(end).Ele = [outp(end).Ele; xyz(:, 3)]; 223 if (! isempty (times)) 224 outp(end).Time = [outp(end).Time; times]; 225 endif 226 else 227 outp(end).Lon = [outp(end).Lon; NaN; xyz(:, 1)]; 228 outp(end).Lat = [outp(end).Lat; NaN; xyz(:, 2)]; 229 outp(end).Ele = [outp(end).Ele; NaN; xyz(:, 3)]; 230 if (! isempty (times)) 231 outp(end).Time = [outp(end).Time; NaN; times]; 232 endif 233 endif 234 outp(end).BoundingBox = [min(outp(end).Lon), max(outp(end).Lon); ... 235 min(outp(end).Lat), max(outp(end).Lat); ... 236 min(outp(end).Ele), max(outp(end).Ele)]; 237 endif 238 239 endwhile 240 241endfunction 242