1 /* -*-c++-*- */
2 /* osgEarth - Geospatial SDK for OpenSceneGraph
3 * Copyright 2019 Pelican Mapping
4 * http://osgearth.org
5 *
6 * osgEarth is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
12 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
13 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
14 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
15 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
16 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
17 * IN THE SOFTWARE.
18 *
19 * You should have received a copy of the GNU Lesser General Public License
20 * along with this program.  If not, see <http://www.gnu.org/licenses/>
21 */
22 #include <osgEarth/Notify>
23 #include <osgEarth/TDTiles>
24 #include <osgEarthFeatures/FeatureSource>
25 #include <osgEarthFeatures/FeatureCursor>
26 #include <osgEarthFeatures/ResampleFilter>
27 #include <osgEarthDrivers/feature_ogr/OGRFeatureOptions>
28 #include <osg/ArgumentParser>
29 #include <osgDB/FileUtils>
30 #include <fstream>
31 
32 #define LC "[featuretiler] "
33 
34 using namespace osgEarth;
35 using namespace osgEarth::Features;
36 using namespace osgEarth::Drivers;
37 
38 int
39 usage(const char* msg)
40 {
41     OE_NOTICE
42         << "\n" << msg
43         << "\nUsage: osgearth_featuretiler"
44         << "\n    --in [location]        ; filename of source data (e.g. a shapefile)"
45         << "\n    --out [directory]      ; output folder name"
46         << "\n    --errors [...]         ; Comma-delimited geometric error per level (e.g. 1,150,500)"
usage(const std::string & app)47         << std::endl;
48     return -1;
49 }
50 
51 struct Env
52 {
53     osg::ref_ptr<FeatureSource> input;
54     double* geometricError;
55     unsigned maxDepth;
56     URIContext uriContext;
57     unsigned counter;
58 
59     ~Env() { delete [] geometricError; }
60 };
61 
62 struct FeatureData {
63     double x, y;
64     FeatureID fid;
main(int argc,char ** argv)65 };
66 
67 struct SortByX {
68     bool operator()(const FeatureData& lhs, const FeatureData& rhs) const {
69         return lhs.x < rhs.x;
70     }
71 };
72 
73 struct SortByY {
74     bool operator()(const FeatureData& lhs, const FeatureData& rhs) const {
75         return lhs.y < rhs.y;
76     }
77 };
78 
79 int
80 split(const GeoExtent& extent, TDTiles::Tile* parentTile, unsigned depth, Env& env)
81 {
82     const Profile* profile = 0L;
83     const FeatureProfile* inputFP = env.input->getFeatureProfile();
84     if (inputFP && inputFP->getExtent().isValid())
85     {
86         const GeoExtent& fex = inputFP->getExtent();
87         profile = Profile::create(fex.getSRS(), fex.xMin(), fex.yMin(), fex.xMax(), fex.yMax());
88     }
89 
90     OGRFeatureOptions outputConf[2];
91     osg::ref_ptr<FeatureSource> output[2];
92 
93     for (unsigned i = 0; i < 2; ++i)
94     {
95         outputConf[i].url() = URI(Stringify() << "test_" << depth << "_" << (env.counter++) << "_" << i << ".shp", env.uriContext);
96         outputConf[i].openWrite() = true;
97         if (profile)
98             outputConf[i].profile() = profile->toProfileOptions();
99 
100         output[i] = FeatureSourceFactory::create(outputConf[i]);
101         if (!output[i].valid())
102             return usage("Failed to open output dataset");
103 
104         const Status& outputStatus = output[i]->create(
105             env.input->getFeatureProfile(),
106             env.input->getSchema(),
107             env.input->getGeometryType(),
108             NULL);
109 
110         if (outputStatus.isError())
111             return usage(outputStatus.message().c_str());
112     }
113 
114     std::vector<FeatureData> data;
115     int count = env.input->getFeatureCount();
116     if (count > 0)
117     {
118         //OE_INFO << "Feature count = " << count << std::endl;
119         data.reserve(count);
120     }
121     else
122     {
123         OE_INFO << "Feature count not available" << std::endl;
124     }
125 
126     // Just divide the error by 10 each level. Placeholder for something smarter.
127     double error = env.geometricError[depth];
128 
129     FeatureList features;
130     unsigned i = 0;
131     Query query;
132     query.bounds() = extent.bounds();
133     osg::ref_ptr<FeatureCursor> cursor = env.input->createFeatureCursor(query, 0L);
134     if (!cursor.valid())
135     {
136         OE_WARN << "Feature cursor error" << std::endl;
137         return 0;
138     }
139 
140     cursor->fill(features);
141     if (features.empty())
142         return 0;
143 
144     ResampleFilter resample(error, DBL_MAX);
145     FilterContext fc(0L, env.input->getFeatureProfile(), extent);
146     resample.push(features, fc);
147 
148     for(FeatureList::iterator i = features.begin(); i != features.end(); ++i)
149     {
150         Feature* f = i->get();
151         const GeoExtent& fex = f->getExtent();
152         if (fex.width() < error && fex.height() < error)
153             continue;
154 
155         FeatureData d;
156         fex.getCentroid(d.x, d.y);
157         d.fid = f->getFID();
158         data.push_back(d);
159     }
160 
161     bool isWide = extent.width() > extent.height();
162 
163     OE_INFO << "Depth = " << depth << ", features = " << data.size() << std::endl;
164 
165     double median;
166 
167     if (data.size() > 0)
168     {
169         if (isWide)
170         {
171             SortByX sortByX;
172             std::sort(data.begin(), data.end(), sortByX);
173             median = ((data.size() & 0x1) == 0) ?
174                 0.5 * (data[data.size() / 2].x + data[data.size() / 2].x) :
175                 data[data.size() / 2].x;
176             OE_INFO << "  Median X = " << median << std::endl;
177         }
178         else
179         {
180             SortByY sortByY;
181             std::sort(data.begin(), data.end(), sortByY);
182             median = ((data.size() & 0x1) == 0) ?
183                 0.5 * (data[data.size() / 2].y + data[data.size() / 2].y) :
184                 data[data.size() / 2].y;
185             OE_INFO << "  Median Y = " << median << std::endl;
186         }
187     }
188     else
189     {
190         if (isWide)
191             median = extent.xMin() + extent.width()/2.0;
192         else
193             median = extent.yMin() + extent.height()/2.0;
194     }
195 
196     osg::ref_ptr<TDTiles::Tile> child[2];
197 
198     child[0] = new TDTiles::Tile();
199     child[0]->refine() = TDTiles::REFINE_REPLACE;
200     parentTile->children().push_back(child[0].get());
201 
202     child[1] = new TDTiles::Tile();
203     child[0]->refine() = TDTiles::REFINE_REPLACE;
204     parentTile->children().push_back(child[1].get());
205 
206     unsigned afeatures = 0u, bfeatures = 0u;
207 
208     for (unsigned i = 0; i < data.size(); ++i)
209     {
210         const FeatureData& d = data[i];
211         Feature* f = env.input->getFeature(d.fid);
212         double x, y;
213         f->getExtent().getCentroid(x, y);
214         double side = isWide ? x : y;
215         if (side < median)
216         {
217             afeatures++;
218             output[0]->insertFeature(f);
219         }
220         else
221         {
222             bfeatures++;
223             output[1]->insertFeature(f);
224         }
225     }
226 
227     GeoExtent a(extent.getSRS()), b(extent.getSRS());
228 
229     if (isWide)
230     {
231         a.set(extent.xMin(), extent.yMin(), median, extent.yMax());
232         b.set(median, extent.yMin(), extent.xMax(), extent.yMax());
233     }
234     else
235     {
236         a.set(extent.xMin(), extent.yMin(), extent.xMax(), median);
237         b.set(extent.xMin(), median, extent.xMax(), extent.yMax());
238     }
239 
240     if (depth < env.maxDepth)
241     {
242         split(a, child[0].get(), depth + 1, env);
243         split(b, child[1].get(), depth + 1, env);
244     }
245 
246     output[0] = 0L;
247     output[1] = 0L;
248 
249     //TODO: heights
250     const float zmin = 0.0f;
251     const float zmax = 1.0f;
252 
253     const SpatialReference* epsg4979 = SpatialReference::get("epsg:4979");
254     a = a.transform(epsg4979);
255     b = b.transform(epsg4979);
256 
257     if (afeatures > 0u)
258         child[0]->content()->uri() = outputConf[0].url().get();
259 
260     child[0]->boundingVolume()->region()->set(
261         osg::DegreesToRadians(a.xMin()), osg::DegreesToRadians(a.yMin()), zmin,
262         osg::DegreesToRadians(a.xMax()), osg::DegreesToRadians(a.yMax()), zmax);
263     child[0]->geometricError() = error;
264 
265     if (bfeatures > 0u)
266         child[1]->content()->uri() = outputConf[1].url().get();
267 
268     child[1]->boundingVolume()->region()->set(
269         osg::DegreesToRadians(b.xMin()), osg::DegreesToRadians(b.yMin()), zmin,
270         osg::DegreesToRadians(b.xMax()), osg::DegreesToRadians(b.yMax()), zmax);
271     child[1]->geometricError() = error;
272 
273     return 0;
274 }
275 
276 int
277 main(int argc, char** argv)
278 {
279     osg::ArgumentParser arguments(&argc,argv);
280 
281     if (arguments.read("--help"))
282         return usage("Help!");
283 
284     std::string source;
285     if (!arguments.read("--in", source))
286         return usage("Missing required --in argument");
287 
288     std::string outputLocation;
289     if (!arguments.read("--out", outputLocation))
290         return usage("Missing required --out argument");
291 
292     std::string errors("1,80,200");
293     arguments.read("--errors", errors);
294 
295     if (!osgDB::makeDirectory(outputLocation))
296         return usage("Unable to create/find output location");
297 
298     osg::ref_ptr<TDTiles::Tile> root = new TDTiles::Tile();
299 
300     std::string outFile = Stringify() << outputLocation << "/tileset.json";
301     URIContext uriContext(outFile);
302 
303     Config inputConf;
304     inputConf.set("driver", "ogr");
305     inputConf.set("url", source);
306     osg::ref_ptr<FeatureSource> input = FeatureSourceFactory::create(ConfigOptions(inputConf));
307     if (!input.valid())
308         return usage("Failed to create the input feature source");
309 
310     if (input->open().isError())
311         return usage("Failed to open the input feature source");
312 
313     std::vector<std::string> errorStrings;
314     StringTokenizer(errors, errorStrings, ",");
315     if (errorStrings.size() < 1)
316         return usage("Illegal errors input");
317 
318     Env env;
319     env.input = input.get();
320     env.counter = 0;
321     env.uriContext = uriContext;
322     env.geometricError = new double[errorStrings.size()];
323     env.maxDepth = errorStrings.size()-1;
324 
325     // generate level errors:
326     OE_NOTICE << "Geometric errors: ";
327     for(int i=0; i<errorStrings.size(); ++i)
328     {
329         int j = errorStrings.size() - 1 - i;
330         env.geometricError[j] = osgEarth::as<double>(errorStrings[i], -1.0);
331         if (env.geometricError[j] <= 0.0)
332             return usage("Illegal geometric error in input");
333 
334         OE_NOTIFY(osg::NOTICE, env.geometricError[j] << " ");
335     }
336     OE_NOTIFY(osg::NOTICE, std::endl);
337 
338     // make it so
339     split(input->getFeatureProfile()->getExtent(), root, 0, env);
340 
341     osg::ref_ptr<TDTiles::Tileset> tileset = new TDTiles::Tileset();
342     tileset->root() = root.get();
343 
344     TDTiles::Asset asset;
345     tileset->asset()->version() = "1.0";
346 
347     std::ofstream out(outFile.c_str());
348     Json::Value tilesetJSON = tileset->getJSON();
349     Json::StyledStreamWriter writer;
350     writer.write(out, tilesetJSON);
351     out.close();
352 
353     return 0;
354 }
355