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