1 /******************************************************************************
2 * $Id: ogrdodsgrid.cpp 28375 2015-01-30 12:06:11Z rouault $
3 *
4 * Project: OGR/DODS Interface
5 * Purpose: Implements OGRDODSGridLayer class, which implements the
6 * "Grid/Array" access strategy.
7 * Author: Frank Warmerdam, warmerdam@pobox.com
8 *
9 ******************************************************************************
10 * Copyright (c) 2004, Frank Warmerdam <warmerdam@pobox.com>
11 *
12 * Permission is hereby granted, free of charge, to any person obtaining a
13 * copy of this software and associated documentation files (the "Software"),
14 * to deal in the Software without restriction, including without limitation
15 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 * and/or sell copies of the Software, and to permit persons to whom the
17 * Software is furnished to do so, subject to the following conditions:
18 *
19 * The above copyright notice and this permission notice shall be included
20 * in all copies or substantial portions of the Software.
21 *
22 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 * DEALINGS IN THE SOFTWARE.
29 ****************************************************************************/
30
31 #include "cpl_conv.h"
32 #include "ogr_dods.h"
33 #include "cpl_string.h"
34
35 CPL_CVSID("$Id: ogrdodsgrid.cpp 28375 2015-01-30 12:06:11Z rouault $");
36
37 /************************************************************************/
38 /* OGRDODSGridLayer() */
39 /************************************************************************/
40
OGRDODSGridLayer(OGRDODSDataSource * poDSIn,const char * pszTargetIn,AttrTable * poOGRLayerInfoIn)41 OGRDODSGridLayer::OGRDODSGridLayer( OGRDODSDataSource *poDSIn,
42 const char *pszTargetIn,
43 AttrTable *poOGRLayerInfoIn )
44
45 : OGRDODSLayer( poDSIn, pszTargetIn, poOGRLayerInfoIn )
46
47 {
48 pRawData = NULL;
49
50 /* -------------------------------------------------------------------- */
51 /* What is the layer name? */
52 /* -------------------------------------------------------------------- */
53 string oLayerName;
54 const char *pszLayerName = pszTargetIn;
55
56 if( poOGRLayerInfo != NULL )
57 {
58 oLayerName = poOGRLayerInfo->get_attr( "layer_name" );
59 if( strlen(oLayerName.c_str()) > 0 )
60 pszLayerName = oLayerName.c_str();
61 }
62
63 poFeatureDefn = new OGRFeatureDefn( pszLayerName );
64 poFeatureDefn->Reference();
65
66 /* -------------------------------------------------------------------- */
67 /* Fetch the target variable. */
68 /* -------------------------------------------------------------------- */
69 BaseType *poTargVar = poDS->poDDS->var( pszTargetIn );
70
71 if( poTargVar->type() == dods_grid_c )
72 {
73 poTargetGrid = dynamic_cast<Grid *>( poTargVar );
74 poTargetArray = dynamic_cast<Array *>(poTargetGrid->array_var());
75 }
76 else if( poTargVar->type() == dods_array_c )
77 {
78 poTargetGrid = NULL;
79 poTargetArray = dynamic_cast<Array *>( poTargVar );
80 }
81 else
82 {
83 CPLAssert( FALSE );
84 return;
85 }
86
87 /* -------------------------------------------------------------------- */
88 /* Count arrays in use. */
89 /* -------------------------------------------------------------------- */
90 AttrTable *poExtraContainers = NULL;
91 nArrayRefCount = 1; // primary target.
92
93 if( poOGRLayerInfo != NULL )
94 poExtraContainers = poOGRLayerInfo->find_container("extra_containers");
95
96 if( poExtraContainers != NULL )
97 {
98 AttrTable::Attr_iter dv_i;
99
100 for( dv_i = poExtraContainers->attr_begin();
101 dv_i != poExtraContainers->attr_end(); dv_i++ )
102 {
103 nArrayRefCount++;
104 }
105 }
106
107 /* -------------------------------------------------------------------- */
108 /* Collect extra_containers. */
109 /* -------------------------------------------------------------------- */
110 paoArrayRefs = new OGRDODSArrayRef[nArrayRefCount];
111 paoArrayRefs[0].pszName = CPLStrdup( pszTargetIn );
112 paoArrayRefs[0].poArray = poTargetArray;
113
114 nArrayRefCount = 1;
115
116 if( poExtraContainers != NULL )
117 {
118 AttrTable::Attr_iter dv_i;
119
120 for( dv_i = poExtraContainers->attr_begin();
121 dv_i != poExtraContainers->attr_end(); dv_i++ )
122 {
123 const char *pszTargetName=poExtraContainers->get_attr(dv_i).c_str();
124 BaseType *poExtraTarget = poDS->poDDS->var( pszTargetName );
125
126 if( poExtraTarget == NULL )
127 {
128 CPLError( CE_Warning, CPLE_AppDefined,
129 "Unable to find extra_container '%s', skipping.",
130 pszTargetName );
131 continue;
132 }
133
134 if( poExtraTarget->type() == dods_array_c )
135 paoArrayRefs[nArrayRefCount].poArray =
136 dynamic_cast<Array *>( poExtraTarget );
137 else if( poExtraTarget->type() == dods_grid_c )
138 {
139 Grid *poGrid = dynamic_cast<Grid *>( poExtraTarget );
140 paoArrayRefs[nArrayRefCount].poArray =
141 dynamic_cast<Array *>( poGrid->array_var() );
142 }
143 else
144 {
145 CPLError( CE_Warning, CPLE_AppDefined,
146 "Target container '%s' is not grid or array, skipping.",
147 pszTargetName );
148 continue;
149 }
150
151 paoArrayRefs[nArrayRefCount++].pszName = CPLStrdup(pszTargetName);
152 }
153 }
154
155 /* -------------------------------------------------------------------- */
156 /* Collect dimension information. */
157 /* -------------------------------------------------------------------- */
158 int iDim;
159 Array::Dim_iter iterDim;
160
161 nDimCount = poTargetArray->dimensions();
162 paoDimensions = new OGRDODSDim[nDimCount];
163 nMaxRawIndex = 1;
164
165 for( iterDim = poTargetArray->dim_begin(), iDim = 0;
166 iterDim != poTargetArray->dim_end();
167 iterDim++, iDim++ )
168 {
169 paoDimensions[iDim].pszDimName =
170 CPLStrdup(poTargetArray->dimension_name(iterDim).c_str());
171 paoDimensions[iDim].nDimStart =
172 poTargetArray->dimension_start(iterDim);
173 paoDimensions[iDim].nDimEnd =
174 poTargetArray->dimension_stop(iterDim);
175 paoDimensions[iDim].nDimStride =
176 poTargetArray->dimension_stride(iterDim);
177 paoDimensions[iDim].poMap = NULL;
178
179 paoDimensions[iDim].nDimEntries =
180 (paoDimensions[iDim].nDimEnd + 1 - paoDimensions[iDim].nDimStart
181 + paoDimensions[iDim].nDimStride - 1)
182 / paoDimensions[iDim].nDimStride;
183
184 nMaxRawIndex *= paoDimensions[iDim].nDimEntries;
185 }
186
187 /* -------------------------------------------------------------------- */
188 /* If we are working with a grid, collect the maps. */
189 /* -------------------------------------------------------------------- */
190 if( poTargetGrid != NULL )
191 {
192 int iMap;
193 Grid::Map_iter iterMap;
194
195 for( iterMap = poTargetGrid->map_begin(), iMap = 0;
196 iterMap != poTargetGrid->map_end();
197 iterMap++, iMap++ )
198 {
199 paoDimensions[iMap].poMap = dynamic_cast<Array *>(*iterMap);
200 }
201
202 CPLAssert( iMap == nDimCount );
203 }
204
205 /* -------------------------------------------------------------------- */
206 /* Setup field definitions. The first nDimCount will be the */
207 /* dimension attributes, and after that comes the actual target */
208 /* array. */
209 /* -------------------------------------------------------------------- */
210 for( iDim = 0; iDim < nDimCount; iDim++ )
211 {
212 OGRFieldDefn oField( paoDimensions[iDim].pszDimName, OFTInteger );
213
214 if( EQUAL(oField.GetNameRef(), poTargetArray->name().c_str()) )
215 oField.SetName(CPLSPrintf("%s_i",paoDimensions[iDim].pszDimName));
216
217 if( paoDimensions[iDim].poMap != NULL )
218 {
219 switch( paoDimensions[iDim].poMap->var()->type() )
220 {
221 case dods_byte_c:
222 case dods_int16_c:
223 case dods_uint16_c:
224 case dods_int32_c:
225 case dods_uint32_c:
226 oField.SetType( OFTInteger );
227 break;
228
229 case dods_float32_c:
230 case dods_float64_c:
231 oField.SetType( OFTReal );
232 break;
233
234 case dods_str_c:
235 case dods_url_c:
236 oField.SetType( OFTString );
237 break;
238
239 default:
240 // Ignore
241 break;
242 }
243 }
244
245 poFeatureDefn->AddFieldDefn( &oField );
246 }
247
248 /* -------------------------------------------------------------------- */
249 /* Setup the array attributes themselves. */
250 /* -------------------------------------------------------------------- */
251 int iArray;
252 for( iArray=0; iArray < nArrayRefCount; iArray++ )
253 {
254 OGRDODSArrayRef *poRef = paoArrayRefs + iArray;
255 OGRFieldDefn oArrayField( poRef->poArray->name().c_str(), OFTInteger );
256
257 switch( poRef->poArray->var()->type() )
258 {
259 case dods_byte_c:
260 case dods_int16_c:
261 case dods_uint16_c:
262 case dods_int32_c:
263 case dods_uint32_c:
264 oArrayField.SetType( OFTInteger );
265 break;
266
267 case dods_float32_c:
268 case dods_float64_c:
269 oArrayField.SetType( OFTReal );
270 break;
271
272 case dods_str_c:
273 case dods_url_c:
274 oArrayField.SetType( OFTString );
275 break;
276
277 default:
278 // Ignore
279 break;
280 }
281
282 poFeatureDefn->AddFieldDefn( &oArrayField );
283 poRef->iFieldIndex = poFeatureDefn->GetFieldCount() - 1;
284 }
285
286 /* -------------------------------------------------------------------- */
287 /* X/Y/Z fields. */
288 /* -------------------------------------------------------------------- */
289 if( poOGRLayerInfo != NULL )
290 {
291 AttrTable *poField = poOGRLayerInfo->find_container("x_field");
292 if( poField != NULL )
293 {
294 oXField.Initialize( poField );
295 oXField.iFieldIndex =
296 poFeatureDefn->GetFieldIndex( oXField.pszFieldName );
297 }
298
299 poField = poOGRLayerInfo->find_container("y_field");
300 if( poField != NULL )
301 {
302 oYField.Initialize( poField );
303 oYField.iFieldIndex =
304 poFeatureDefn->GetFieldIndex( oYField.pszFieldName );
305 }
306
307 poField = poOGRLayerInfo->find_container("z_field");
308 if( poField != NULL )
309 {
310 oZField.Initialize( poField );
311 oZField.iFieldIndex =
312 poFeatureDefn->GetFieldIndex( oZField.pszFieldName );
313 }
314
315 }
316
317 /* -------------------------------------------------------------------- */
318 /* If we have no layerinfo, then check if there are obvious x/y */
319 /* fields. */
320 /* -------------------------------------------------------------------- */
321 else
322 {
323 if( poFeatureDefn->GetFieldIndex( "lat" ) != -1
324 && poFeatureDefn->GetFieldIndex( "lon" ) != -1 )
325 {
326 oXField.Initialize( "lon", "dds" );
327 oXField.iFieldIndex = poFeatureDefn->GetFieldIndex( "lon" );
328 oYField.Initialize( "lat", "dds" );
329 oYField.iFieldIndex = poFeatureDefn->GetFieldIndex( "lat" );
330 }
331 else if( poFeatureDefn->GetFieldIndex( "latitude" ) != -1
332 && poFeatureDefn->GetFieldIndex( "longitude" ) != -1 )
333 {
334 oXField.Initialize( "longitude", "dds" );
335 oXField.iFieldIndex = poFeatureDefn->GetFieldIndex( "longitude" );
336 oYField.Initialize( "latitude", "dds" );
337 oYField.iFieldIndex = poFeatureDefn->GetFieldIndex( "latitude" );
338 }
339 }
340
341 /* -------------------------------------------------------------------- */
342 /* Set the layer geometry type if we have point inputs. */
343 /* -------------------------------------------------------------------- */
344 if( oZField.iFieldIndex >= 0 )
345 poFeatureDefn->SetGeomType( wkbPoint25D );
346 else if( oXField.iFieldIndex >= 0 && oYField.iFieldIndex >= 0 )
347 poFeatureDefn->SetGeomType( wkbPoint );
348 else
349 poFeatureDefn->SetGeomType( wkbNone );
350 }
351
352 /************************************************************************/
353 /* ~OGRDODSGridLayer() */
354 /************************************************************************/
355
~OGRDODSGridLayer()356 OGRDODSGridLayer::~OGRDODSGridLayer()
357
358 {
359 delete[] paoArrayRefs;
360 delete[] paoDimensions;
361 }
362
363 /************************************************************************/
364 /* ArrayEntryToField() */
365 /************************************************************************/
366
ArrayEntryToField(Array * poArray,void * pRawData,int iArrayIndex,OGRFeature * poFeature,int iField)367 int OGRDODSGridLayer::ArrayEntryToField( Array *poArray, void *pRawData,
368 int iArrayIndex,
369 OGRFeature *poFeature, int iField)
370
371 {
372 switch( poArray->var()->type() )
373 {
374 case dods_byte_c:
375 {
376 GByte *pabyRawData = (GByte *) pRawData;
377 poFeature->SetField( iField, pabyRawData[iArrayIndex] );
378 }
379 break;
380
381 case dods_int16_c:
382 {
383 GInt16 *panRawData = (GInt16 *) pRawData;
384 poFeature->SetField( iField, panRawData[iArrayIndex] );
385 }
386 break;
387
388 case dods_uint16_c:
389 {
390 GUInt16 *panRawData = (GUInt16 *) pRawData;
391 poFeature->SetField( iField, panRawData[iArrayIndex] );
392 }
393 break;
394
395 case dods_int32_c:
396 {
397 GInt32 *panRawData = (GInt32 *) pRawData;
398 poFeature->SetField( iField, panRawData[iArrayIndex] );
399 }
400 break;
401
402 case dods_uint32_c:
403 {
404 GUInt32 *panRawData = (GUInt32 *) pRawData;
405 poFeature->SetField( iField, (int) panRawData[iArrayIndex] );
406 }
407 break;
408
409 case dods_float32_c:
410 {
411 float * pafRawData = (float *) pRawData;
412 poFeature->SetField( iField, pafRawData[iArrayIndex] );
413 }
414 break;
415
416 case dods_float64_c:
417 {
418 double * padfRawData = (double *) pRawData;
419 poFeature->SetField( iField, padfRawData[iArrayIndex] );
420 }
421 break;
422
423 default:
424 return FALSE;
425 }
426
427 return TRUE;
428 }
429
430 /************************************************************************/
431 /* GetFeature() */
432 /************************************************************************/
433
GetFeature(GIntBig nFeatureId)434 OGRFeature *OGRDODSGridLayer::GetFeature( GIntBig nFeatureId )
435
436 {
437 if( nFeatureId < 0 || nFeatureId >= nMaxRawIndex )
438 return NULL;
439
440 /* -------------------------------------------------------------------- */
441 /* Ensure we have the dataset. */
442 /* -------------------------------------------------------------------- */
443 if( !ProvideDataDDS() )
444 return NULL;
445
446 /* -------------------------------------------------------------------- */
447 /* Create the feature being read. */
448 /* -------------------------------------------------------------------- */
449 OGRFeature *poFeature;
450
451 poFeature = new OGRFeature( poFeatureDefn );
452 poFeature->SetFID( nFeatureId );
453 m_nFeaturesRead++;
454
455 /* -------------------------------------------------------------------- */
456 /* Establish the values for the various dimension indices. */
457 /* -------------------------------------------------------------------- */
458 int iDim;
459 int nRemainder = nFeatureId;
460
461 for( iDim = nDimCount-1; iDim >= 0; iDim-- )
462 {
463 paoDimensions[iDim].iLastValue =
464 (nRemainder % paoDimensions[iDim].nDimEntries)
465 * paoDimensions[iDim].nDimStride
466 + paoDimensions[iDim].nDimStart;
467 nRemainder = nRemainder / paoDimensions[iDim].nDimEntries;
468
469 if( poTargetGrid == NULL )
470 poFeature->SetField( iDim, paoDimensions[iDim].iLastValue );
471 }
472 CPLAssert( nRemainder == 0 );
473
474 /* -------------------------------------------------------------------- */
475 /* For grids, we need to apply the values of the dimensions */
476 /* looked up in the corresponding map. */
477 /* -------------------------------------------------------------------- */
478 if( poTargetGrid != NULL )
479 {
480 for( iDim = 0; iDim < nDimCount; iDim++ )
481 {
482 ArrayEntryToField( paoDimensions[iDim].poMap,
483 paoDimensions[iDim].pRawData,
484 paoDimensions[iDim].iLastValue,
485 poFeature, iDim );
486 }
487 }
488
489 /* -------------------------------------------------------------------- */
490 /* Process all the regular data fields. */
491 /* -------------------------------------------------------------------- */
492 int iArray;
493 for( iArray = 0; iArray < nArrayRefCount; iArray++ )
494 {
495 OGRDODSArrayRef *poRef = paoArrayRefs + iArray;
496
497 ArrayEntryToField( poRef->poArray, poRef->pRawData, nFeatureId,
498 poFeature, poRef->iFieldIndex );
499 }
500
501 /* -------------------------------------------------------------------- */
502 /* Do we have geometry information? */
503 /* -------------------------------------------------------------------- */
504 if( oXField.iFieldIndex >= 0 && oYField.iFieldIndex >= 0 )
505 {
506 OGRPoint *poPoint = new OGRPoint();
507
508 poPoint->setX( poFeature->GetFieldAsDouble( oXField.iFieldIndex ) );
509 poPoint->setY( poFeature->GetFieldAsDouble( oYField.iFieldIndex ) );
510
511 if( oZField.iFieldIndex >= 0 )
512 poPoint->setZ( poFeature->GetFieldAsDouble(oZField.iFieldIndex) );
513
514 poFeature->SetGeometryDirectly( poPoint );
515 }
516
517 return poFeature;
518 }
519
520 /************************************************************************/
521 /* ProvideDataDDS() */
522 /************************************************************************/
523
ProvideDataDDS()524 int OGRDODSGridLayer::ProvideDataDDS()
525
526 {
527 if( bDataLoaded )
528 return poTargetVar != NULL;
529
530 int bResult = OGRDODSLayer::ProvideDataDDS();
531
532 if( !bResult )
533 return bResult;
534
535 int iArray;
536 for( iArray=0; iArray < nArrayRefCount; iArray++ )
537 {
538 OGRDODSArrayRef *poRef = paoArrayRefs + iArray;
539 BaseType *poTarget = poDataDDS->var( poRef->pszName );
540
541 // Reset ref array pointer to point in DataDDS result.
542 if( poTarget->type() == dods_grid_c )
543 {
544 Grid *poGrid = dynamic_cast<Grid *>( poTarget );
545 poRef->poArray = dynamic_cast<Array *>(poGrid->array_var());
546
547 if( iArray == 0 )
548 poTargetGrid = poGrid;
549 }
550 else if( poTarget->type() == dods_array_c )
551 {
552 poRef->poArray = dynamic_cast<Array *>( poTarget );
553 }
554 else
555 {
556 CPLAssert( FALSE );
557 return FALSE;
558 }
559
560 if( iArray == 0 )
561 poTargetArray = poRef->poArray;
562
563 // Allocate appropriate raw data array, and pull out data into it.
564 poRef->pRawData = CPLMalloc( poRef->poArray->width() );
565 poRef->poArray->buf2val( &(poRef->pRawData) );
566 }
567
568 // Setup pointers to each of the map objects.
569 if( poTargetGrid != NULL )
570 {
571 int iMap;
572 Grid::Map_iter iterMap;
573
574 for( iterMap = poTargetGrid->map_begin(), iMap = 0;
575 iterMap != poTargetGrid->map_end();
576 iterMap++, iMap++ )
577 {
578 paoDimensions[iMap].poMap = dynamic_cast<Array *>(*iterMap);
579 paoDimensions[iMap].pRawData =
580 CPLMalloc( paoDimensions[iMap].poMap->width() );
581 paoDimensions[iMap].poMap->buf2val( &(paoDimensions[iMap].pRawData) );
582 }
583 }
584
585 return bResult;
586 }
587
588 /************************************************************************/
589 /* GetFeatureCount() */
590 /************************************************************************/
591
GetFeatureCount(int bForce)592 GIntBig OGRDODSGridLayer::GetFeatureCount( int bForce )
593
594 {
595 if( m_poFilterGeom == NULL && m_poAttrQuery == NULL )
596 return nMaxRawIndex;
597 else
598 return OGRDODSLayer::GetFeatureCount( bForce );
599 }
600