Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: SOSI Data Source
4 : * Purpose: Provide SOSI Data to OGR.
5 : * Author: Thomas Hirsch, <thomas.hirsch statkart no>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2010, Thomas Hirsch
9 : * Copyright (c) 2010-2013, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "ogr_sosi.h"
31 : #include <map>
32 : #include <math.h>
33 :
34 : /* This is the most common encoding for SOSI files. Let's at least try if
35 : * it is supported, or generate a meaningful error message. */
36 : #ifndef CPL_ENC_ISO8859_10
37 : #define CPL_ENC_ISO8859_10 "ISO8859-10"
38 : #endif
39 :
40 : #ifdef WRITE_SUPPORT
41 : /************************************************************************/
42 : /* utility methods */
43 : /************************************************************************/
44 :
45 : static int epsg2sosi(int nEPSG)
46 : {
47 : int nSOSI = 23;
48 : switch (nEPSG)
49 : {
50 : case 27391: /* NGO 1984 Axis I-VIII */
51 : case 27392:
52 : case 27393:
53 : case 27394:
54 : case 27395:
55 : case 27396:
56 : case 27397:
57 : case 27398:
58 : {
59 : nSOSI = nEPSG - 27390;
60 : break;
61 : }
62 : case 3043: /* UTM ZONE 31-36 */
63 : case 3044:
64 : case 3045:
65 : case 3046:
66 : case 3047:
67 : case 3048:
68 : {
69 : nSOSI = nEPSG - 3022;
70 : break;
71 : }
72 : case 23031: /* UTM ZONE 31-36 / ED50 */
73 : case 23032:
74 : case 23033:
75 : case 23034:
76 : case 23035:
77 : case 23036:
78 : {
79 : nSOSI = nEPSG - 23000;
80 : break;
81 : }
82 : case 4326:
83 : { /* WSG84 */
84 : nSOSI = 84;
85 : break;
86 : }
87 : default:
88 : {
89 : CPLError(CE_Warning, CPLE_AppDefined,
90 : "(Yet) unsupported coordinate system writing to SOSI "
91 : "file: %i. Defaulting to EPSG:4326/SOSI 84.",
92 : nEPSG);
93 : }
94 : }
95 : return nSOSI;
96 : }
97 : #endif
98 :
99 3 : static int sosi2epsg(int nSOSI)
100 : {
101 3 : int nEPSG = 4326;
102 3 : switch (nSOSI)
103 : {
104 0 : case 1: /* NGO 1984 Axis I-VIII */
105 : case 2:
106 : case 3:
107 : case 4:
108 : case 5:
109 : case 6:
110 : case 7:
111 : case 8:
112 : {
113 0 : nEPSG = 27390 + nSOSI;
114 0 : break;
115 : }
116 3 : case 21: /* UTM ZONE 31-36 */
117 : case 22:
118 : case 23:
119 : case 24:
120 : case 25:
121 : case 26:
122 : {
123 3 : nEPSG = 3022 + nSOSI;
124 3 : break;
125 : }
126 0 : case 31: /* UTM ZONE 31-36 / ED50 */
127 : case 32:
128 : case 33:
129 : case 34:
130 : case 35:
131 : case 36:
132 : {
133 0 : nEPSG = 23000 + nSOSI;
134 0 : break;
135 : }
136 0 : case 84:
137 : { /* WSG84 */
138 0 : nEPSG = 4326;
139 0 : break;
140 : }
141 0 : default:
142 : {
143 0 : CPLError(CE_Warning, CPLE_AppDefined,
144 : "(Yet) unsupported coordinate system in SOSI-file: %i. "
145 : "Defaulting to EPSG:4326.",
146 : nSOSI);
147 : }
148 : }
149 3 : return nEPSG;
150 : }
151 :
152 : /************************************************************************/
153 : /* OGRSOSIDataSource() */
154 : /************************************************************************/
155 :
156 3 : OGRSOSIDataSource::OGRSOSIDataSource()
157 : {
158 3 : nLayers = 0;
159 :
160 3 : poFileadm = nullptr;
161 3 : poBaseadm = nullptr;
162 3 : papoBuiltGeometries = nullptr;
163 3 : papoLayers = nullptr;
164 3 : pszName = nullptr;
165 3 : poSRS = nullptr;
166 :
167 3 : poPolyHeaders = nullptr;
168 3 : poTextHeaders = nullptr;
169 3 : poPointHeaders = nullptr;
170 3 : poCurveHeaders = nullptr;
171 :
172 3 : pszEncoding = CPL_ENC_UTF8;
173 3 : nNumFeatures = 0;
174 :
175 3 : nMode = MODE_READING;
176 3 : }
177 :
178 : /************************************************************************/
179 : /* ~OGRSOSIDataSource() */
180 : /************************************************************************/
181 :
182 6 : OGRSOSIDataSource::~OGRSOSIDataSource()
183 : {
184 3 : if (papoBuiltGeometries != nullptr)
185 : {
186 60 : for (unsigned int i = 0; i < nNumFeatures; i++)
187 : {
188 57 : if (papoBuiltGeometries[i] != nullptr)
189 : {
190 51 : delete papoBuiltGeometries[i];
191 51 : papoBuiltGeometries[i] = nullptr;
192 : }
193 : }
194 3 : CPLFree(papoBuiltGeometries);
195 3 : papoBuiltGeometries = nullptr;
196 : }
197 :
198 3 : if (poPolyHeaders != nullptr)
199 3 : delete poPolyHeaders;
200 3 : if (poTextHeaders != nullptr)
201 0 : delete poTextHeaders;
202 3 : if (poPointHeaders != nullptr)
203 0 : delete poPointHeaders;
204 3 : if (poCurveHeaders != nullptr)
205 3 : delete poCurveHeaders;
206 :
207 3 : if (nMode == MODE_WRITING)
208 : {
209 0 : if (poFileadm != nullptr)
210 0 : LC_CloseSos(poFileadm, RESET_IDX);
211 0 : if (poBaseadm != nullptr)
212 0 : LC_CloseBase(poBaseadm, RESET_IDX);
213 : }
214 : else
215 : {
216 3 : if (poFileadm != nullptr)
217 3 : LC_CloseSos(poFileadm, SAVE_IDX);
218 3 : if (poBaseadm != nullptr)
219 3 : LC_CloseBase(poBaseadm, SAVE_IDX);
220 : }
221 3 : poFileadm = nullptr;
222 3 : poBaseadm = nullptr;
223 :
224 3 : if (papoLayers != nullptr)
225 : {
226 9 : for (int i = 0; i < nLayers; i++)
227 : {
228 6 : delete papoLayers[i];
229 : }
230 3 : CPLFree(papoLayers);
231 : }
232 :
233 3 : if (poSRS != nullptr)
234 3 : poSRS->Release();
235 3 : if (pszName != nullptr)
236 3 : CPLFree(pszName);
237 6 : }
238 :
239 6 : static OGRFeatureDefn *defineLayer(const char *szName,
240 : OGRwkbGeometryType szType, S2I *poHeaders,
241 : S2I **ppoHeadersNew)
242 : {
243 6 : OGRFeatureDefn *poFeatureDefn = new OGRFeatureDefn(szName);
244 6 : poFeatureDefn->SetGeomType(szType);
245 6 : S2I *poHeadersNew = *ppoHeadersNew;
246 :
247 57 : for (S2I::iterator i = poHeaders->begin(); i != poHeaders->end(); ++i)
248 : {
249 51 : OGRSOSIDataType *poType = SOSIGetType(i->first);
250 51 : OGRSOSISimpleDataType *poElements = poType->getElements();
251 132 : for (int k = 0; k < poType->getElementCount(); k++)
252 : {
253 81 : if (strcmp(poElements[k].GetName(), "") == 0)
254 9 : continue;
255 72 : OGRFieldDefn oFieldTemplate(poElements[k].GetName(),
256 144 : poElements[k].GetType());
257 72 : (*poHeadersNew)[CPLString(poElements[k].GetName())] =
258 72 : poFeatureDefn->GetFieldCount();
259 72 : poFeatureDefn->AddFieldDefn(&oFieldTemplate);
260 : }
261 : }
262 6 : return poFeatureDefn;
263 : }
264 :
265 : /************************************************************************/
266 : /* Open() */
267 : /************************************************************************/
268 :
269 3 : int OGRSOSIDataSource::Open(const char *pszFilename, int bUpdate)
270 : {
271 3 : papoBuiltGeometries = nullptr;
272 3 : poFileadm = nullptr;
273 3 : poBaseadm = nullptr;
274 : char *pszPos;
275 :
276 3 : if (bUpdate)
277 : {
278 0 : CPLError(CE_Failure, CPLE_OpenFailed,
279 : "Update access not supported by the SOSI driver.");
280 0 : return FALSE;
281 : }
282 :
283 : /* Check that the file exists otherwise HO_TestSOSI() emits an error */
284 : VSIStatBuf sStat;
285 3 : if (VSIStat(pszFilename, &sStat) != 0)
286 0 : return FALSE;
287 :
288 3 : pszName = CPLStrdup(pszFilename);
289 : /* We ignore any layer parameters for now. */
290 3 : pszPos = strchr(pszName, ',');
291 3 : if (pszPos != nullptr)
292 : {
293 0 : pszPos[0] = '\0';
294 : }
295 :
296 : /* Confirm that we are dealing with a SOSI file. Used also by data
297 : * format auto-detection in some ogr utilities. */
298 3 : UT_INT64 nEnd = 0;
299 3 : int bIsSosi = HO_TestSOSI(pszName, &nEnd);
300 3 : if (bIsSosi == UT_FALSE)
301 : {
302 0 : return FALSE; /* No error message: This is used by file format
303 : auto-detection */
304 : }
305 :
306 3 : short nStatus = 0, nDetStatus = 0; /* immediate status, detailed status */
307 :
308 : /* open index base and sosi file */
309 3 : poBaseadm = LC_OpenBase(LC_BASE);
310 3 : nStatus = LC_OpenSos(pszName, LC_BASE_FRAMGR, LC_GML_IDX, LC_INGEN_STATUS,
311 : &poFileadm, &nDetStatus);
312 3 : if (nStatus == UT_FALSE)
313 : {
314 : char *pszErrorMessage;
315 0 : LC_StrError(nDetStatus, &pszErrorMessage);
316 0 : CPLError(CE_Failure, CPLE_OpenFailed,
317 : "File %s could not be opened by SOSI Driver: %s", pszName,
318 : pszErrorMessage);
319 0 : return FALSE;
320 : }
321 :
322 : /* --------------------------------------------------------------------*
323 : * Prefetch all the information needed to determine layers *
324 : * and prebuild LineString features for later assembly. *
325 : * --------------------------------------------------------------------*/
326 :
327 : /* allocate room for one pointer per feature */
328 3 : nNumFeatures = static_cast<unsigned int>(poFileadm->lAntGr);
329 3 : papoBuiltGeometries = static_cast<OGRGeometry **>(
330 3 : VSI_MALLOC2_VERBOSE(nNumFeatures, sizeof(OGRGeometry *)));
331 3 : if (papoBuiltGeometries == nullptr)
332 : {
333 0 : nNumFeatures = 0;
334 0 : return FALSE;
335 : }
336 60 : for (unsigned int i = 0; i < nNumFeatures; i++)
337 57 : papoBuiltGeometries[i] = nullptr;
338 :
339 : /* Various iterators and return values used to iterate through SOSI features
340 : */
341 : short nName, nNumLines;
342 : long nNumCoo;
343 : unsigned short nInfo;
344 : LC_SNR_ADM oSnradm;
345 : LC_BGR oNextSerial;
346 : LC_BGR *poNextSerial;
347 3 : poNextSerial = &oNextSerial;
348 :
349 3 : bool bPointLayer =
350 : FALSE; /* Initialize four layers for the different geometry types */
351 3 : bool bCurveLayer = FALSE;
352 3 : bool bPolyLayer = FALSE;
353 3 : bool bTextLayer = FALSE;
354 3 : poPolyHeaders = new S2I();
355 3 : poPointHeaders = new S2I();
356 3 : poCurveHeaders = new S2I();
357 3 : poTextHeaders = new S2I();
358 :
359 3 : LC_SBSn(&oSnradm, poFileadm, 0, nNumFeatures); /* Set FYBA search limits */
360 3 : LC_InitNextBgr(poNextSerial);
361 :
362 : /* Prebuilding simple features and extracting layer information. */
363 60 : while (LC_NextBgr(poNextSerial, LC_FRAMGR))
364 : {
365 : /* Fetch next group information */
366 : nName =
367 57 : LC_RxGr(poNextSerial, LES_OPTIMALT, &nNumLines, &nNumCoo, &nInfo);
368 :
369 57 : S2S oHeaders;
370 57 : S2S::iterator iHeaders;
371 : int iH;
372 : /* Extract all strings from group header. */
373 510 : for (short i = 1; i <= nNumLines; i++)
374 : {
375 453 : char *pszLine = LC_GetGi(i); /* Get one header line */
376 453 : if ((pszLine[0] == ':') || (pszLine[0] == '('))
377 3 : continue; /* If we have a continued REF line, skip it. */
378 450 : if (pszLine[0] == '!')
379 0 : continue; /* If we have a comment line, skip it. */
380 :
381 450 : char *pszUTFLine = CPLRecode(
382 : pszLine, pszEncoding, CPL_ENC_UTF8); /* switch to UTF encoding
383 : here, if it is known. */
384 450 : char *pszUTFLineIter = pszUTFLine;
385 1308 : while (pszUTFLineIter[0] == '.')
386 858 : pszUTFLineIter++; /* Skipping the dots at the beginning of a
387 : SOSI line */
388 : char *pszPos2 =
389 450 : strstr(pszUTFLineIter, " "); /* Split header and value */
390 450 : if (pszPos2 != nullptr)
391 : {
392 882 : const std::string osKey(pszUTFLineIter, pszPos2);
393 441 : oHeaders[osKey] =
394 882 : CPLString(pszPos2 + 1); /* Add to header map */
395 441 : switch (nName)
396 : { /* Add to header list for the corresponding layer, if it is
397 : not */
398 27 : case L_FLATE:
399 : { /* in there already */
400 27 : if (poPolyHeaders->find(osKey) == poPolyHeaders->end())
401 : {
402 24 : iH = static_cast<int>(poPolyHeaders->size());
403 24 : (*poPolyHeaders)[osKey] = iH;
404 : }
405 27 : break;
406 : }
407 387 : case L_KURVE:
408 : case L_LINJE:
409 : case L_BUEP:
410 : { /* FIXME: maybe not use the same headers for both */
411 387 : if (poCurveHeaders->find(osKey) ==
412 774 : poCurveHeaders->end())
413 : {
414 27 : iH = static_cast<int>(poCurveHeaders->size());
415 27 : (*poCurveHeaders)[osKey] = iH;
416 : }
417 387 : break;
418 : }
419 0 : case L_PUNKT:
420 : case L_SYMBOL:
421 : {
422 0 : if (poPointHeaders->find(osKey) ==
423 0 : poPointHeaders->end())
424 : {
425 0 : iH = static_cast<int>(poPointHeaders->size());
426 0 : (*poPointHeaders)[osKey] = iH;
427 : }
428 0 : break;
429 : }
430 0 : case L_TEKST:
431 : {
432 0 : if (poTextHeaders->find(osKey) == poTextHeaders->end())
433 : {
434 0 : iH = static_cast<int>(poTextHeaders->size());
435 0 : (*poTextHeaders)[osKey] = iH;
436 : }
437 0 : break;
438 : }
439 : }
440 : }
441 450 : CPLFree(pszUTFLine);
442 : }
443 :
444 : /* Feature-specific tasks */
445 57 : switch (nName)
446 : {
447 0 : case L_PUNKT:
448 : {
449 : /* Pre-build a point feature. Activate point layer. */
450 0 : bPointLayer = TRUE;
451 0 : buildOGRPoint(oNextSerial.lNr);
452 0 : break;
453 : }
454 3 : case L_FLATE:
455 : {
456 : /* Activate polygon layer. */
457 3 : bPolyLayer = TRUE;
458 : /* cannot build geometries that reference others yet */
459 3 : break;
460 : }
461 51 : case L_KURVE:
462 : case L_LINJE:
463 : {
464 : /* Pre-build a line feature. Activate line/curve layer. */
465 51 : bCurveLayer = TRUE;
466 51 : buildOGRLineString(static_cast<int>(nNumCoo), oNextSerial.lNr);
467 51 : break;
468 : }
469 0 : case L_BUEP:
470 : {
471 : /* Pre-build a line feature as interpolation from an arc.
472 : * Activate line/curve layer. */
473 0 : bCurveLayer = TRUE;
474 0 : buildOGRLineStringFromArc(oNextSerial.lNr);
475 0 : break;
476 : }
477 0 : case L_TEKST:
478 : {
479 : /* Pre-build a text line contour feature. Activate text layer.
480 : */
481 : /* Todo: observe only points 2ff if more than one point is given
482 : * for follow mode */
483 0 : bTextLayer = TRUE;
484 0 : buildOGRMultiPoint(static_cast<int>(nNumCoo), oNextSerial.lNr);
485 0 : break;
486 : }
487 3 : case L_HODE:
488 : {
489 : /* Get SRS from SOSI header. */
490 3 : unsigned short nMask = LC_TR_ALLT;
491 : LC_TRANSPAR oTrans;
492 3 : if (LC_GetTransEx(&nMask, &oTrans) == UT_FALSE)
493 : {
494 0 : CPLError(CE_Failure, CPLE_OpenFailed,
495 : "TRANSPAR section not found - No reference system "
496 : "information available.");
497 0 : return FALSE;
498 : }
499 3 : poSRS = new OGRSpatialReference();
500 3 : poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
501 :
502 : /* Get coordinate system from SOSI header. */
503 3 : int nEPSG = sosi2epsg(oTrans.sKoordsys);
504 3 : if (poSRS->importFromEPSG(nEPSG) != OGRERR_NONE)
505 : {
506 0 : CPLError(CE_Failure, CPLE_OpenFailed,
507 : "OGR could not load coordinate system definition "
508 : "EPSG:%i.",
509 : nEPSG);
510 0 : return FALSE;
511 : }
512 :
513 : /* Get character encoding from SOSI header. */
514 3 : iHeaders = oHeaders.find("TEGNSETT");
515 3 : if (iHeaders != oHeaders.end())
516 : {
517 6 : CPLString osLine = iHeaders->second;
518 3 : if (osLine.compare("ISO8859-1") == 0)
519 : {
520 3 : pszEncoding = CPL_ENC_ISO8859_1;
521 : }
522 0 : else if (osLine.compare("ISO8859-10") == 0)
523 : {
524 0 : pszEncoding = CPL_ENC_ISO8859_10;
525 : }
526 0 : else if (osLine.compare("UTF-8") == 0)
527 : {
528 0 : pszEncoding = CPL_ENC_UTF8;
529 : }
530 : }
531 :
532 3 : break;
533 : }
534 0 : default:
535 : {
536 0 : break;
537 : }
538 : }
539 : }
540 :
541 : /* -------------------------------------------------------------------- *
542 : * Create a corresponding layers. One per geometry type *
543 : * -------------------------------------------------------------------- */
544 3 : int l_nLayers = 0;
545 3 : if (bPolyLayer)
546 3 : l_nLayers++;
547 3 : if (bCurveLayer)
548 3 : l_nLayers++;
549 3 : if (bPointLayer)
550 0 : l_nLayers++;
551 3 : if (bTextLayer)
552 0 : l_nLayers++;
553 3 : this->nLayers = l_nLayers;
554 : /* allocate some memory for up to three layers */
555 3 : papoLayers =
556 3 : (OGRSOSILayer **)VSI_MALLOC2_VERBOSE(sizeof(void *), l_nLayers);
557 3 : if (papoLayers == nullptr)
558 : {
559 0 : this->nLayers = 0;
560 0 : return FALSE;
561 : }
562 :
563 : /* Define each layer, using a proper feature definition, geometry type,
564 : * and adding every SOSI header encountered in the file as field. */
565 3 : if (bPolyLayer)
566 : {
567 3 : S2I *poHeadersNew = new S2I();
568 : OGRFeatureDefn *poFeatureDefn =
569 3 : defineLayer("polygons", wkbPolygon, poPolyHeaders, &poHeadersNew);
570 3 : delete poPolyHeaders;
571 3 : poPolyHeaders = poHeadersNew;
572 3 : poFeatureDefn->Reference();
573 3 : papoLayers[--l_nLayers] =
574 3 : new OGRSOSILayer(this, poFeatureDefn, poFileadm, poPolyHeaders);
575 : }
576 : else
577 : {
578 0 : delete poPolyHeaders;
579 0 : poPolyHeaders = nullptr;
580 : }
581 3 : if (bCurveLayer)
582 : {
583 3 : S2I *poHeadersNew = new S2I();
584 : OGRFeatureDefn *poFeatureDefn =
585 3 : defineLayer("lines", wkbLineString, poCurveHeaders, &poHeadersNew);
586 3 : delete poCurveHeaders;
587 3 : poCurveHeaders = poHeadersNew;
588 3 : poFeatureDefn->Reference();
589 3 : papoLayers[--l_nLayers] =
590 3 : new OGRSOSILayer(this, poFeatureDefn, poFileadm, poCurveHeaders);
591 : }
592 : else
593 : {
594 0 : delete poCurveHeaders;
595 0 : poCurveHeaders = nullptr;
596 : }
597 3 : if (bPointLayer)
598 : {
599 0 : S2I *poHeadersNew = new S2I();
600 : OGRFeatureDefn *poFeatureDefn =
601 0 : defineLayer("points", wkbPoint, poPointHeaders, &poHeadersNew);
602 0 : delete poPointHeaders;
603 0 : poPointHeaders = poHeadersNew;
604 0 : poFeatureDefn->Reference();
605 0 : papoLayers[--l_nLayers] =
606 0 : new OGRSOSILayer(this, poFeatureDefn, poFileadm, poPointHeaders);
607 : }
608 : else
609 : {
610 3 : delete poPointHeaders;
611 3 : poPointHeaders = nullptr;
612 : }
613 3 : if (bTextLayer)
614 : {
615 0 : S2I *poHeadersNew = new S2I();
616 : OGRFeatureDefn *poFeatureDefn =
617 0 : defineLayer("text", wkbMultiPoint, poTextHeaders, &poHeadersNew);
618 0 : delete poTextHeaders;
619 0 : poTextHeaders = poHeadersNew;
620 0 : poFeatureDefn->Reference();
621 0 : papoLayers[--l_nLayers] =
622 0 : new OGRSOSILayer(this, poFeatureDefn, poFileadm, poTextHeaders);
623 : }
624 : else
625 : {
626 3 : delete poTextHeaders;
627 3 : poTextHeaders = nullptr;
628 : }
629 3 : return TRUE;
630 : }
631 :
632 : #ifdef WRITE_SUPPORT
633 : /************************************************************************/
634 : /* Create() */
635 : /************************************************************************/
636 :
637 : int OGRSOSIDataSource::Create(const char *pszFilename)
638 : {
639 : short nStatus;
640 : short nDetStatus;
641 :
642 : poBaseadm = LC_OpenBase(LC_KLADD);
643 : nStatus = LC_OpenSos(pszFilename, LC_SEKV_SKRIV, LC_NY_IDX, LC_INGEN_STATUS,
644 : &poFileadm, &nDetStatus);
645 : if (nStatus == UT_FALSE)
646 : {
647 : CPLError(CE_Failure, CPLE_OpenFailed,
648 : "Could not open SOSI file for writing (Status %i).",
649 : nDetStatus);
650 : return FALSE;
651 : }
652 :
653 : LC_NyttHode(); /* Create new file header, will be written to file when all
654 : header information elements are set. */
655 :
656 : return TRUE;
657 : }
658 :
659 : /************************************************************************/
660 : /* ICreateLayer() */
661 : /************************************************************************/
662 :
663 : OGRLayer *OGRSOSIDataSource::ICreateLayer(
664 : const char *pszNameIn, const OGRSpatialReference *poSpatialRef,
665 : OGRwkbGeometryType eGType, CPL_UNUSED char **papszOptions)
666 : {
667 : /* SOSI does not really support layers - so let's first see that the global
668 : * settings are consistent */
669 : if (poSRS == NULL)
670 : {
671 : if (poSpatialRef != NULL)
672 : {
673 : poSRS = poSpatialRef->Clone();
674 :
675 : const char *pszKoosys = poSRS->GetAuthorityCode("PROJCS");
676 : if (pszKoosys == NULL)
677 : {
678 : OGRErr err = poSRS->AutoIdentifyEPSG();
679 : if (err == OGRERR_UNSUPPORTED_SRS)
680 : {
681 : CPLError(CE_Failure, CPLE_OpenFailed,
682 : "Could not identify EPSG code for spatial "
683 : "reference system");
684 : return NULL;
685 : }
686 : pszKoosys = poSRS->GetAuthorityCode("PROJCS");
687 : }
688 :
689 : if (pszKoosys != NULL)
690 : {
691 : int nKoosys = epsg2sosi(atoi(pszKoosys));
692 : CPLDebug("[CreateLayer]", "Projection set to SOSI %i", nKoosys);
693 : LC_PutTrans(nKoosys, 0, 0, 0.01, 0.01, 0.01);
694 : }
695 : else
696 : {
697 : pszKoosys = poSRS->GetAuthorityCode("GEOGCS");
698 : if (pszKoosys != NULL)
699 : {
700 : int nKoosys = epsg2sosi(atoi(pszKoosys));
701 : LC_PutTrans(nKoosys, 0, 0, 0.01, 0.01, 0.01);
702 : }
703 : else
704 : {
705 : CPLError(CE_Failure, CPLE_OpenFailed,
706 : "Could not retrieve EPSG code for spatial "
707 : "reference system");
708 : return NULL;
709 : }
710 : }
711 : }
712 : LC_WsGr(poFileadm); /* Writing the header here! */
713 : }
714 : else
715 : {
716 : if (!poSRS->IsSame(poSpatialRef))
717 : {
718 : CPLError(CE_Failure, CPLE_AppDefined,
719 : "SOSI driver does not support different spatial reference "
720 : "systems in one file.");
721 : }
722 : }
723 :
724 : OGRFeatureDefn *poFeatureDefn = new OGRFeatureDefn(pszNameIn);
725 : poFeatureDefn->Reference();
726 : poFeatureDefn->SetGeomType(eGType);
727 : OGRSOSILayer *poLayer =
728 : new OGRSOSILayer(this, poFeatureDefn, poFileadm, NULL /*poHeaderDefn*/);
729 : /* todo: where do we delete poLayer and poFeatureDefn? */
730 : return poLayer;
731 : }
732 : #endif
733 :
734 : /************************************************************************/
735 : /* GetLayer() */
736 : /************************************************************************/
737 :
738 6 : OGRLayer *OGRSOSIDataSource::GetLayer(int iLayer)
739 : {
740 6 : if (iLayer < 0 || iLayer >= nLayers)
741 0 : return nullptr;
742 : else
743 6 : return papoLayers[iLayer];
744 : }
745 :
746 0 : void OGRSOSIDataSource::buildOGRMultiPoint(int nNumCoo, long iSerial)
747 : {
748 0 : if (papoBuiltGeometries[iSerial] != nullptr)
749 : {
750 0 : return;
751 : }
752 :
753 0 : OGRMultiPoint *poMP = new OGRMultiPoint();
754 :
755 : long i;
756 0 : double dfEast = 0.0;
757 0 : double dfNorth = 0.0;
758 0 : for (i = (nNumCoo > 1) ? 2 : 1; i <= nNumCoo; i++)
759 : {
760 0 : LC_GetTK(i, &dfEast, &dfNorth);
761 0 : OGRPoint poP = OGRPoint(dfEast, dfNorth);
762 0 : poMP->addGeometry(&poP); /*poP will be cloned before returning*/
763 : }
764 0 : papoBuiltGeometries[iSerial] = poMP;
765 : }
766 :
767 51 : void OGRSOSIDataSource::buildOGRLineString(int nNumCoo, long iSerial)
768 : {
769 51 : if (papoBuiltGeometries[iSerial] != nullptr)
770 : {
771 0 : return;
772 : }
773 :
774 51 : OGRLineString *poLS = new OGRLineString();
775 51 : poLS->setNumPoints(nNumCoo);
776 :
777 : int i;
778 51 : double dfEast = 0, dfNorth = 0;
779 516 : for (i = 1; i <= nNumCoo; i++)
780 : {
781 465 : LC_GetTK(i, &dfEast, &dfNorth);
782 465 : poLS->setPoint(i - 1, dfEast, dfNorth);
783 : }
784 51 : papoBuiltGeometries[iSerial] = poLS;
785 : }
786 :
787 0 : static double sqr(double x)
788 : {
789 0 : return x * x;
790 : }
791 :
792 0 : void OGRSOSIDataSource::buildOGRLineStringFromArc(long iSerial)
793 : {
794 0 : if (papoBuiltGeometries[iSerial] != nullptr)
795 : {
796 0 : return;
797 : }
798 :
799 0 : OGRLineString *poLS = new OGRLineString();
800 :
801 : /* fetch reference points on circle (easting, northing) */
802 0 : double e1 = 0, e2 = 0, e3 = 0;
803 0 : double n1 = 0, n2 = 0, n3 = 0;
804 0 : LC_GetTK(1, &e1, &n1);
805 0 : LC_GetTK(2, &e2, &n2);
806 0 : LC_GetTK(3, &e3, &n3);
807 :
808 : /* helper constants */
809 0 : double p12 = (e1 * e1 - e2 * e2 + n1 * n1 - n2 * n2) / 2;
810 0 : double p13 = (e1 * e1 - e3 * e3 + n1 * n1 - n3 * n3) / 2;
811 :
812 0 : double dE12 = e1 - e2;
813 0 : double dE13 = e1 - e3;
814 0 : double dN12 = n1 - n2;
815 0 : double dN13 = n1 - n3;
816 :
817 : /* center of the circle */
818 0 : double cE = (dN13 * p12 - dN12 * p13) / (dE12 * dN13 - dN12 * dE13);
819 0 : double cN = (dE13 * p12 - dE12 * p13) / (dN12 * dE13 - dE12 * dN13);
820 :
821 : /* radius of the circle */
822 0 : double r = sqrt(sqr(e1 - cE) + sqr(n1 - cN));
823 :
824 : /* angles of points A and B (1 and 3) */
825 0 : double th1 = atan2(n1 - cN, e1 - cE);
826 0 : double th3 = atan2(n3 - cN, e3 - cE);
827 :
828 : /* interpolation step in radians */
829 0 : double dth = th3 - th1;
830 0 : if (dth < 0)
831 : {
832 0 : dth += 2 * M_PI;
833 : }
834 0 : if (dth > M_PI)
835 : {
836 0 : dth = -2 * M_PI + dth;
837 : }
838 0 : int npt = (int)(ARC_INTERPOLATION_FULL_CIRCLE * dth / 2 * M_PI);
839 0 : if (npt < 0)
840 0 : npt = -npt;
841 0 : if (npt < 3)
842 0 : npt = 3;
843 0 : poLS->setNumPoints(npt);
844 0 : dth = dth / (npt - 1);
845 :
846 0 : for (int i = 0; i < npt; i++)
847 : {
848 0 : const double dfEast = cE + r * cos(th1 + dth * i);
849 0 : const double dfNorth = cN + r * sin(th1 + dth * i);
850 0 : if (dfEast != dfEast)
851 : { /* which is a wonderful property of nans */
852 0 : CPLError(CE_Warning, CPLE_AppDefined,
853 : "Calculated %lf for point %d of %d in curve %li.", dfEast,
854 : i, npt, iSerial);
855 : }
856 0 : poLS->setPoint(i, dfEast, dfNorth);
857 : }
858 0 : papoBuiltGeometries[iSerial] = poLS;
859 : }
860 :
861 0 : void OGRSOSIDataSource::buildOGRPoint(long iSerial)
862 : {
863 0 : double dfEast = 0, dfNorth = 0;
864 0 : LC_GetTK(1, &dfEast, &dfNorth);
865 0 : papoBuiltGeometries[iSerial] = new OGRPoint(dfEast, dfNorth);
866 0 : }
867 :
868 : /************************************************************************/
869 : /* TestCapability() */
870 : /************************************************************************/
871 :
872 0 : int OGRSOSIDataSource::TestCapability(CPL_UNUSED const char *pszCap)
873 : {
874 : #ifdef WRITE_SUPPORT
875 : if (strcmp("CreateLayer", pszCap) == 0)
876 : {
877 : return TRUE;
878 : }
879 : #endif
880 0 : return FALSE;
881 : }
|
