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