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