File: sdts_tut.dox

package info (click to toggle)
gdal 1.10.1%2Bdfsg-8
  • links: PTS, VCS
  • area: main
  • in suites: jessie, jessie-kfreebsd
  • size: 84,320 kB
  • ctags: 74,726
  • sloc: cpp: 677,199; ansic: 162,820; python: 13,816; cs: 11,163; sh: 10,446; java: 5,279; perl: 4,429; php: 2,971; xml: 1,500; yacc: 934; makefile: 494; sql: 112
file content (469 lines) | stat: -rw-r--r-- 17,079 bytes parent folder | download | duplicates (6)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
/*! \page SDTS_AL_TUT

<center>
<title>SDTS Abstraction Library Tutorial</title>
</center>

This page is a walk through of the polygon layer portion of the 
<a href="sdts2shp.cpp.html">sdts2shp.cpp</a> example application.  It is should
give sufficient information to utilize the SDTS_AL library to read SDTS
files.<p>

<h2>Opening the Transfer</h2>

The following statements will open an SDTS transfer.  The filename
passed to SDTSTransfer::Open() should be the name of the catalog file,
such as <tt>palo_alto/SC01CATD.DDF</tt>.  The Open() method returns FALSE
if it fails for any reason.  In addition to the message we print out ourselves,
the SDTSTransfer::Open() method will also emit it's own error message using
CPLError().  See the cpl_error.h page for more information on how to 
capture and control CPLError() style error reporting.<p>

<pre>
#include "stds_al.h"

...

    SDTSTransfer oTransfer;

    if( !oTransfer.Open( pszCATDFilename ) )
    {
        fprintf( stderr,
                 "Failed to read CATD file `%s'\n",
                 pszCATDFilename );
        exit( 100 );
    }
</pre>

<h2>Getting a Layer List</h2>

Once an SDTSTransfer has been opened, it is possible to establish what
layers are available.  The sdts2shp example problem includes a -v argument
to dump a list of available layers.  It isn't normally necessary to use the
SDTS_CATD (catalog) from an application to access SDTS files; however, in 
this example we use it to fetch a module name, and description for each
of the available layers.<p>

In particular, the SDTSTransfer::GetLayerCount() method returns the
number of feature layers in the transfer and the 
SDTSTransfer::GetLayerCATDEntry() is used to translate layer indexes into
SDTS_CATD compatible CATD indexes.<p>

<pre>
        printf( "Layers:\n" );
        for( i = 0; i < oTransfer.GetLayerCount(); i++ )
        {
            int		iCATDEntry = oTransfer.GetLayerCATDEntry(i);
            
            printf( "  %s: `%s'\n",
                    oTransfer.GetCATD()->GetEntryModule(iCATDEntry),
                    oTransfer.GetCATD()->GetEntryTypeDesc(iCATDEntry) );
        }
        printf( "\n" );
</pre>

The following would be a typical layer list.  Note that there are many
other modules (files) registered with the catalog, but only these ones
are considered to be feature layers by the SDTSTransfer object.  The
rest are supporting information, much of it, like data quality, is ignored
by the SDTS_AL library.<p>

<pre>
warmerda@cs46980-c[113]% sdts2shp data/SC01CATD.DDF -v
Layers:
  ASCF: `Attribute Primary         '
  AHDR: `Attribute Primary         '
  NP01: `Point-Node                '
  NA01: `Point-Node                '
  NO01: `Point-Node                '
  LE01: `Line                      '
  PC01: `Polygon                   '
</pre>

<h2>Getting a Reader</h2>

In order to read polygon features, it is necessary to instantiate a polygon
reader on the desired layer.  The sdts2shp.cpp program allow the user to
select a module name (such as PC01, stored in pszMODN) to write to shape 
format.  Other application might just search for, and operate on all known
layers of a desired type.<p>

The SDTSTransfer::GetLayerIndexedReader() method instantiates a reader of
the desired type.  In this case we know we are instantiating a 
SDTSPolygonReader so we can safely cast the returned SDTSIndexedReader 
pointer to the more specific type SDTSPolygonReader.<p>

<pre>
    SDTSPolygonReader *poPolyReader;

    poPolyReader = (SDTSPolygonReader *) 
        poTransfer->GetLayerIndexedReader( poTransfer->FindLayer( pszMODN ) );
    
    if( poPolyReader == NULL )
    {
        fprintf( stderr, "Failed to open %s.\n",
                 poTransfer->GetCATD()->GetModuleFilePath( pszMODN ) );
        return;
    }
</pre>

Note that readers returned by SDTSTransfer::GetLayerIndexedReader() are
managed by the SDTSTransfer, and should not be deleted by the application.<p>

<h2>Collecting Polygon Geometry</h2>

The SDTS TVP format does not directly associate a polygons geometry (the
points forming it's boundary) with the polygon feature.  Instead it is
stored in separate line layers, and the lines contain references to the
right, and left polygons that the lines border.<p>

The SDTS_AL library provides a convenient method for forming the polygon
geometry.  Basically just call the SDTSPolygonReader::AssemblePolygons()
method.  This method will scan all SLTLine layers in the transfer, indexing
them and attaching their line work to the polygons.  Then it assembles the
line work into rings.  It also ensures that the outer ring comes first, that
the outer ring is counter-clockwise and that the inner ring(s) are
clockwise. 

<pre>
    poPolyReader->AssembleRings( poTransfer );
</pre>

Upon completion the SDTSPolygonReader will have been "indexed".  That means
that all the polygon information will have been read from disk, and the
polygon objects will now have information stored with them indicating the
list of edges that form their border.<p>

<h2>Identifying Attributes</h2>

In order to create the schema for the output shapefile dataset, it is
necessary to identify the attributes associated with the polygons.  There
are two types of attributes which can occur.  The first are hardcoded 
attributes specific to the feature type, and the second are generic 
user attributes stored in a separate primary attribute layer.<p>

In the case of SDTSRawPolygon, there is only one attribute of interest,
and that is the record number of the polygon.  This is actually stored within
the oModId data member of the SDTSIndexedFeature base class, as will be seen
in later examples when we write it to disk.  For now we create a DBF
field for the record number.  This record number is a unique identifier of
the polygon within this module/layer.<p>

<pre>
    nSDTSRecordField = DBFAddField( hDBF, "SDTSRecId", FTInteger, 8, 0 );
</pre>

Identification of user attributes is more complicated.  Any feature in a
layer can have associates with 0, 1, 2 or potentially more attribute records
in other primary attribute layers.  In order to establish a schema for the
layer it is necessary to build up a list of all attribute layers (tables) 
to which references appear.  The SDTSIndexedReader::ScanModuleReferences()
method can be used to scan a whole module for references to attribute modules
via the ATID field.  The return result is a list of referenced modules in the
form of a string list.  In a typical case this is one or two modules, such
as "ASCF".<p>  

<pre>
    char  **papszModRefs = poPolyReader->ScanModuleReferences();
</pre>

In sdts2shp.cpp, a subroutine (AddPrimaryAttrToDBFSchema()) is defined
to add all the fields of all references attribute layers to the DBF file. 
For each module in the list the following steps are executed.<p>

<h3>Fetch an Attribute Module Reader</h3>

The following code is similar to our code for create a polygon layer
reader.  It creates a reader on one of the attribute layers referenced.
We explicitly rewind it since it may have been previously opened and
read by another part of the application.<p>

<pre>
        SDTSAttrReader	*poAttrReader;

        poAttrReader = (SDTSAttrReader *)
            poTransfer->GetLayerIndexedReader(
                poTransfer->FindLayer( papszModuleList[iModule] ) );

        if( poAttrReader == NULL )
        {
            printf( "Unable to open attribute module %s, skipping.\n" ,
                    papszModuleList[iModule] );
            continue;
        }

        poAttrReader->Rewind();
</pre>

<h3>Get a Prototype Record</h3>

In order to get access to field definitions, and in order to establish
some sort of reasonable default lengths for field without fixed lengths
the sdts2shp program fetches a prototype record from the attribute module.

<pre>
        SDTSAttrRecord 	*poAttrFeature;

        poAttrFeature = (SDTSAttrRecord *) poAttrReader->GetNextFeature();
        if( poAttrFeature == NULL )
        {
            fprintf( stderr,
                     "Didn't find any meaningful attribute records in %s.\n",
                     papszModuleList[iModule] );
        
            continue;
        }
</pre>

When no longer needed, the attribute record may need to be explicitly 
deleted if it is not part of an indexed cached.<p>

<pre>
        if( !poAttrReader->IsIndexed() )
            delete poAttrFeature;
</pre>

<h3>Extract Field Definitions</h3>


The Shapefile DBF fields are defined based on the information available for
each of the subfields of the attribute records ATTR DDFField (the poATTR
data member).  The following code loops over each of the subfields, 
getting a pointer to the DDBSubfieldDefn containing information about that
subfield.<p>

<pre>
        DDFFieldDefn 	*poFDefn = poAttrFeature->poATTR->GetFieldDefn();
        int		iSF;
        DDFField	*poSR = poAttrFeature->poATTR;
    
        for( iSF=0; iSF < poFDefn->GetSubfieldCount(); iSF++ )
        {
            DDFSubfieldDefn	*poSFDefn = poFDefn->GetSubfield( iSF );
</pre>

Then each of the significant ISO8211 field types is translated to an
appropriate DBF field type.  In cases where the nWidth field is zero, 
indicating that the field is variable width, we use the length of the
field in the prototype record.  Ideally we would scan the whole file to find
the longest value for each field, but that would be a significant amount of
work. <p>

<pre>
            int		nWidth = poSFDefn->GetWidth();

            switch( poSFDefn->GetType() )
            {
              case DDFString:
                if( nWidth == 0 )
                {
                    int		nMaxBytes;
                
                    const char * pachData = poSR->GetSubfieldData(poSFDefn,
                                                                  &nMaxBytes);

                    nWidth = strlen(poSFDefn->ExtractStringData(pachData,
                                                                nMaxBytes, NULL ));
                }
            
                DBFAddField( hDBF, poSFDefn->GetName(), FTString, nWidth, 0 );
                break;

              case DDFInt:
                if( nWidth == 0 )
                    nWidth = 9;

                DBFAddField( hDBF, poSFDefn->GetName(), FTInteger, nWidth, 0 );
                break;

              case DDFFloat:
                DBFAddField( hDBF, poSFDefn->GetName(), FTDouble, 18, 6 );
                break;

              default:
                fprintf( stderr,
                         "Dropping attribute `%s' of module `%s'.  "
                         "Type unsupported\n",
                         poSFDefn->GetName(),
                         papszModuleList[iModule] );
                break;
            }
        }
</pre>

<h2>Reading Polygon Features</h2>

With definition of the attribute schema out of the way, we return to the
main event, reading polygons from the polygon layer.  We have already
instantiated the SDTSPolygonReader (poPolyReader), and now we loop reading
features from it.  Note that we Rewind() the reader to ensure we are
starting at the beginning.  After we are done process the polygon we
delete it, if and only if the layer does not have an index cache.<p>

<pre>
    SDTSRawPolygon	*poRawPoly;
        
    poPolyReader->Rewind();
    while( (poRawPoly = (SDTSRawPolygon *) poPolyReader->GetNextFeature())
           != NULL )
    {
        ... process and write polygon ...

        if( !poPolyReader->IsIndexed() )
            delete poRawPoly;
    }
</pre>

<h2>Translate Geometry</h2>

In an earlier step we used the SDTSPolygonReader::AssembleRings() method to
build ring geometry on the polygons from the linework in the line layers.<p>

Coincidently (well, ok, maybe it isn't a coincidence) it so happens that the 
ring organization exactly matches what is needed for the shapefile api.  
The following call creates a polygon from the ring information in the
SDTSRawPolygon.  See the SDTSRawPolygon reference help for a fuller 
definition of the nRings, panRingStart, nVertices, and vertex fields.<p>

<pre>
        psShape = SHPCreateObject( SHPT_POLYGON, -1, poRawPoly->nRings,
                                   poRawPoly->panRingStart, NULL,
                                   poRawPoly->nVertices,
                                   poRawPoly->padfX, 
                                   poRawPoly->padfY, 
                                   poRawPoly->padfZ,
                                   NULL );
</pre>

<h2>Write Record Number</h2>

The following call is used to write out the record number of the polygon,
fetched from the SDTSIndexedFeature::oModId data member.  The szModule value
in this data field will always match the module name for the whole layer. 
While not shown here, there is also an szOBRP field on oModId which have
different values depending on whether the polygon is a universe or regular
polygon.<p>

<pre>
        DBFWriteIntegerAttribute( hDBF, iShape, nSDTSRecordField,
                                  poRawPoly->oModId.nRecord );
</pre>

<h2>Fetch Associated User Records</h2>

In keeping with the setting up of the schema, accessing the user records
is somewhat complicated.  In sdts2shp, the primary attribute records associated
with any feature (including SDTSRawPolygons) can be fetched with the 
WriteAttrRecordToDBF() function defined as follows.<p>

In particular, the poFeature->nAttributes member indicates how many 
associated attribute records there are.  The poFeature->aoATID[] array
contains the SDTSModId's for each record.  This SDTSModId can be passed
to SDTSTransfer::GetAttr() to fetch the DDFField pointer for the user
attributes.  The WriteAttrRecordToDBF() method is specific to sdts2shp
and will be define later.<p>

<pre>
    int		iAttrRecord;
    
    for( iAttrRecord = 0; iAttrRecord < poFeature->nAttributes; iAttrRecord++)
    {
        DDFField	*poSR;

        poSR = poTransfer->GetAttr( poFeature->aoATID+iAttrRecord );
          
        WriteAttrRecordToDBF( hDBF, iRecord, poTransfer, poSR );
    }
</pre>

<h2>Write User Attributes</h2>

In a manner analygous to the definition of the fields from the prototype
attribute record, the following code loops over the subfields, and fetches
the data for each.  The data extraction via poSR->GetSubfieldData() is
a bit involved, and more information can be found on the DDFField reference
page.<p>

<pre>
/* -------------------------------------------------------------------- */
/*      Process each subfield in the record.                            */
/* -------------------------------------------------------------------- */
    DDFFieldDefn 	*poFDefn = poSR->GetFieldDefn();
        
    for( int iSF=0; iSF < poFDefn->GetSubfieldCount(); iSF++ )
    {
        DDFSubfieldDefn	*poSFDefn = poFDefn->GetSubfield( iSF );
        int			iField;
        int			nMaxBytes;
        const char * 	pachData = poSR->GetSubfieldData(poSFDefn,
                                                         &nMaxBytes);

/* -------------------------------------------------------------------- */
/*      Identify the related DBF field, if any.                         */
/* -------------------------------------------------------------------- */
        for( iField = 0; iField < hDBF->nFields; iField++ )
        {
            if( EQUALN(poSFDefn->GetName(),
                       hDBF->pszHeader+iField*32,10) )
                break;
        }

        if( iField == hDBF->nFields )
            iField = -1;
            
/* -------------------------------------------------------------------- */
/*      Handle each of the types.                                       */
/* -------------------------------------------------------------------- */
        switch( poSFDefn->GetType() )
        {
          case DDFString:
            const char	*pszValue;

            pszValue = poSFDefn->ExtractStringData(pachData, nMaxBytes,
                                                   NULL);

            if( iField != -1 )
                DBFWriteStringAttribute(hDBF, iRecord, iField, pszValue );
            break;

          case DDFFloat:
            double	dfValue;

            dfValue = poSFDefn->ExtractFloatData(pachData, nMaxBytes,
                                                 NULL);

            if( iField != -1 )
                DBFWriteDoubleAttribute( hDBF, iRecord, iField, dfValue );
            break;

          case DDFInt:
            int		nValue;

            nValue = poSFDefn->ExtractIntData(pachData, nMaxBytes, NULL);

            if( iField != -1 )
                DBFWriteIntegerAttribute( hDBF, iRecord, iField, nValue );
            break;

          default:
            break;
        }
    } /* next subfield */
</pre>

<h2>Cleanup</h2>

In the case of sdts2shp, the SDTSTransfer is created on the stack.  When it
falls out of scope it is destroyed, and all the indexed readers, and their
indexed features caches are also cleaned up.<p>

*/

/*!
\page sdts2shp.cpp
<center>
<title>SDTS To Shape Example Application</title>
</center>

\include sdts2shp.cpp
*/