1 : /******************************************************************************
2 : * $Id: netcdfdataset.cpp 18175 2009-12-05 00:03:36Z rouault $
3 : *
4 : * Project: netCDF read/write Driver
5 : * Purpose: GDAL bindings over netCDF library.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2004, Frank Warmerdam
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 "netcdfdataset.h"
31 : #include "cpl_error.h"
32 : CPL_CVSID("$Id: netcdfdataset.cpp 18175 2009-12-05 00:03:36Z rouault $");
33 :
34 :
35 : /************************************************************************/
36 : /* ==================================================================== */
37 : /* netCDFRasterBand */
38 : /* ==================================================================== */
39 : /************************************************************************/
40 :
41 : class netCDFRasterBand : public GDALPamRasterBand
42 : {
43 : nc_type nc_datatype;
44 : int nZId;
45 : int nZDim;
46 : int nLevel;
47 : int nBandXPos;
48 : int nBandYPos;
49 : int *panBandZPos;
50 : int *panBandZLev;
51 : int bNoDataSet;
52 : double dfNoDataValue;
53 : CPLErr CreateBandMetadata( );
54 :
55 : public:
56 :
57 : netCDFRasterBand( netCDFDataset *poDS,
58 : int nZId,
59 : int nZDim,
60 : int nLevel,
61 : int *panBandZLen,
62 : int *panBandPos,
63 : int nBand );
64 : ~netCDFRasterBand( );
65 : virtual double GetNoDataValue( int * );
66 : virtual CPLErr SetNoDataValue( double );
67 : virtual CPLErr IReadBlock( int, int, void * );
68 :
69 :
70 : };
71 :
72 : /************************************************************************/
73 : /* GetMetadata() */
74 : /************************************************************************/
75 12 : char **netCDFDataset::GetMetadata( const char *pszDomain )
76 : {
77 12 : if( pszDomain != NULL && EQUALN( pszDomain, "SUBDATASETS", 11 ) )
78 0 : return papszSubDatasets;
79 : else
80 12 : return GDALDataset::GetMetadata( pszDomain );
81 : }
82 :
83 :
84 : /************************************************************************/
85 : /* GetProjectionRef() */
86 : /************************************************************************/
87 :
88 25 : const char * netCDFDataset::GetProjectionRef()
89 : {
90 25 : if( bGotGeoTransform )
91 17 : return pszProjection;
92 : else
93 8 : return GDALPamDataset::GetProjectionRef();
94 : }
95 :
96 : /************************************************************************/
97 : /* GetNoDataValue() */
98 : /************************************************************************/
99 :
100 1 : double netCDFRasterBand::GetNoDataValue( int * pbSuccess )
101 :
102 : {
103 1 : if( pbSuccess )
104 1 : *pbSuccess = bNoDataSet;
105 :
106 1 : if( bNoDataSet )
107 1 : return dfNoDataValue;
108 : else
109 0 : return GDALPamRasterBand::GetNoDataValue( pbSuccess );
110 : }
111 :
112 : /************************************************************************/
113 : /* SetNoDataValue() */
114 : /************************************************************************/
115 :
116 27 : CPLErr netCDFRasterBand::SetNoDataValue( double dfNoData )
117 :
118 : {
119 27 : bNoDataSet = TRUE;
120 27 : dfNoDataValue = dfNoData;
121 :
122 27 : return CE_None;
123 : }
124 :
125 : /************************************************************************/
126 : /* ~netCDFRasterBand() */
127 : /************************************************************************/
128 :
129 54 : netCDFRasterBand::~netCDFRasterBand()
130 : {
131 27 : if( panBandZPos )
132 16 : CPLFree( panBandZPos );
133 27 : if( panBandZLev )
134 16 : CPLFree( panBandZLev );
135 54 : }
136 :
137 : /************************************************************************/
138 : /* CreateBandMetadata() */
139 : /************************************************************************/
140 :
141 27 : CPLErr netCDFRasterBand::CreateBandMetadata( )
142 : {
143 : char szVarName[NC_MAX_NAME];
144 : char szMetaName[NC_MAX_NAME];
145 : char szMetaTemp[8192];
146 : int nd;
147 : int i,j;
148 27 : int Sum = 1;
149 27 : int Taken = 0;
150 27 : int result = 0;
151 : int status;
152 27 : int nVarID = -1;
153 : int nDims;
154 : size_t start[1];
155 : size_t count[1];
156 : char szTemp[NC_MAX_NAME];
157 : const char *pszValue;
158 :
159 : nc_type nVarType;
160 : netCDFDataset *poDS;
161 :
162 27 : poDS = (netCDFDataset *) this->poDS;
163 : /* -------------------------------------------------------------------- */
164 : /* Compute all dimensions from Band number and save in Metadata */
165 : /* -------------------------------------------------------------------- */
166 27 : nc_inq_varname( poDS->cdfid, nZId, szVarName );
167 27 : nc_inq_varndims( poDS->cdfid, nZId, &nd );
168 : /* -------------------------------------------------------------------- */
169 : /* Compute multidimention band position */
170 : /* */
171 : /* BandPosition = (Total - sum(PastBandLevels) - 1)/sum(remainingLevels)*/
172 : /* if Data[2,3,4,x,y] */
173 : /* */
174 : /* BandPos0 = (nBand ) / (3*4) */
175 : /* BandPos1 = (nBand - BandPos0*(3*4) ) / (4) */
176 : /* BandPos2 = (nBand - BandPos0*(3*4) ) % (4) */
177 : /* -------------------------------------------------------------------- */
178 :
179 27 : sprintf( szMetaName,"NETCDF_VARNAME");
180 27 : sprintf( szMetaTemp,"%s",szVarName);
181 27 : SetMetadataItem( szMetaName, szMetaTemp );
182 27 : if( nd == 3 ) {
183 0 : Sum *= panBandZLev[0];
184 : }
185 :
186 75 : for( i=0; i < nd-2 ; i++ ) {
187 48 : if( i != nd - 2 -1 ) {
188 32 : Sum = 1;
189 80 : for( j=i+1; j < nd-2; j++ ) {
190 48 : Sum *= panBandZLev[j];
191 : }
192 32 : result = (int) ( ( nLevel-Taken) / Sum );
193 : }
194 : else {
195 16 : result = (int) ( ( nLevel-Taken) % Sum );
196 : }
197 :
198 : strcpy(szVarName, poDS->papszDimName[poDS->paDimIds[
199 48 : panBandZPos[i]]] );
200 :
201 48 : sprintf( szMetaName,"NETCDF_DIMENSION_%s", szVarName );
202 :
203 : status=nc_inq_varid(poDS->cdfid,
204 : szVarName,
205 48 : &nVarID );
206 :
207 : /* -------------------------------------------------------------------- */
208 : /* Try to uppercase the first letter of the variable */
209 : /* -------------------------------------------------------------------- */
210 :
211 48 : if( status != NC_NOERR ) {
212 48 : szVarName[0]=toupper(szVarName[0]);
213 : status=nc_inq_varid(poDS->cdfid,
214 : szVarName,
215 48 : &nVarID );
216 : }
217 :
218 48 : status = nc_inq_vartype( poDS->cdfid, nVarID, &nVarType );
219 :
220 48 : nDims = 0;
221 48 : status = nc_inq_varndims( poDS->cdfid, nVarID, &nDims );
222 :
223 48 : if( nDims == 1 ) {
224 0 : count[0]=1;
225 0 : start[0]=result;
226 0 : switch( nVarType ) {
227 : case NC_SHORT:
228 : short sData;
229 : status = nc_get_vara_short( poDS->cdfid, nVarID,
230 : start,
231 0 : count, &sData );
232 0 : sprintf( szMetaTemp,"%d", sData );
233 0 : break;
234 : case NC_INT:
235 : int nData;
236 : status = nc_get_vara_int( poDS->cdfid, nVarID,
237 : start,
238 0 : count, &nData );
239 0 : sprintf( szMetaTemp,"%d", nData );
240 0 : break;
241 : case NC_FLOAT:
242 : float fData;
243 : status = nc_get_vara_float( poDS->cdfid, nVarID,
244 : start,
245 0 : count, &fData );
246 0 : sprintf( szMetaTemp,"%f", fData );
247 0 : break;
248 : case NC_DOUBLE:
249 : double dfData;
250 : status = nc_get_vara_double( poDS->cdfid, nVarID,
251 : start,
252 0 : count, &dfData);
253 0 : sprintf( szMetaTemp,"%g", dfData );
254 : break;
255 : default:
256 : break;
257 : }
258 : }
259 : else
260 48 : sprintf( szMetaTemp,"%d", result+1);
261 :
262 48 : SetMetadataItem( szMetaName, szMetaTemp );
263 :
264 : /* -------------------------------------------------------------------- */
265 : /* Fetch dimension units */
266 : /* -------------------------------------------------------------------- */
267 :
268 48 : strcpy( szTemp, szVarName );
269 48 : strcat( szTemp, "#units" );
270 48 : pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp);
271 48 : if( pszValue != NULL ) {
272 0 : if( EQUAL( pszValue, "T") ) {
273 0 : strcpy( szTemp, szVarName );
274 0 : strcat( szTemp, "#original_units" );
275 0 : pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp);
276 0 : strcpy( szTemp, "NETCDF_");
277 0 : strcat( szTemp, szVarName );
278 0 : strcat( szTemp, "_original_units" );
279 0 : SetMetadataItem( szTemp, pszValue );
280 : }
281 : else {
282 0 : strcpy( szTemp, "NETCDF_");
283 0 : strcat( szTemp, szVarName );
284 0 : strcat( szTemp, "_units" );
285 0 : SetMetadataItem( szTemp, pszValue );
286 : }
287 : }
288 48 : Taken += result * Sum;
289 : }
290 :
291 :
292 :
293 27 : return CE_None;
294 : }
295 :
296 : /************************************************************************/
297 : /* netCDFRasterBand() */
298 : /************************************************************************/
299 :
300 27 : netCDFRasterBand::netCDFRasterBand( netCDFDataset *poDS,
301 : int nZId,
302 : int nZDim,
303 : int nLevel,
304 : int *panBandZLev,
305 : int *panBandDimPos,
306 27 : int nBand)
307 :
308 : {
309 : double dfNoData;
310 27 : int bNoDataSet = FALSE;
311 27 : nc_type vartype=NC_NAT;
312 27 : nc_type atttype=NC_NAT;
313 : size_t attlen;
314 : int status;
315 : char szNoValueName[8192];
316 :
317 :
318 27 : this->panBandZPos = NULL;
319 27 : this->panBandZLev = NULL;
320 27 : this->poDS = poDS;
321 27 : this->nBand = nBand;
322 27 : this->nZId = nZId;
323 27 : this->nZDim = nZDim;
324 27 : this->nLevel = nLevel;
325 27 : this->nBandXPos = panBandDimPos[0];
326 27 : this->nBandYPos = panBandDimPos[1];
327 :
328 : /* -------------------------------------------------------------------- */
329 : /* Take care of all other dimmensions */
330 : /* ------------------------------------------------------------------ */
331 27 : if( nZDim > 2 ) {
332 : this->panBandZPos =
333 16 : (int *) CPLCalloc( nZDim-1, sizeof( int ) );
334 : this->panBandZLev =
335 16 : (int *) CPLCalloc( nZDim-1, sizeof( int ) );
336 :
337 64 : for ( int i=0; i < nZDim - 2; i++ ){
338 48 : this->panBandZPos[i] = panBandDimPos[i+2];
339 48 : this->panBandZLev[i] = panBandZLev[i];
340 : }
341 : }
342 27 : CreateBandMetadata();
343 27 : bNoDataSet = FALSE;
344 27 : dfNoDataValue = -9999.0;
345 :
346 27 : nBlockXSize = poDS->GetRasterXSize( );
347 27 : nBlockYSize = 1;
348 :
349 : /* -------------------------------------------------------------------- */
350 : /* Get the type of the "z" variable, our target raster array. */
351 : /* -------------------------------------------------------------------- */
352 27 : if( nc_inq_var( poDS->cdfid, nZId, NULL, &nc_datatype, NULL, NULL,
353 : NULL ) != NC_NOERR ){
354 : CPLError( CE_Failure, CPLE_AppDefined,
355 0 : "Error in nc_var_inq() on 'z'." );
356 0 : return;
357 : }
358 :
359 27 : if( (nc_datatype == NC_BYTE) )
360 3 : eDataType = GDT_Byte;
361 24 : else if( nc_datatype == NC_SHORT )
362 2 : eDataType = GDT_Int16;
363 22 : else if( nc_datatype == NC_INT )
364 2 : eDataType = GDT_Int32;
365 20 : else if( nc_datatype == NC_FLOAT )
366 3 : eDataType = GDT_Float32;
367 17 : else if( nc_datatype == NC_DOUBLE )
368 17 : eDataType = GDT_Float64;
369 : else
370 : {
371 0 : if( nBand == 1 )
372 : CPLError( CE_Warning, CPLE_AppDefined,
373 : "Unsupported netCDF datatype (%d), treat as Float32.",
374 0 : (int) nc_datatype );
375 0 : eDataType = GDT_Float32;
376 : }
377 : /* -------------------------------------------------------------------- */
378 : /* Find out what is No Data for this variable */
379 : /* -------------------------------------------------------------------- */
380 :
381 : status = nc_inq_att( poDS->cdfid, nZId,
382 27 : _FillValue, &atttype, &attlen);
383 :
384 : /* -------------------------------------------------------------------- */
385 : /* Look for either Missing_Value or _FillValue attributes */
386 : /* -------------------------------------------------------------------- */
387 :
388 27 : if( status == NC_NOERR ) {
389 7 : strcpy(szNoValueName, _FillValue );
390 : }
391 : else {
392 : status = nc_inq_att( poDS->cdfid, nZId,
393 20 : "missing_value", &atttype, &attlen );
394 20 : if( status == NC_NOERR ) {
395 :
396 1 : strcpy( szNoValueName, "missing_value" );
397 : }
398 : }
399 :
400 27 : nc_inq_vartype( poDS->cdfid, nZId, &vartype );
401 :
402 27 : if( status == NC_NOERR ) {
403 8 : switch( atttype ) {
404 : case NC_CHAR:
405 : char *fillc;
406 0 : fillc = (char *) CPLCalloc( attlen+1, sizeof(char) );
407 : status=nc_get_att_text( poDS->cdfid, nZId,
408 0 : szNoValueName, fillc );
409 0 : dfNoData = atof( fillc );
410 0 : CPLFree(fillc);
411 0 : break;
412 : case NC_SHORT:
413 : short sNoData;
414 : status = nc_get_att_short( poDS->cdfid, nZId,
415 2 : szNoValueName, &sNoData );
416 2 : dfNoData = (double) sNoData;
417 2 : break;
418 : case NC_INT:
419 : int nNoData;
420 : status = nc_get_att_int( poDS->cdfid, nZId,
421 2 : szNoValueName, &nNoData );
422 2 : dfNoData = (double) nNoData;
423 2 : break;
424 : case NC_FLOAT:
425 : float fNoData;
426 : status = nc_get_att_float( poDS->cdfid, nZId,
427 2 : szNoValueName, &fNoData );
428 2 : dfNoData = (double) fNoData;
429 2 : break;
430 : case NC_DOUBLE:
431 : status = nc_get_att_double( poDS->cdfid, nZId,
432 2 : szNoValueName, &dfNoData );
433 : break;
434 : default:
435 : break;
436 : }
437 : status = nc_get_att_double( poDS->cdfid, nZId,
438 8 : szNoValueName, &dfNoData );
439 :
440 : } else {
441 19 : switch( vartype ) {
442 : case NC_BYTE:
443 : /* don't do default fill-values for bytes, too risky */
444 3 : dfNoData = 0.0;
445 3 : break;
446 : case NC_CHAR:
447 0 : dfNoData = NC_FILL_CHAR;
448 0 : break;
449 : case NC_SHORT:
450 0 : dfNoData = NC_FILL_SHORT;
451 0 : break;
452 : case NC_INT:
453 0 : dfNoData = NC_FILL_INT;
454 0 : break;
455 : case NC_FLOAT:
456 0 : dfNoData = NC_FILL_FLOAT;
457 0 : break;
458 : case NC_DOUBLE:
459 16 : dfNoData = NC_FILL_DOUBLE;
460 16 : break;
461 : default:
462 0 : dfNoData = 0.0;
463 : break;
464 : }
465 19 : bNoDataSet = TRUE;
466 : }
467 27 : SetNoDataValue( dfNoData );
468 :
469 :
470 0 : }
471 :
472 : /************************************************************************/
473 : /* IReadBlock() */
474 : /************************************************************************/
475 :
476 205 : CPLErr netCDFRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
477 : void * pImage )
478 :
479 : {
480 205 : int nErr=-1;
481 205 : int cdfid = ( ( netCDFDataset * ) poDS )->cdfid;
482 : size_t start[ MAX_NC_DIMS ];
483 : size_t edge[ MAX_NC_DIMS ];
484 : char pszName[ MAX_STR_LEN ];
485 : int i,j;
486 205 : int Sum=-1;
487 205 : int Taken=-1;
488 : int nd;
489 :
490 205 : *pszName='\0';
491 205 : memset( start, 0, sizeof( start ) );
492 205 : memset( edge, 0, sizeof( edge ) );
493 205 : nc_inq_varndims ( cdfid, nZId, &nd );
494 :
495 : /* -------------------------------------------------------------------- */
496 : /* Locate X, Y and Z position in the array */
497 : /* -------------------------------------------------------------------- */
498 :
499 205 : start[nBandXPos] = 0; // x dim can move arround in array
500 : // check y order
501 205 : if( ( ( netCDFDataset *) poDS )->bBottomUp ) {
502 0 : start[nBandYPos] = ( ( netCDFDataset * ) poDS )->ydim - 1 - nBlockYOff;
503 : } else {
504 205 : start[nBandYPos] = nBlockYOff; // y
505 : }
506 :
507 205 : edge[nBandXPos] = nBlockXSize;
508 205 : edge[nBandYPos] = 1;
509 :
510 205 : if( nd == 3 ) {
511 0 : start[panBandZPos[0]] = nLevel; // z
512 0 : edge [panBandZPos[0]] = 1;
513 : }
514 :
515 : /* -------------------------------------------------------------------- */
516 : /* Compute multidimention band position */
517 : /* */
518 : /* BandPosition = (Total - sum(PastBandLevels) - 1)/sum(remainingLevels)*/
519 : /* if Data[2,3,4,x,y] */
520 : /* */
521 : /* BandPos0 = (nBand ) / (3*4) */
522 : /* BandPos1 = (nBand - (3*4) ) / (4) */
523 : /* BandPos2 = (nBand - (3*4) ) % (4) */
524 : /* -------------------------------------------------------------------- */
525 205 : if (nd > 3)
526 : {
527 20 : Taken = 0;
528 80 : for( i=0; i < nd-2 ; i++ )
529 : {
530 60 : if( i != nd - 2 -1 ) {
531 40 : Sum = 1;
532 100 : for( j=i+1; j < nd-2; j++ ) {
533 60 : Sum *= panBandZLev[j];
534 : }
535 40 : start[panBandZPos[i]] = (int) ( ( nLevel-Taken) / Sum );
536 40 : edge[panBandZPos[i]] = 1;
537 : } else {
538 20 : start[panBandZPos[i]] = (int) ( ( nLevel-Taken) % Sum );
539 20 : edge[panBandZPos[i]] = 1;
540 : }
541 60 : Taken += start[panBandZPos[i]] * Sum;
542 : }
543 : }
544 :
545 205 : if( eDataType == GDT_Byte )
546 : nErr = nc_get_vara_uchar( cdfid, nZId, start, edge,
547 20 : (unsigned char *) pImage );
548 185 : else if( eDataType == GDT_Int16 )
549 : nErr = nc_get_vara_short( cdfid, nZId, start, edge,
550 0 : (short int *) pImage );
551 185 : else if( eDataType == GDT_Int32 )
552 : {
553 : if( sizeof(long) == 4 )
554 : nErr = nc_get_vara_long( cdfid, nZId, start, edge,
555 0 : (long *) pImage );
556 : else
557 : nErr = nc_get_vara_int( cdfid, nZId, start, edge,
558 : (int *) pImage );
559 : }
560 185 : else if( eDataType == GDT_Float32 ){
561 : nErr = nc_get_vara_float( cdfid, nZId, start, edge,
562 165 : (float *) pImage );
563 18558 : for( i=0; i<nBlockXSize; i++ ){
564 18393 : if( CPLIsNan( ( (float *) pImage )[i] ) )
565 0 : ( (float *)pImage )[i] = dfNoDataValue;
566 : }
567 : }
568 20 : else if( eDataType == GDT_Float64 ){
569 : nErr = nc_get_vara_double( cdfid, nZId, start, edge,
570 20 : (double *) pImage );
571 220 : for( i=0; i<nBlockXSize; i++ ){
572 200 : if( CPLIsNan( ( (double *) pImage)[i] ) )
573 0 : ( (double *)pImage )[i] = dfNoDataValue;
574 : }
575 :
576 : }
577 :
578 205 : if( nErr != NC_NOERR )
579 : {
580 : CPLError( CE_Failure, CPLE_AppDefined,
581 : "netCDF scanline fetch failed: %s",
582 0 : nc_strerror( nErr ) );
583 0 : return CE_Failure;
584 : }
585 : else
586 205 : return CE_None;
587 : }
588 :
589 : /************************************************************************/
590 : /* ==================================================================== */
591 : /* netCDFDataset */
592 : /* ==================================================================== */
593 : /************************************************************************/
594 :
595 : /************************************************************************/
596 : /* netCDFDataset() */
597 : /************************************************************************/
598 :
599 21 : netCDFDataset::netCDFDataset()
600 :
601 : {
602 21 : papszMetadata = NULL;
603 21 : papszSubDatasets = NULL;
604 21 : bGotGeoTransform = FALSE;
605 21 : pszProjection = NULL;
606 21 : cdfid = 0;
607 21 : bBottomUp = FALSE;
608 21 : }
609 :
610 :
611 : /************************************************************************/
612 : /* ~netCDFDataset() */
613 : /************************************************************************/
614 :
615 42 : netCDFDataset::~netCDFDataset()
616 :
617 : {
618 :
619 21 : FlushCache();
620 :
621 21 : CSLDestroy( papszMetadata );
622 21 : CSLDestroy( papszSubDatasets );
623 :
624 21 : CPLFree( pszProjection );
625 :
626 21 : if( cdfid )
627 21 : nc_close( cdfid );
628 42 : }
629 :
630 : /************************************************************************/
631 : /* FetchCopyParm() */
632 : /************************************************************************/
633 :
634 10 : double netCDFDataset::FetchCopyParm( const char *pszGridMappingValue,
635 : const char *pszParm, double dfDefault )
636 :
637 : {
638 : char szTemp[ MAX_NC_NAME ];
639 : const char *pszValue;
640 :
641 10 : strcpy(szTemp,pszGridMappingValue);
642 10 : strcat( szTemp, "#" );
643 10 : strcat( szTemp, pszParm );
644 10 : pszValue = CSLFetchNameValue(papszMetadata, szTemp);
645 :
646 10 : if( pszValue )
647 : {
648 10 : return CPLAtofM(pszValue);
649 : }
650 : else
651 0 : return dfDefault;
652 : }
653 :
654 :
655 : /************************************************************************/
656 : /* SetProjection() */
657 : /************************************************************************/
658 13 : void netCDFDataset::SetProjection( int var )
659 : {
660 : /* -------------------------------------------------------------------- */
661 : /* Set Projection */
662 : /* -------------------------------------------------------------------- */
663 :
664 : size_t start[2], edge[2];
665 : int status;
666 : unsigned int i;
667 : const char *pszValue;
668 : int nVarProjectionID;
669 : char szVarName[ MAX_NC_NAME ];
670 : char szTemp[ MAX_NC_NAME ];
671 : char szGridMappingName[ MAX_NC_NAME ];
672 : char szGridMappingValue[ MAX_NC_NAME ];
673 :
674 : double dfStdP1;
675 : double dfStdP2;
676 : double dfCenterLat;
677 : double dfCenterLon;
678 : double dfScale;
679 : double dfFalseEasting;
680 : double dfFalseNorthing;
681 : double dfCentralMeridian;
682 :
683 13 : OGRSpatialReference oSRS;
684 13 : int nVarDimXID = -1;
685 13 : int nVarDimYID = -1;
686 : double *pdfXCoord;
687 : double *pdfYCoord;
688 : char szDimNameX[ MAX_NC_NAME ];
689 : char szDimNameY[ MAX_NC_NAME ];
690 : int nSpacingBegin;
691 : int nSpacingMiddle;
692 : int nSpacingLast;
693 :
694 : const char *pszWKT;
695 : const char *pszGeoTransform;
696 13 : char **papszGeoTransform=NULL;
697 :
698 :
699 : netCDFDataset * poDS;
700 13 : poDS = this;
701 :
702 : /* -------------------------------------------------------------------- */
703 : /* Get x/y range information. */
704 : /* -------------------------------------------------------------------- */
705 :
706 13 : poDS->adfGeoTransform[0] = 0.0;
707 13 : poDS->adfGeoTransform[1] = 1.0;
708 13 : poDS->adfGeoTransform[2] = 0.0;
709 13 : poDS->adfGeoTransform[3] = 0.0;
710 13 : poDS->adfGeoTransform[4] = 0.0;
711 13 : poDS->adfGeoTransform[5] = 1.0;
712 13 : poDS->pszProjection = NULL;
713 :
714 :
715 : /* -------------------------------------------------------------------- */
716 : /* Look for grid_mapping metadata */
717 : /* -------------------------------------------------------------------- */
718 :
719 13 : strcpy( szGridMappingValue, "" );
720 13 : strcpy( szGridMappingName, "" );
721 :
722 13 : nc_inq_varname( cdfid, var, szVarName );
723 13 : strcpy(szTemp,szVarName);
724 13 : strcat(szTemp,"#grid_mapping");
725 13 : pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp);
726 13 : if( pszValue ) {
727 9 : strcpy(szGridMappingName,szTemp);
728 9 : strcpy(szGridMappingValue,pszValue);
729 : }
730 :
731 : /* -------------------------------------------------------------------- */
732 : /* Look for dimension: lon */
733 : /* -------------------------------------------------------------------- */
734 :
735 32 : for( i = 0; (i < strlen( poDS->papszDimName[ poDS->nDimXid ] ) &&
736 : i < 3 ); i++ ) {
737 19 : szDimNameX[i] = tolower( ( poDS->papszDimName[poDS->nDimXid] )[i] );
738 : }
739 13 : szDimNameX[3] = '\0';
740 32 : for( i = 0; (i < strlen( poDS->papszDimName[ poDS->nDimYid ] ) &&
741 : i < 3 ); i++ ) {
742 19 : szDimNameY[i] = tolower( ( poDS->papszDimName[poDS->nDimYid] )[i] );
743 : }
744 13 : szDimNameY[3] = '\0';
745 :
746 : /* -------------------------------------------------------------------- */
747 : /* Read grid_mappinginformation and set projections */
748 : /* -------------------------------------------------------------------- */
749 :
750 13 : if( !( EQUAL(szGridMappingName,"" ) ) ) {
751 9 : nc_inq_varid( cdfid, szGridMappingValue, &nVarProjectionID );
752 9 : poDS->ReadAttributes( cdfid, nVarProjectionID );
753 :
754 9 : strcpy( szTemp, szGridMappingValue );
755 9 : strcat( szTemp, "#" );
756 9 : strcat( szTemp, GRD_MAPPING_NAME );
757 9 : pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp);
758 :
759 9 : if( pszValue != NULL ) {
760 : /* -------------------------------------------------------------------- */
761 : /* Transverse Mercator */
762 : /* -------------------------------------------------------------------- */
763 :
764 9 : if( EQUAL( pszValue, TM ) ) {
765 :
766 : dfScale =
767 : poDS->FetchCopyParm( szGridMappingValue,
768 2 : SCALE_FACTOR, 1.0 );
769 :
770 : dfCenterLon =
771 : poDS->FetchCopyParm( szGridMappingValue,
772 2 : LONG_CENTRAL_MERIDIAN, 0.0 );
773 :
774 : dfCenterLat =
775 : poDS->FetchCopyParm( szGridMappingValue,
776 2 : LAT_PROJ_ORIGIN, 0.0 );
777 :
778 : dfFalseEasting =
779 : poDS->FetchCopyParm( szGridMappingValue,
780 2 : FALSE_EASTING, 0.0 );
781 :
782 : dfFalseNorthing =
783 : poDS->FetchCopyParm( szGridMappingValue,
784 2 : FALSE_NORTHING, 0.0 );
785 :
786 : oSRS.SetTM( dfCenterLat,
787 : dfCenterLon,
788 : dfScale,
789 : dfFalseEasting,
790 2 : dfFalseNorthing );
791 2 : oSRS.SetWellKnownGeogCS( "WGS84" );
792 : }
793 :
794 : /* -------------------------------------------------------------------- */
795 : /* Cylindrical Equal Area */
796 : /* -------------------------------------------------------------------- */
797 :
798 7 : else if( EQUAL( pszValue, CEA ) ) {
799 : dfStdP1 =
800 : poDS->FetchCopyParm( szGridMappingValue,
801 0 : STD_PARALLEL_1, 0.0 );
802 : dfCentralMeridian =
803 : poDS->FetchCopyParm( szGridMappingValue,
804 0 : LONG_CENTRAL_MERIDIAN, 0.0 );
805 :
806 : dfFalseEasting =
807 : poDS->FetchCopyParm( szGridMappingValue,
808 0 : FALSE_EASTING, 0.0 );
809 :
810 : dfFalseNorthing =
811 : poDS->FetchCopyParm( szGridMappingValue,
812 0 : FALSE_NORTHING, 0.0 );
813 :
814 : oSRS.SetCEA( dfStdP1, dfCentralMeridian,
815 0 : dfFalseEasting, dfFalseNorthing );
816 :
817 0 : oSRS.SetWellKnownGeogCS( "WGS84" );
818 : }
819 : /* -------------------------------------------------------------------- */
820 : /* lambert_azimuthal_equal_area */
821 : /* -------------------------------------------------------------------- */
822 7 : else if( EQUAL( pszValue, LAEA ) ) {
823 : dfCenterLon =
824 : poDS->FetchCopyParm( szGridMappingValue,
825 0 : LON_PROJ_ORIGIN, 0.0 );
826 :
827 : dfCenterLat =
828 : poDS->FetchCopyParm( szGridMappingValue,
829 0 : LAT_PROJ_ORIGIN, 0.0 );
830 :
831 : dfFalseEasting =
832 : poDS->FetchCopyParm( szGridMappingValue,
833 0 : FALSE_EASTING, 0.0 );
834 :
835 : dfFalseNorthing =
836 : poDS->FetchCopyParm( szGridMappingValue,
837 0 : FALSE_NORTHING, 0.0 );
838 :
839 : /*
840 : dfLonOrig =
841 : poDS->FetchCopyParm( szGridMappingValue,
842 : LON_PROJ_ORIGIN, 0.0 );
843 :
844 : dfLatOrig =
845 : poDS->FetchCopyParm( szGridMappingValue,
846 : LAT_PROJ_ORIGIN, 0.0 );
847 :
848 : dfScaleFactorOrig =
849 : poDS->FetchCopyParm( szGridMappingValue,
850 : SCALE_FACTOR_ORIGN, 0.0 );
851 :
852 : dfProjXOrig =
853 : poDS->FetchCopyParm( szGridMappingValue,
854 : PROJ_X_ORIGIN, 0.0 );
855 :
856 : dfProjYOrig =
857 : poDS->FetchCopyParm( szGridMappingValue,
858 : PROJ_Y_ORIGIN, 0.0 );
859 :
860 : dfFalseEasting =
861 : poDS->FetchCopyParm( szGridMappingValue,
862 : FALSE_EASTING, 0.0 );
863 :
864 : dfFalseNorthing =
865 : poDS->FetchCopyParm( szGridMappingValue,
866 : FALSE_NORTHING, 0.0 );
867 : */
868 0 : oSRS.SetProjCS( "LAEA (WGS84) " );
869 0 : oSRS.SetWellKnownGeogCS( "WGS84" );
870 :
871 : oSRS.SetLAEA( dfCenterLat, dfCenterLon,
872 0 : dfFalseEasting, dfFalseNorthing );
873 :
874 :
875 : }
876 :
877 :
878 : /* -------------------------------------------------------------------- */
879 : /* Lambert conformal conic */
880 : /* -------------------------------------------------------------------- */
881 7 : else if( EQUAL( pszValue, L_C_CONIC ) ) {
882 :
883 : dfStdP1 =
884 : poDS->FetchCopyParm( szGridMappingValue,
885 0 : STD_PARALLEL_1, 0.0 );
886 :
887 : dfStdP2 =
888 : poDS->FetchCopyParm( szGridMappingValue,
889 0 : STD_PARALLEL_2, 0.0 );
890 :
891 : dfCenterLon =
892 : poDS->FetchCopyParm( szGridMappingValue,
893 0 : LONG_CENTRAL_MERIDIAN, 0.0 );
894 :
895 : dfCenterLat =
896 : poDS->FetchCopyParm( szGridMappingValue,
897 0 : LAT_PROJ_ORIGIN, 0.0 );
898 :
899 : dfFalseEasting =
900 : poDS->FetchCopyParm( szGridMappingValue,
901 0 : FALSE_EASTING, 0.0 );
902 :
903 : dfFalseNorthing =
904 : poDS->FetchCopyParm( szGridMappingValue,
905 0 : FALSE_NORTHING, 0.0 );
906 :
907 0 : oSRS.SetWellKnownGeogCS( "WGS84" );
908 : oSRS.SetLCC( dfStdP1, dfStdP2, dfCenterLat, dfCenterLon,
909 0 : dfFalseEasting, dfFalseNorthing );
910 :
911 : }
912 : /* -------------------------------------------------------------------- */
913 : /* Is this Latitude/Longitude Grid */
914 : /* -------------------------------------------------------------------- */
915 :
916 0 : } else if( EQUAL( szDimNameX,"lon" ) ) {
917 0 : oSRS.SetWellKnownGeogCS( "WGS84" );
918 :
919 : } else {
920 : // This would be too indiscrimant. But we should set
921 : // it if we know the data is geographic.
922 : //oSRS.SetWellKnownGeogCS( "WGS84" );
923 : }
924 : }
925 : /* -------------------------------------------------------------------- */
926 : /* Read projection coordinates */
927 : /* -------------------------------------------------------------------- */
928 :
929 13 : nc_inq_varid( cdfid, poDS->papszDimName[nDimXid], &nVarDimXID );
930 13 : nc_inq_varid( cdfid, poDS->papszDimName[nDimYid], &nVarDimYID );
931 :
932 13 : if( ( nVarDimXID != -1 ) && ( nVarDimYID != -1 ) ) {
933 2 : pdfXCoord = (double *) CPLCalloc( xdim, sizeof(double) );
934 2 : pdfYCoord = (double *) CPLCalloc( ydim, sizeof(double) );
935 :
936 : /* -------------------------------------------------------------------- */
937 : /* Is pixel spacing is uniform accross the map? */
938 : /* -------------------------------------------------------------------- */
939 2 : start[0] = 0;
940 2 : edge[0] = xdim;
941 :
942 : status = nc_get_vara_double( cdfid, nVarDimXID,
943 2 : start, edge, pdfXCoord);
944 2 : edge[0] = ydim;
945 : status = nc_get_vara_double( cdfid, nVarDimYID,
946 2 : start, edge, pdfYCoord);
947 :
948 : /* -------------------------------------------------------------------- */
949 : /* Check Longitude */
950 : /* -------------------------------------------------------------------- */
951 :
952 2 : nSpacingBegin = (int) poDS->rint((pdfXCoord[1]-pdfXCoord[0]) * 1000);
953 :
954 2 : nSpacingMiddle = (int) poDS->rint((pdfXCoord[xdim / 2] -
955 2 : pdfXCoord[(xdim / 2) + 1]) * 1000);
956 :
957 2 : nSpacingLast = (int) poDS->rint((pdfXCoord[xdim - 2] -
958 2 : pdfXCoord[xdim-1]) * 1000);
959 :
960 2 : if( ( abs( nSpacingBegin ) == abs( nSpacingLast ) ) &&
961 : ( abs( nSpacingBegin ) == abs( nSpacingMiddle ) ) &&
962 : ( abs( nSpacingMiddle ) == abs( nSpacingLast ) ) ) {
963 :
964 : /* -------------------------------------------------------------------- */
965 : /* Longitude is equaly spaced, check lattitde */
966 : /* -------------------------------------------------------------------- */
967 1 : nSpacingBegin = (int) poDS->rint((pdfYCoord[1]-pdfYCoord[0]) *
968 1 : 1000);
969 :
970 1 : nSpacingMiddle = (int) poDS->rint((pdfYCoord[ydim / 2] -
971 1 : pdfYCoord[(ydim / 2) + 1]) *
972 1 : 1000);
973 :
974 1 : nSpacingLast = (int) poDS->rint((pdfYCoord[ydim - 2] -
975 1 : pdfYCoord[ydim-1]) *
976 1 : 1000);
977 :
978 :
979 : /* -------------------------------------------------------------------- */
980 : /* For Latitude we allow an error of 0.1 degrees for gaussion */
981 : /* gridding */
982 : /* -------------------------------------------------------------------- */
983 :
984 1 : if((( abs( abs(nSpacingBegin) - abs(nSpacingLast) ) ) < 100 ) &&
985 : (( abs( abs(nSpacingBegin) - abs(nSpacingMiddle) ) ) < 100 ) &&
986 : (( abs( abs(nSpacingMiddle) - abs(nSpacingLast) ) ) < 100) ) {
987 :
988 1 : if( ( abs( nSpacingBegin ) != abs( nSpacingLast ) ) ||
989 : ( abs( nSpacingBegin ) != abs( nSpacingMiddle ) ) ||
990 : ( abs( nSpacingMiddle ) != abs( nSpacingLast ) ) ) {
991 :
992 0 : CPLError(CE_Warning, 1,"Latitude grid not spaced evenly.\nSeting projection for grid spacing is within 0.1 degrees threshold.\n");
993 :
994 : }
995 : /* -------------------------------------------------------------------- */
996 : /* We have gridded data s we can set the Gereferencing info. */
997 : /* -------------------------------------------------------------------- */
998 :
999 : /* -------------------------------------------------------------------- */
1000 : /* Enable GeoTransform */
1001 : /* -------------------------------------------------------------------- */
1002 : /* ----------------------------------------------------------*/
1003 : /* In the following "actual_range" and "node_offset" */
1004 : /* are attributes used by netCDF files created by GMT. */
1005 : /* If we find them we know how to proceed. Else, use */
1006 : /* the original algorithm. */
1007 : /* --------------------------------------------------------- */
1008 : double dummy[2], xMinMax[2], yMinMax[2];
1009 1 : int node_offset = 0;
1010 :
1011 1 : poDS->bGotGeoTransform = TRUE;
1012 :
1013 1 : nc_get_att_int (cdfid, NC_GLOBAL, "node_offset", &node_offset);
1014 :
1015 1 : if (!nc_get_att_double (cdfid, nVarDimXID, "actual_range", dummy)) {
1016 1 : xMinMax[0] = dummy[0];
1017 1 : xMinMax[1] = dummy[1];
1018 : }
1019 : else {
1020 0 : xMinMax[0] = pdfXCoord[0];
1021 0 : xMinMax[1] = pdfXCoord[xdim-1];
1022 0 : node_offset = 0;
1023 : }
1024 :
1025 1 : if (!nc_get_att_double (cdfid, nVarDimYID, "actual_range", dummy)) {
1026 1 : yMinMax[0] = dummy[0];
1027 1 : yMinMax[1] = dummy[1];
1028 : }
1029 : else {
1030 0 : yMinMax[0] = pdfYCoord[0];
1031 0 : yMinMax[1] = pdfYCoord[ydim-1];
1032 0 : node_offset = 0;
1033 : }
1034 : /* for CF-1 conventions, assume bottom first */
1035 1 : if( EQUAL( szDimNameY, "lat" ) && pdfYCoord[0] < pdfYCoord[1] )
1036 0 : poDS->bBottomUp = TRUE;
1037 :
1038 : // Check for reverse order of y-coordinate
1039 1 : if ( yMinMax[0] > yMinMax[1] ) {
1040 0 : dummy[0] = yMinMax[1];
1041 0 : dummy[1] = yMinMax[0];
1042 0 : yMinMax[0] = dummy[0];
1043 0 : yMinMax[1] = dummy[1];
1044 : }
1045 :
1046 1 : poDS->adfGeoTransform[0] = xMinMax[0];
1047 1 : poDS->adfGeoTransform[2] = 0;
1048 1 : poDS->adfGeoTransform[3] = yMinMax[1];
1049 1 : poDS->adfGeoTransform[4] = 0;
1050 1 : poDS->adfGeoTransform[1] = ( xMinMax[1] - xMinMax[0] ) /
1051 1 : ( poDS->nRasterXSize + (node_offset - 1) );
1052 1 : poDS->adfGeoTransform[5] = ( yMinMax[0] - yMinMax[1] ) /
1053 1 : ( poDS->nRasterYSize + (node_offset - 1) );
1054 :
1055 : /* -------------------------------------------------------------------- */
1056 : /* Compute the center of the pixel */
1057 : /* -------------------------------------------------------------------- */
1058 1 : if ( !node_offset ) { // Otherwise its already the pixel center
1059 1 : poDS->adfGeoTransform[0] -= (poDS->adfGeoTransform[1] / 2);
1060 1 : poDS->adfGeoTransform[3] -= (poDS->adfGeoTransform[5] / 2);
1061 : }
1062 :
1063 1 : oSRS.exportToWkt( &(poDS->pszProjection) );
1064 :
1065 : }
1066 : }
1067 :
1068 2 : CPLFree( pdfXCoord );
1069 2 : CPLFree( pdfYCoord );
1070 : }
1071 :
1072 :
1073 : /* -------------------------------------------------------------------- */
1074 : /* Is this a netCDF file created by GDAL? */
1075 : /* -------------------------------------------------------------------- */
1076 13 : if( !EQUAL( szGridMappingValue, "" ) ) {
1077 9 : strcpy( szTemp,szGridMappingValue);
1078 9 : strcat( szTemp, "#" );
1079 9 : strcat( szTemp, "spatial_ref");
1080 9 : pszWKT = CSLFetchNameValue(poDS->papszMetadata, szTemp);
1081 :
1082 9 : if( pszWKT != NULL ) {
1083 :
1084 9 : pszProjection = CPLStrdup( pszWKT );
1085 :
1086 9 : strcpy(szTemp,szGridMappingValue);
1087 9 : strcat( szTemp, "#" );
1088 9 : strcat( szTemp, "GeoTransform");
1089 :
1090 9 : pszGeoTransform = CSLFetchNameValue(poDS->papszMetadata, szTemp);
1091 :
1092 : /* -------------------------------------------------------------------- */
1093 : /* Look for GeoTransform Array */
1094 : /* -------------------------------------------------------------------- */
1095 :
1096 9 : if( pszGeoTransform != NULL ) {
1097 : papszGeoTransform = CSLTokenizeString2( pszGeoTransform,
1098 : " ",
1099 9 : CSLT_HONOURSTRINGS );
1100 9 : poDS->bGotGeoTransform = TRUE;
1101 :
1102 9 : poDS->adfGeoTransform[0] = atof( papszGeoTransform[0] );
1103 9 : poDS->adfGeoTransform[1] = atof( papszGeoTransform[1] );
1104 9 : poDS->adfGeoTransform[2] = atof( papszGeoTransform[2] );
1105 9 : poDS->adfGeoTransform[3] = atof( papszGeoTransform[3] );
1106 9 : poDS->adfGeoTransform[4] = atof( papszGeoTransform[4] );
1107 9 : poDS->adfGeoTransform[5] = atof( papszGeoTransform[5] );
1108 : /* -------------------------------------------------------------------- */
1109 : /* Look for corner array values */
1110 : /* -------------------------------------------------------------------- */
1111 : } else {
1112 : double dfNN, dfSN, dfEE, dfWE;
1113 0 : strcpy(szTemp,szGridMappingValue);
1114 0 : strcat( szTemp, "#" );
1115 0 : strcat( szTemp, "Northernmost_Northing");
1116 0 : pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp);
1117 :
1118 0 : if( pszValue != NULL ) {
1119 0 : dfNN = atof( pszValue );
1120 : }
1121 0 : strcpy(szTemp,szGridMappingValue);
1122 0 : strcat( szTemp, "#" );
1123 0 : strcat( szTemp, "Southernmost_Northing");
1124 0 : pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp);
1125 :
1126 0 : if( pszValue != NULL ) {
1127 0 : dfSN = atof( pszValue );
1128 : }
1129 :
1130 0 : strcpy(szTemp,szGridMappingValue);
1131 0 : strcat( szTemp, "#" );
1132 0 : strcat( szTemp, "Easternmost_Easting");
1133 0 : pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp);
1134 :
1135 0 : if( pszValue != NULL ) {
1136 0 : dfEE = atof( pszValue );
1137 : }
1138 :
1139 0 : strcpy(szTemp,szGridMappingValue);
1140 0 : strcat( szTemp, "#" );
1141 0 : strcat( szTemp, "Westernmost_Easting");
1142 0 : pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp);
1143 :
1144 0 : if( pszValue != NULL ) {
1145 0 : dfWE = atof( pszValue );
1146 : }
1147 :
1148 0 : adfGeoTransform[0] = dfWE;
1149 : adfGeoTransform[1] = (dfEE - dfWE) /
1150 0 : ( poDS->GetRasterXSize() - 1 );
1151 0 : adfGeoTransform[2] = 0.0;
1152 0 : adfGeoTransform[3] = dfNN;
1153 0 : adfGeoTransform[4] = 0.0;
1154 : adfGeoTransform[5] = (dfSN - dfNN) /
1155 0 : ( poDS->GetRasterYSize() - 1 );
1156 : /* -------------------------------------------------------------------- */
1157 : /* Compute the center of the pixel */
1158 : /* -------------------------------------------------------------------- */
1159 : adfGeoTransform[0] = dfWE
1160 0 : - (adfGeoTransform[1] / 2);
1161 :
1162 : adfGeoTransform[3] = dfNN
1163 0 : - (adfGeoTransform[5] / 2);
1164 :
1165 :
1166 0 : bGotGeoTransform = TRUE;
1167 : }
1168 9 : CSLDestroy( papszGeoTransform );
1169 : }
1170 13 : }
1171 :
1172 13 : }
1173 : /************************************************************************/
1174 : /* GetGeoTransform() */
1175 : /************************************************************************/
1176 :
1177 12 : CPLErr netCDFDataset::GetGeoTransform( double * padfTransform )
1178 :
1179 : {
1180 12 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
1181 12 : if( bGotGeoTransform )
1182 8 : return CE_None;
1183 : else
1184 4 : return GDALPamDataset::GetGeoTransform( padfTransform );;
1185 : }
1186 :
1187 : /************************************************************************/
1188 : /* rint() */
1189 : /************************************************************************/
1190 :
1191 9 : double netCDFDataset::rint( double dfX)
1192 : {
1193 9 : if( dfX > 0 ) {
1194 3 : int nX = (int) (dfX+0.5);
1195 3 : if( nX % 2 ) {
1196 0 : double dfDiff = dfX - (double)nX;
1197 0 : if( dfDiff == -0.5 )
1198 0 : return double( nX-1 );
1199 : }
1200 3 : return double( nX );
1201 : } else {
1202 6 : int nX= (int) (dfX-0.5);
1203 6 : if( nX % 2 ) {
1204 1 : double dfDiff = dfX - (double)nX;
1205 1 : if( dfDiff == 0.5 )
1206 0 : return double(nX+1);
1207 : }
1208 6 : return double(nX);
1209 : }
1210 : }
1211 :
1212 : /************************************************************************/
1213 : /* ReadAttributes() */
1214 : /************************************************************************/
1215 67 : CPLErr netCDFDataset::SafeStrcat(char** ppszDest, char* pszSrc, size_t* nDestSize)
1216 : {
1217 : /* Reallocate the data string until the content fits */
1218 235 : while(*nDestSize < (strlen(*ppszDest) + strlen(pszSrc) + 1)) {
1219 101 : (*nDestSize) *= 2;
1220 101 : *ppszDest = (char*) CPLRealloc((void*) *ppszDest, *nDestSize);
1221 : }
1222 67 : strcat(*ppszDest, pszSrc);
1223 :
1224 67 : return CE_None;
1225 : }
1226 :
1227 49 : CPLErr netCDFDataset::ReadAttributes( int cdfid, int var)
1228 :
1229 : {
1230 : char szAttrName[ NC_MAX_NAME ];
1231 : char szVarName [ NC_MAX_NAME ];
1232 : char szMetaName[ NC_MAX_NAME * 2 ];
1233 49 : char *pszMetaTemp = NULL;
1234 : size_t nMetaTempSize;
1235 : nc_type nAttrType;
1236 : size_t nAttrLen, m;
1237 : int nbAttr;
1238 : char szTemp[ MAX_STR_LEN ];
1239 :
1240 49 : nc_inq_varnatts( cdfid, var, &nbAttr );
1241 49 : if( var == NC_GLOBAL ) {
1242 21 : strcpy( szVarName,"NC_GLOBAL" );
1243 : }
1244 : else {
1245 28 : nc_inq_varname( cdfid, var, szVarName );
1246 : }
1247 :
1248 216 : for( int l=0; l < nbAttr; l++) {
1249 :
1250 167 : nc_inq_attname( cdfid, var, l, szAttrName);
1251 167 : sprintf( szMetaName, "%s#%s", szVarName, szAttrName );
1252 167 : nc_inq_att( cdfid, var, szAttrName, &nAttrType, &nAttrLen );
1253 :
1254 : /* Allocate guaranteed minimum size */
1255 167 : nMetaTempSize = nAttrLen + 1;
1256 167 : pszMetaTemp = (char *) CPLCalloc( nMetaTempSize, sizeof( char ));
1257 167 : *pszMetaTemp = '\0';
1258 :
1259 167 : switch (nAttrType) {
1260 : case NC_CHAR:
1261 103 : nc_get_att_text( cdfid, var, szAttrName, pszMetaTemp );
1262 103 : pszMetaTemp[nAttrLen]='\0';
1263 103 : break;
1264 : case NC_SHORT:
1265 : short *psTemp;
1266 2 : psTemp = (short *) CPLCalloc( nAttrLen, sizeof( short ) );
1267 2 : nc_get_att_short( cdfid, var, szAttrName, psTemp );
1268 2 : for(m=0; m < nAttrLen-1; m++) {
1269 0 : sprintf( szTemp, "%hd, ", psTemp[m] );
1270 0 : SafeStrcat(&pszMetaTemp, szTemp, &nMetaTempSize);
1271 : }
1272 2 : sprintf( szTemp, "%hd", psTemp[m] );
1273 2 : SafeStrcat(&pszMetaTemp, szTemp, &nMetaTempSize);
1274 2 : CPLFree(psTemp);
1275 2 : break;
1276 : case NC_INT:
1277 : int *pnTemp;
1278 5 : pnTemp = (int *) CPLCalloc( nAttrLen, sizeof( int ) );
1279 5 : nc_get_att_int( cdfid, var, szAttrName, pnTemp );
1280 5 : for(m=0; m < nAttrLen-1; m++) {
1281 0 : sprintf( szTemp, "%d, ", pnTemp[m] );
1282 0 : SafeStrcat(&pszMetaTemp, szTemp, &nMetaTempSize);
1283 : }
1284 5 : sprintf( szTemp, "%d", pnTemp[m] );
1285 5 : SafeStrcat(&pszMetaTemp, szTemp, &nMetaTempSize);
1286 5 : CPLFree(pnTemp);
1287 5 : break;
1288 : case NC_FLOAT:
1289 : float *pfTemp;
1290 13 : pfTemp = (float *) CPLCalloc( nAttrLen, sizeof( float ) );
1291 13 : nc_get_att_float( cdfid, var, szAttrName, pfTemp );
1292 14 : for(m=0; m < nAttrLen-1; m++) {
1293 1 : sprintf( szTemp, "%e, ", pfTemp[m] );
1294 1 : SafeStrcat(&pszMetaTemp, szTemp, &nMetaTempSize);
1295 : }
1296 13 : sprintf( szTemp, "%e", pfTemp[m] );
1297 13 : SafeStrcat(&pszMetaTemp,szTemp, &nMetaTempSize);
1298 13 : CPLFree(pfTemp);
1299 13 : break;
1300 : case NC_DOUBLE:
1301 : double *pdfTemp;
1302 44 : pdfTemp = (double *) CPLCalloc(nAttrLen, sizeof(double));
1303 44 : nc_get_att_double( cdfid, var, szAttrName, pdfTemp );
1304 46 : for(m=0; m < nAttrLen-1; m++) {
1305 2 : sprintf( szTemp, "%g, ", pdfTemp[m] );
1306 2 : SafeStrcat(&pszMetaTemp, szTemp, &nMetaTempSize);
1307 : }
1308 44 : sprintf( szTemp, "%g", pdfTemp[m] );
1309 44 : SafeStrcat(&pszMetaTemp, szTemp, &nMetaTempSize);
1310 44 : CPLFree(pdfTemp);
1311 : break;
1312 : default:
1313 : break;
1314 : }
1315 :
1316 : papszMetadata = CSLSetNameValue(papszMetadata,
1317 : szMetaName,
1318 167 : pszMetaTemp);
1319 167 : CPLFree(pszMetaTemp);
1320 : }
1321 :
1322 :
1323 49 : return CE_None;
1324 :
1325 : }
1326 :
1327 :
1328 : /************************************************************************/
1329 : /* netCDFDataset::CreateSubDatasetList() */
1330 : /************************************************************************/
1331 4 : void netCDFDataset::CreateSubDatasetList( )
1332 : {
1333 :
1334 : char szDim[ MAX_NC_NAME ];
1335 : char szTemp[ MAX_NC_NAME ];
1336 : char szType[ MAX_NC_NAME ];
1337 : char szName[ MAX_NC_NAME ];
1338 : char szVarStdName[ MAX_NC_NAME ];
1339 : int nDims;
1340 : int nVar;
1341 : int nVarCount;
1342 : int i;
1343 : nc_type nVarType;
1344 : int *ponDimIds;
1345 : size_t nDimLen;
1346 : int nSub;
1347 : nc_type nAttype;
1348 : size_t nAttlen;
1349 :
1350 : netCDFDataset *poDS;
1351 4 : poDS = this;
1352 :
1353 4 : nSub=1;
1354 4 : nc_inq_nvars ( cdfid, &nVarCount );
1355 22 : for ( nVar = 0; nVar < nVarCount; nVar++ ) {
1356 :
1357 18 : nc_inq_varndims ( cdfid, nVar, &nDims );
1358 18 : if( nDims >= 2 ) {
1359 14 : ponDimIds = (int *) CPLCalloc( nDims, sizeof( int ) );
1360 14 : nc_inq_vardimid ( cdfid, nVar, ponDimIds );
1361 :
1362 : /* -------------------------------------------------------------------- */
1363 : /* Create Sub dataset list */
1364 : /* -------------------------------------------------------------------- */
1365 14 : szDim[0]='\0';
1366 42 : for( i = 0; i < nDims; i++ ) {
1367 28 : nc_inq_dimlen ( cdfid, ponDimIds[i], &nDimLen );
1368 28 : sprintf(szTemp, "%d", (int) nDimLen);
1369 28 : strcat(szTemp, "x" );
1370 28 : strcat(szDim, szTemp);
1371 : }
1372 :
1373 14 : nc_inq_vartype( cdfid, nVar, &nVarType );
1374 : /* -------------------------------------------------------------------- */
1375 : /* Get rid of the last "x" character */
1376 : /* -------------------------------------------------------------------- */
1377 14 : szDim[strlen(szDim) - 1] = '\0';
1378 14 : switch( nVarType ) {
1379 :
1380 : case NC_BYTE:
1381 : case NC_CHAR:
1382 14 : strcpy(szType, "8-bit character");
1383 14 : break;
1384 :
1385 : case NC_SHORT:
1386 0 : strcpy(szType, "8-bit integer");
1387 0 : break;
1388 : case NC_INT:
1389 0 : strcpy(szType, "16-bit integer");
1390 0 : break;
1391 : case NC_FLOAT:
1392 0 : strcpy(szType, "32-bit floating-point");
1393 0 : break;
1394 : case NC_DOUBLE:
1395 0 : strcpy(szType, "64-bit floating-point");
1396 : break;
1397 :
1398 : default:
1399 : break;
1400 : }
1401 14 : nc_inq_varname ( cdfid, nVar, szName);
1402 14 : nc_inq_att( cdfid, nVar, "standard_name", &nAttype, &nAttlen);
1403 14 : if( nc_get_att_text ( cdfid, nVar, "standard_name",
1404 : szVarStdName ) == NC_NOERR ) {
1405 0 : szVarStdName[nAttlen] = '\0';
1406 : }
1407 : else {
1408 14 : strcpy( szVarStdName, szName );
1409 : }
1410 :
1411 14 : sprintf( szTemp, "SUBDATASET_%d_NAME", nSub) ;
1412 :
1413 : poDS->papszSubDatasets =
1414 : CSLSetNameValue( poDS->papszSubDatasets, szTemp,
1415 : CPLSPrintf( "NETCDF:\"%s\":%s",
1416 : poDS->osFilename.c_str(),
1417 14 : szName) ) ;
1418 :
1419 14 : sprintf( szTemp, "SUBDATASET_%d_DESC", nSub++ );
1420 :
1421 : poDS->papszSubDatasets =
1422 : CSLSetNameValue( poDS->papszSubDatasets, szTemp,
1423 : CPLSPrintf( "[%s] %s (%s)",
1424 : szDim,
1425 : szVarStdName,
1426 14 : szType ) );
1427 14 : CPLFree(ponDimIds);
1428 : }
1429 : }
1430 :
1431 4 : }
1432 :
1433 : /************************************************************************/
1434 : /* Open() */
1435 : /************************************************************************/
1436 :
1437 9496 : GDALDataset *netCDFDataset::Open( GDALOpenInfo * poOpenInfo )
1438 :
1439 : {
1440 : int j;
1441 : unsigned int k;
1442 : int nd;
1443 : int cdfid, dim_count, var, var_count;
1444 9496 : int i = 0;
1445 : size_t lev_count;
1446 9496 : size_t nTotLevCount = 1;
1447 9496 : int nDim = 2;
1448 : int status;
1449 : int nDimID;
1450 : char attname[NC_MAX_NAME];
1451 : int ndims, nvars, ngatts, unlimdimid;
1452 9496 : int nCount=0;
1453 9496 : int nVarID=-1;
1454 :
1455 : /* -------------------------------------------------------------------- */
1456 : /* Does this appear to be a netcdf file? */
1457 : /* -------------------------------------------------------------------- */
1458 9496 : if( !EQUALN(poOpenInfo->pszFilename,"NETCDF:",7)
1459 : && ( poOpenInfo->nHeaderBytes < 5
1460 : || !EQUALN((const char *) (poOpenInfo->pabyHeader),"CDF\001",5)))
1461 9475 : return NULL;
1462 :
1463 : /* -------------------------------------------------------------------- */
1464 : /* Check if filename start with NETCDF: tag */
1465 : /* -------------------------------------------------------------------- */
1466 : netCDFDataset *poDS;
1467 21 : poDS = new netCDFDataset();
1468 21 : poDS->SetDescription( poOpenInfo->pszFilename );
1469 :
1470 21 : if( EQUALN( poOpenInfo->pszFilename,"NETCDF:",7) )
1471 : {
1472 : char **papszName =
1473 : CSLTokenizeString2( poOpenInfo->pszFilename,
1474 3 : ":", CSLT_HONOURSTRINGS|CSLT_PRESERVEESCAPES );
1475 :
1476 : /* -------------------------------------------------------------------- */
1477 : /* Check for drive name in windows NETCDF:"D:\... */
1478 : /* -------------------------------------------------------------------- */
1479 3 : if ( CSLCount(papszName) == 4 &&
1480 0 : strlen(papszName[1]) == 1 &&
1481 0 : (papszName[2][0] == '/' || papszName[2][0] == '\\') )
1482 : {
1483 0 : poDS->osFilename = papszName[1];
1484 0 : poDS->osFilename += ':';
1485 0 : poDS->osFilename += papszName[2];
1486 0 : poDS->osSubdatasetName = papszName[3];
1487 0 : poDS->bTreatAsSubdataset = TRUE;
1488 0 : CSLDestroy( papszName );
1489 : }
1490 3 : else if( CSLCount(papszName) == 3 )
1491 : {
1492 3 : poDS->osFilename = papszName[1];
1493 6 : poDS->osSubdatasetName = papszName[2];
1494 3 : poDS->bTreatAsSubdataset = TRUE;
1495 3 : CSLDestroy( papszName );
1496 : }
1497 : else
1498 : {
1499 0 : CSLDestroy( papszName );
1500 0 : delete poDS;
1501 : CPLError( CE_Failure, CPLE_AppDefined,
1502 0 : "Failed to parse NETCDF: prefix string into expected three fields." );
1503 0 : return NULL;
1504 : }
1505 : }
1506 : else
1507 : {
1508 18 : poDS->osFilename = poOpenInfo->pszFilename;
1509 18 : poDS->bTreatAsSubdataset = FALSE;
1510 : }
1511 :
1512 : /* -------------------------------------------------------------------- */
1513 : /* Try opening the dataset. */
1514 : /* -------------------------------------------------------------------- */
1515 21 : if( nc_open( poDS->osFilename, NC_NOWRITE, &cdfid ) != NC_NOERR ) {
1516 0 : delete poDS;
1517 0 : return NULL;
1518 : }
1519 : /* -------------------------------------------------------------------- */
1520 : /* Is this a real netCDF file? */
1521 : /* -------------------------------------------------------------------- */
1522 21 : status = nc_inq(cdfid, &ndims, &nvars, &ngatts, &unlimdimid);
1523 21 : if( status != NC_NOERR ) {
1524 0 : delete poDS;
1525 0 : return NULL;
1526 : }
1527 :
1528 : /* -------------------------------------------------------------------- */
1529 : /* Confirm the requested access is supported. */
1530 : /* -------------------------------------------------------------------- */
1531 21 : if( poOpenInfo->eAccess == GA_Update )
1532 : {
1533 : CPLError( CE_Failure, CPLE_NotSupported,
1534 : "The NETCDF driver does not support update access to existing"
1535 0 : " datasets.\n" );
1536 0 : nc_close( cdfid );
1537 0 : delete poDS;
1538 0 : return NULL;
1539 : }
1540 :
1541 : /* -------------------------------------------------------------------- */
1542 : /* Does the request variable exist? */
1543 : /* -------------------------------------------------------------------- */
1544 21 : if( poDS->bTreatAsSubdataset )
1545 : {
1546 3 : status = nc_inq_varid( cdfid, poDS->osSubdatasetName, &var);
1547 3 : if( status != NC_NOERR ) {
1548 : CPLError( CE_Warning, CPLE_AppDefined,
1549 : "%s is a netCDF file, but %s is not a variable.",
1550 : poOpenInfo->pszFilename,
1551 0 : poDS->osSubdatasetName.c_str() );
1552 :
1553 0 : nc_close( cdfid );
1554 0 : delete poDS;
1555 0 : return NULL;
1556 : }
1557 : }
1558 :
1559 21 : if( nc_inq_ndims( cdfid, &dim_count ) != NC_NOERR || dim_count < 2 )
1560 : {
1561 : CPLError( CE_Warning, CPLE_AppDefined,
1562 : "%s is a netCDF file, but not in GMT configuration.",
1563 0 : poOpenInfo->pszFilename );
1564 :
1565 0 : nc_close( cdfid );
1566 0 : delete poDS;
1567 0 : return NULL;
1568 : }
1569 :
1570 21 : CPLDebug( "GDAL_netCDF", "dim_count = %d\n", dim_count );
1571 :
1572 21 : if( (status = nc_get_att_text( cdfid, NC_GLOBAL, "Conventions",
1573 : attname )) != NC_NOERR ) {
1574 : CPLError( CE_Warning, CPLE_AppDefined,
1575 3 : "No UNIDATA NC_GLOBAL:Conventions attribute");
1576 : /* note that 'Conventions' is always capital 'C' in CF spec*/
1577 : }
1578 :
1579 :
1580 : /* -------------------------------------------------------------------- */
1581 : /* Create band information objects. */
1582 : /* -------------------------------------------------------------------- */
1583 21 : if ( nc_inq_nvars ( cdfid, &var_count) != NC_NOERR )
1584 : {
1585 0 : delete poDS;
1586 0 : return NULL;
1587 : }
1588 :
1589 21 : CPLDebug( "GDAL_netCDF", "var_count = %d\n", var_count );
1590 :
1591 : /* -------------------------------------------------------------------- */
1592 : /* Create a corresponding GDALDataset. */
1593 : /* Create Netcdf Subdataset if filename as NETCDF tag */
1594 : /* -------------------------------------------------------------------- */
1595 21 : poDS->cdfid = cdfid;
1596 :
1597 21 : poDS->ReadAttributes( cdfid, NC_GLOBAL );
1598 :
1599 : /* -------------------------------------------------------------------- */
1600 : /* Verify if only one variable has 2 dimensions */
1601 : /* -------------------------------------------------------------------- */
1602 73 : for ( j = 0; j < nvars; j++ ) {
1603 :
1604 52 : nc_inq_varndims ( cdfid, j, &ndims );
1605 52 : if( ndims >= 2 ) {
1606 29 : nVarID=j;
1607 29 : nCount++;
1608 : }
1609 : }
1610 :
1611 : /* -------------------------------------------------------------------- */
1612 : /* We have more than one variable with 2 dimensions in the */
1613 : /* file, then treat this as a subdataset container dataset. */
1614 : /* -------------------------------------------------------------------- */
1615 21 : if( (nCount > 1) && !poDS->bTreatAsSubdataset )
1616 : {
1617 4 : poDS->CreateSubDatasetList();
1618 4 : poDS->SetMetadata( poDS->papszMetadata );
1619 4 : poDS->TryLoadXML();
1620 4 : return( poDS );
1621 : }
1622 :
1623 : /* -------------------------------------------------------------------- */
1624 : /* If we are not treating things as a subdataset, then capture */
1625 : /* the name of the single available variable as the subdataset. */
1626 : /* -------------------------------------------------------------------- */
1627 17 : if( !poDS->bTreatAsSubdataset ) // nCount must be 1!
1628 : {
1629 : char szVarName[NC_MAX_NAME];
1630 :
1631 14 : nc_inq_varname( cdfid, nVarID, szVarName);
1632 14 : poDS->osSubdatasetName = szVarName;
1633 : }
1634 :
1635 : /* -------------------------------------------------------------------- */
1636 : /* Open the NETCDF subdataset NETCDF:"filename":subdataset */
1637 : /* -------------------------------------------------------------------- */
1638 :
1639 17 : var=-1;
1640 17 : nc_inq_varid( cdfid, poDS->osSubdatasetName, &var);
1641 17 : nd = 0;
1642 17 : nc_inq_varndims ( cdfid, var, &nd );
1643 :
1644 17 : poDS->paDimIds = (int *)CPLCalloc(nd, sizeof( int ) );
1645 17 : poDS->panBandDimPos = ( int * ) CPLCalloc( nd, sizeof( int ) );
1646 :
1647 17 : nc_inq_vardimid( cdfid, var, poDS->paDimIds );
1648 :
1649 :
1650 : /* -------------------------------------------------------------------- */
1651 : /* Check fi somebody tried to pass a variable with less than 2D */
1652 : /* -------------------------------------------------------------------- */
1653 :
1654 17 : if ( nd < 2 ) {
1655 4 : CPLFree( poDS->paDimIds );
1656 4 : CPLFree( poDS->panBandDimPos );
1657 4 : delete poDS;
1658 4 : return NULL;
1659 : }
1660 :
1661 : /* -------------------------------------------------------------------- */
1662 : /* CF-1 Convention */
1663 : /* dimensions to appear in the relative order T, then Z, then Y, */
1664 : /* then X to the file. All other dimensions should, whenever */
1665 : /* possible, be placed to the left of the spatiotemporal */
1666 : /* dimensions. */
1667 : /* -------------------------------------------------------------------- */
1668 :
1669 : /* -------------------------------------------------------------------- */
1670 : /* Get X dimensions information */
1671 : /* -------------------------------------------------------------------- */
1672 13 : poDS->nDimXid = poDS->paDimIds[nd-1];
1673 13 : nc_inq_dimlen ( cdfid, poDS->nDimXid, &poDS->xdim );
1674 13 : poDS->nRasterXSize = poDS->xdim;
1675 :
1676 : /* -------------------------------------------------------------------- */
1677 : /* Get Y dimension information */
1678 : /* -------------------------------------------------------------------- */
1679 13 : poDS->nDimYid = poDS->paDimIds[nd-2];
1680 13 : nc_inq_dimlen ( cdfid, poDS->nDimYid, &poDS->ydim );
1681 13 : poDS->nRasterYSize = poDS->ydim;
1682 :
1683 :
1684 45 : for( j=0,k=0; j < nd; j++ ){
1685 32 : if( poDS->paDimIds[j] == poDS->nDimXid ){
1686 13 : poDS->panBandDimPos[0] = j; // Save Position of XDim
1687 13 : k++;
1688 : }
1689 32 : if( poDS->paDimIds[j] == poDS->nDimYid ){
1690 13 : poDS->panBandDimPos[1] = j; // Save Position of YDim
1691 13 : k++;
1692 : }
1693 : }
1694 : /* -------------------------------------------------------------------- */
1695 : /* X and Y Dimension Ids were not found! */
1696 : /* -------------------------------------------------------------------- */
1697 13 : if( k != 2 ) {
1698 0 : CPLFree( poDS->paDimIds );
1699 0 : CPLFree( poDS->panBandDimPos );
1700 0 : return NULL;
1701 : }
1702 :
1703 : /* -------------------------------------------------------------------- */
1704 : /* Read Metadata for this variable */
1705 : /* -------------------------------------------------------------------- */
1706 13 : poDS->ReadAttributes( cdfid, var );
1707 :
1708 : /* -------------------------------------------------------------------- */
1709 : /* Read Metadata for each dimension */
1710 : /* -------------------------------------------------------------------- */
1711 :
1712 45 : for( j=0; j < dim_count; j++ ){
1713 32 : nc_inq_dimname( cdfid, j, poDS->papszDimName[j] );
1714 32 : status = nc_inq_varid( cdfid, poDS->papszDimName[j], &nDimID );
1715 32 : if( status == NC_NOERR ) {
1716 6 : poDS->ReadAttributes( cdfid, nDimID );
1717 : }
1718 : }
1719 :
1720 13 : poDS->SetProjection( var );
1721 13 : poDS->SetMetadata( poDS->papszMetadata );
1722 :
1723 : /* -------------------------------------------------------------------- */
1724 : /* Create bands */
1725 : /* -------------------------------------------------------------------- */
1726 13 : poDS->panBandZLev = (int *)CPLCalloc( nd-2, sizeof( int ) );
1727 :
1728 13 : nTotLevCount = 1;
1729 13 : if ( dim_count > 2 ) {
1730 2 : nDim=2;
1731 12 : for( j=0; j < nd; j++ ){
1732 18 : if( ( poDS->paDimIds[j] != poDS->nDimXid ) &&
1733 8 : ( poDS->paDimIds[j] != poDS->nDimYid ) ){
1734 6 : nc_inq_dimlen ( cdfid, poDS->paDimIds[j], &lev_count );
1735 6 : nTotLevCount *= lev_count;
1736 6 : poDS->panBandZLev[ nDim-2 ] = lev_count;
1737 6 : poDS->panBandDimPos[ nDim++ ] = j; //Save Position of ZDim
1738 : }
1739 : }
1740 : }
1741 13 : i=0;
1742 :
1743 80 : for ( unsigned int lev = 0; lev < nTotLevCount ; lev++ ) {
1744 : char ** papszToken;
1745 27 : papszToken=NULL;
1746 :
1747 : netCDFRasterBand *poBand =
1748 : new netCDFRasterBand(poDS, var, nDim, lev,
1749 27 : poDS->panBandZLev, poDS->panBandDimPos, i+1 );
1750 :
1751 27 : poDS->SetBand( i+1, poBand );
1752 27 : i++;
1753 : }
1754 :
1755 13 : CPLFree( poDS->paDimIds );
1756 13 : CPLFree( poDS->panBandDimPos );
1757 13 : CPLFree( poDS->panBandZLev );
1758 :
1759 13 : poDS->nBands = i;
1760 :
1761 : // Handle angular geographic coordinates here
1762 :
1763 : /* -------------------------------------------------------------------- */
1764 : /* Initialize any PAM information. */
1765 : /* -------------------------------------------------------------------- */
1766 13 : if( poDS->bTreatAsSubdataset )
1767 : {
1768 3 : poDS->SetPhysicalFilename( poDS->osFilename );
1769 3 : poDS->SetSubdatasetName( poDS->osSubdatasetName );
1770 : }
1771 :
1772 13 : poDS->TryLoadXML();
1773 :
1774 13 : if( poDS->bTreatAsSubdataset )
1775 3 : poDS->oOvManager.Initialize( poDS, ":::VIRTUAL:::" );
1776 : else
1777 10 : poDS->oOvManager.Initialize( poDS, poDS->osFilename );
1778 :
1779 13 : return( poDS );
1780 : }
1781 :
1782 :
1783 : /************************************************************************/
1784 : /* CopyMetadata() */
1785 : /* */
1786 : /* Create a copy of metadata for NC_GLOBAL or a variable */
1787 : /************************************************************************/
1788 :
1789 42 : void CopyMetadata( void *poDS, int fpImage, int CDFVarID ) {
1790 :
1791 : /* -------------------------------------------------------------------- */
1792 : /* Add CF-1.0 Conventions Global attribute */
1793 : /* -------------------------------------------------------------------- */
1794 : char **papszMetadata;
1795 : char **papszFieldData;
1796 : const char *pszField;
1797 : char szMetaName[ MAX_STR_LEN ];
1798 : char szMetaValue[ MAX_STR_LEN ];
1799 : int nDataLength;
1800 : int nItems;
1801 :
1802 : /* -------------------------------------------------------------------- */
1803 : /* Global metadata are set with NC_GLOBAL as the varid */
1804 : /* -------------------------------------------------------------------- */
1805 :
1806 42 : if( CDFVarID == NC_GLOBAL ) {
1807 :
1808 16 : papszMetadata = GDALGetMetadata( (GDALDataset *) poDS,"");
1809 :
1810 : nc_put_att_text( fpImage,
1811 : NC_GLOBAL,
1812 : "Conventions",
1813 : 6,
1814 16 : "CF-1.0" );
1815 : } else {
1816 :
1817 26 : papszMetadata = GDALGetMetadata( (GDALRasterBandH) poDS, NULL );
1818 :
1819 : }
1820 :
1821 42 : nItems = CSLCount( papszMetadata );
1822 :
1823 57 : for(int k=0; k < nItems; k++ ) {
1824 15 : pszField = CSLGetField( papszMetadata, k );
1825 : papszFieldData = CSLTokenizeString2 (pszField, "=",
1826 15 : CSLT_HONOURSTRINGS );
1827 15 : if( papszFieldData[1] != NULL ) {
1828 15 : strcpy( szMetaName, papszFieldData[ 0 ] );
1829 15 : strcpy( szMetaValue, papszFieldData[ 1 ] );
1830 :
1831 : /* -------------------------------------------------------------------- */
1832 : /* netCDF attributes do not like the '#' character. */
1833 : /* -------------------------------------------------------------------- */
1834 :
1835 195 : for( unsigned int h=0; h < strlen( szMetaName ) -1 ; h++ ) {
1836 180 : if( szMetaName[h] == '#' ) szMetaName[h] = '-';
1837 : }
1838 :
1839 15 : nDataLength = strlen( szMetaValue );
1840 : nc_put_att_text( fpImage,
1841 : CDFVarID,
1842 : szMetaName,
1843 : nDataLength,
1844 15 : szMetaValue );
1845 :
1846 :
1847 : }
1848 15 : CSLDestroy( papszFieldData );
1849 :
1850 : }
1851 42 : }
1852 : /************************************************************************/
1853 : /* CreateCopy() */
1854 : /************************************************************************/
1855 :
1856 :
1857 : static GDALDataset*
1858 17 : NCDFCreateCopy( const char * pszFilename, GDALDataset *poSrcDS,
1859 : int bStrict, char ** papszOptions,
1860 : GDALProgressFunc pfnProgress, void * pProgressData )
1861 :
1862 : {
1863 17 : int nBands = poSrcDS->GetRasterCount();
1864 17 : int nXSize = poSrcDS->GetRasterXSize();
1865 17 : int nYSize = poSrcDS->GetRasterYSize();
1866 17 : int bProgressive = FALSE;
1867 :
1868 : int anBandDims[ NC_MAX_DIMS ];
1869 : int anBandMap[ NC_MAX_DIMS ];
1870 :
1871 17 : int bWriteGeoTransform = FALSE;
1872 : char pszNetcdfProjection[ NC_MAX_NAME ];
1873 :
1874 17 : if (nBands == 0)
1875 : {
1876 : CPLError( CE_Failure, CPLE_NotSupported,
1877 1 : "NetCDF driver does not support source dataset with zero band.\n");
1878 1 : return NULL;
1879 : }
1880 :
1881 16 : if( !pfnProgress( 0.0, NULL, pProgressData ) )
1882 0 : return NULL;
1883 :
1884 :
1885 16 : bProgressive = CSLFetchBoolean( papszOptions, "PROGRESSIVE", FALSE );
1886 :
1887 : /* -------------------------------------------------------------------- */
1888 : /* Create the dataset. */
1889 : /* -------------------------------------------------------------------- */
1890 : int fpImage;
1891 : int status;
1892 16 : int nXDimID = 0;
1893 16 : int nYDimID = 0;
1894 :
1895 16 : status = nc_create( pszFilename, NC_CLOBBER, &fpImage );
1896 :
1897 16 : if( status != NC_NOERR )
1898 : {
1899 : CPLError( CE_Failure, CPLE_OpenFailed,
1900 : "Unable to create netCDF file %s.\n",
1901 0 : pszFilename );
1902 0 : return NULL;
1903 : }
1904 :
1905 : GDALDataType eDT;
1906 :
1907 16 : status = nc_def_dim( fpImage, "x", nXSize, &nXDimID );
1908 16 : CPLDebug( "GDAL_netCDF", "status nc_def_dim X = %d\n", status );
1909 :
1910 16 : status = nc_def_dim( fpImage, "y", nYSize, &nYDimID );
1911 16 : CPLDebug( "GDAL_netCDF", "status nc_def_dim Y = %d\n", status );
1912 :
1913 16 : CPLDebug( "GDAL_netCDF", "nYDimID = %d\n", nXDimID );
1914 16 : CPLDebug( "GDAL_netCDF", "nXDimID = %d\n", nYDimID );
1915 16 : CPLDebug( "GDAL_netCDF", "nXSize = %d\n", nXSize );
1916 16 : CPLDebug( "GDAL_netCDF", "nYSize = %d\n", nYSize );
1917 :
1918 :
1919 16 : CopyMetadata((void *) poSrcDS, fpImage, NC_GLOBAL );
1920 :
1921 : /* -------------------------------------------------------------------- */
1922 : /* Set Projection for netCDF data CF-1 Convention */
1923 : /* -------------------------------------------------------------------- */
1924 :
1925 16 : OGRSpatialReference oSRS;
1926 16 : char *pszWKT = (char *) poSrcDS->GetProjectionRef();
1927 :
1928 :
1929 16 : if( pszWKT != NULL )
1930 16 : oSRS.importFromWkt( &pszWKT );
1931 :
1932 16 : if( oSRS.IsGeographic() ) {
1933 :
1934 : int status;
1935 : int i;
1936 : int NCDFVarID;
1937 :
1938 15 : double dfNN=0.0;
1939 15 : double dfSN=0.0;
1940 15 : double dfEE=0.0;
1941 15 : double dfWE=0.0;
1942 : double adfGeoTransform[6];
1943 : char szGeoTransform[ MAX_STR_LEN ];
1944 : char szTemp[ MAX_STR_LEN ];
1945 :
1946 :
1947 : /* -------------------------------------------------------------------- */
1948 : /* Copy GeoTransform array from source */
1949 : /* -------------------------------------------------------------------- */
1950 15 : poSrcDS->GetGeoTransform( adfGeoTransform );
1951 15 : bWriteGeoTransform = TRUE;
1952 :
1953 15 : *szGeoTransform = '\0';
1954 105 : for( i=0; i<6; i++ ) {
1955 : sprintf( szTemp, "%g ",
1956 90 : adfGeoTransform[i] );
1957 90 : strcat( szGeoTransform, szTemp );
1958 : }
1959 15 : CPLDebug( "GDAL_netCDF", "szGeoTranform = %s", szGeoTransform );
1960 15 : strcpy( pszNetcdfProjection, "GDAL_Geographics" );
1961 :
1962 : status = nc_def_var( fpImage,
1963 : pszNetcdfProjection,
1964 : NC_CHAR,
1965 15 : 0, NULL, &NCDFVarID );
1966 :
1967 15 : dfNN = adfGeoTransform[3];
1968 15 : dfSN = ( adfGeoTransform[5] * nYSize ) + dfNN;
1969 15 : dfWE = adfGeoTransform[0];
1970 15 : dfEE = ( adfGeoTransform[1] * nXSize ) + dfWE;
1971 :
1972 : status = nc_put_att_double( fpImage,
1973 : NCDFVarID,
1974 : "Northernmost_Northing",
1975 : NC_DOUBLE,
1976 : 1,
1977 15 : &dfNN );
1978 : status = nc_put_att_double( fpImage,
1979 : NCDFVarID,
1980 : "Southernmost_Northing",
1981 : NC_DOUBLE,
1982 : 1,
1983 15 : &dfSN );
1984 : status = nc_put_att_double( fpImage,
1985 : NCDFVarID,
1986 : "Easternmost_Easting",
1987 : NC_DOUBLE,
1988 : 1,
1989 15 : &dfEE );
1990 : status = nc_put_att_double( fpImage,
1991 : NCDFVarID,
1992 : "Westernmost_Easting",
1993 : NC_DOUBLE,
1994 : 1,
1995 15 : &dfWE );
1996 15 : pszWKT = (char *) poSrcDS->GetProjectionRef() ;
1997 :
1998 : nc_put_att_text( fpImage,
1999 : NCDFVarID,
2000 : "spatial_ref",
2001 : strlen( pszWKT ),
2002 15 : pszWKT );
2003 :
2004 : nc_put_att_text( fpImage,
2005 : NCDFVarID,
2006 : "GeoTransform",
2007 : strlen( szGeoTransform ),
2008 15 : szGeoTransform );
2009 :
2010 : nc_put_att_text( fpImage,
2011 : NCDFVarID,
2012 : GRD_MAPPING_NAME,
2013 : 30,
2014 15 : "Geographics Coordinate System" );
2015 : nc_put_att_text( fpImage,
2016 : NCDFVarID,
2017 : LNG_NAME,
2018 : 14,
2019 15 : "Grid_latitude" );
2020 : }
2021 :
2022 16 : if( oSRS.IsProjected() )
2023 : {
2024 : const char *pszParamStr, *pszParamVal;
2025 1 : const OGR_SRSNode *poPROJCS = oSRS.GetAttrNode( "PROJCS" );
2026 : int status;
2027 : int i;
2028 : int NCDFVarID;
2029 : const char *pszProjection;
2030 1 : double dfNN=0.0;
2031 1 : double dfSN=0.0;
2032 1 : double dfEE=0.0;
2033 1 : double dfWE=0.0;
2034 : double adfGeoTransform[6];
2035 : char szGeoTransform[ MAX_STR_LEN ];
2036 : char szTemp[ MAX_STR_LEN ];
2037 :
2038 1 : poSrcDS->GetGeoTransform( adfGeoTransform );
2039 :
2040 1 : *szGeoTransform = '\0';
2041 7 : for( i=0; i<6; i++ ) {
2042 : sprintf( szTemp, "%g ",
2043 6 : adfGeoTransform[i] );
2044 6 : strcat( szGeoTransform, szTemp );
2045 : }
2046 :
2047 1 : CPLDebug( "GDAL_netCDF", "szGeoTranform = %s", szGeoTransform );
2048 :
2049 1 : pszProjection = oSRS.GetAttrValue( "PROJECTION" );
2050 1 : bWriteGeoTransform = TRUE;
2051 :
2052 :
2053 32 : for(i=0; poNetcdfSRS[i].netCDFSRS != NULL; i++ ) {
2054 32 : if( EQUAL( poNetcdfSRS[i].SRS, pszProjection ) ) {
2055 : CPLDebug( "GDAL_netCDF", "PROJECTION = %s",
2056 1 : poNetcdfSRS[i].netCDFSRS);
2057 1 : strcpy( pszNetcdfProjection, poNetcdfSRS[i].netCDFSRS );
2058 :
2059 1 : break;
2060 : }
2061 : }
2062 :
2063 : status = nc_def_var( fpImage,
2064 : poNetcdfSRS[i].netCDFSRS,
2065 : NC_CHAR,
2066 1 : 0, NULL, &NCDFVarID );
2067 :
2068 1 : dfNN = adfGeoTransform[3];
2069 1 : dfSN = ( adfGeoTransform[5] * nYSize ) + dfNN;
2070 1 : dfWE = adfGeoTransform[0];
2071 1 : dfEE = ( adfGeoTransform[1] * nXSize ) + dfWE;
2072 :
2073 : status = nc_put_att_double( fpImage,
2074 : NCDFVarID,
2075 : "Northernmost_Northing",
2076 : NC_DOUBLE,
2077 : 1,
2078 1 : &dfNN );
2079 : status = nc_put_att_double( fpImage,
2080 : NCDFVarID,
2081 : "Southernmost_Northing",
2082 : NC_DOUBLE,
2083 : 1,
2084 1 : &dfSN );
2085 : status = nc_put_att_double( fpImage,
2086 : NCDFVarID,
2087 : "Easternmost_Easting",
2088 : NC_DOUBLE,
2089 : 1,
2090 1 : &dfEE );
2091 : status = nc_put_att_double( fpImage,
2092 : NCDFVarID,
2093 : "Westernmost_Easting",
2094 : NC_DOUBLE,
2095 : 1,
2096 1 : &dfWE );
2097 1 : pszWKT = (char *) poSrcDS->GetProjectionRef() ;
2098 :
2099 : nc_put_att_text( fpImage,
2100 : NCDFVarID,
2101 : "spatial_ref",
2102 : strlen( pszWKT ),
2103 1 : pszWKT );
2104 :
2105 : nc_put_att_text( fpImage,
2106 : NCDFVarID,
2107 : "GeoTransform",
2108 : strlen( szGeoTransform ),
2109 1 : szGeoTransform );
2110 :
2111 : nc_put_att_text( fpImage,
2112 : NCDFVarID,
2113 : GRD_MAPPING_NAME,
2114 : strlen( pszNetcdfProjection ),
2115 1 : pszNetcdfProjection );
2116 :
2117 :
2118 11 : for( int iChild = 0; iChild < poPROJCS->GetChildCount(); iChild++ )
2119 : {
2120 : const OGR_SRSNode *poNode;
2121 : float fValue;
2122 :
2123 10 : poNode = poPROJCS->GetChild( iChild );
2124 10 : if( !EQUAL(poNode->GetValue(),"PARAMETER")
2125 : || poNode->GetChildCount() != 2 )
2126 5 : continue;
2127 :
2128 : /* -------------------------------------------------------------------- */
2129 : /* Look for projection attributes */
2130 : /* -------------------------------------------------------------------- */
2131 5 : pszParamStr = poNode->GetChild(0)->GetValue();
2132 5 : pszParamVal = poNode->GetChild(1)->GetValue();
2133 :
2134 :
2135 192 : for(i=0; poNetcdfSRS[i].netCDFSRS != NULL; i++ ) {
2136 192 : if( EQUAL( poNetcdfSRS[i].SRS, pszParamStr ) ) {
2137 : CPLDebug( "GDAL_netCDF", "%s = %s",
2138 : poNetcdfSRS[i].netCDFSRS,
2139 5 : pszParamVal );
2140 5 : break;
2141 : }
2142 : }
2143 : /* -------------------------------------------------------------------- */
2144 : /* Write Projection attribute */
2145 : /* -------------------------------------------------------------------- */
2146 5 : sscanf( pszParamVal, "%f", &fValue );
2147 5 : if( poNetcdfSRS[i].netCDFSRS != NULL ) {
2148 : nc_put_att_float( fpImage,
2149 : NCDFVarID,
2150 : poNetcdfSRS[i].netCDFSRS,
2151 : NC_FLOAT,
2152 : 1,
2153 5 : &fValue );
2154 :
2155 : }
2156 : }
2157 : }
2158 :
2159 : /* -------------------------------------------------------------------- */
2160 : /* Initialize Band Map */
2161 : /* -------------------------------------------------------------------- */
2162 :
2163 42 : for(int j=1; j <= nBands; j++ ) {
2164 26 : anBandMap[j-1]=j;
2165 : }
2166 :
2167 : /* -------------------------------------------------------------------- */
2168 : /* Create netCDF variable */
2169 : /* -------------------------------------------------------------------- */
2170 :
2171 42 : for( int i=1; i <= nBands; i++ ) {
2172 :
2173 : char szBandName[ NC_MAX_NAME ];
2174 26 : GByte *pabScanline = NULL;
2175 26 : GInt16 *pasScanline = NULL;
2176 26 : GInt32 *panScanline = NULL;
2177 26 : float *pafScanline = NULL;
2178 26 : double *padScanline = NULL;
2179 : int NCDFVarID;
2180 : size_t start[ GDALNBDIM ];
2181 : size_t count[ GDALNBDIM ];
2182 : size_t nBandNameLen;
2183 : double dfNoDataValue;
2184 : unsigned char cNoDataValue;
2185 : float fNoDataValue;
2186 : int nlNoDataValue;
2187 : short nsNoDataValue;
2188 : GDALRasterBandH hBand;
2189 :
2190 26 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand( i );
2191 26 : hBand = GDALGetRasterBand( poSrcDS, i );
2192 :
2193 26 : sprintf( szBandName, "Band%d", i );
2194 :
2195 26 : eDT = poSrcDS->GetRasterBand(i)->GetRasterDataType();
2196 26 : anBandDims[0] = nYDimID;
2197 26 : anBandDims[1] = nXDimID;
2198 26 : CPLErr eErr = CE_None;
2199 :
2200 26 : dfNoDataValue = poSrcBand->GetNoDataValue(0);
2201 :
2202 26 : if( eDT == GDT_Byte ) {
2203 16 : CPLDebug( "GDAL_netCDF", "%s = GDT_Byte ", szBandName );
2204 :
2205 : status = nc_def_var( fpImage, szBandName, NC_BYTE,
2206 16 : GDALNBDIM, anBandDims, &NCDFVarID );
2207 :
2208 :
2209 : /* -------------------------------------------------------------------- */
2210 : /* Write data line per line */
2211 : /* -------------------------------------------------------------------- */
2212 :
2213 16 : pabScanline = (GByte *) CPLMalloc( nBands * nXSize );
2214 186 : for( int iLine = 0; iLine < nYSize && eErr == CE_None; iLine++ ) {
2215 :
2216 : /* -------------------------------------------------------------------- */
2217 : /* Read data from band i */
2218 : /* -------------------------------------------------------------------- */
2219 : eErr = poSrcBand->RasterIO( GF_Read, 0, iLine, nXSize, 1,
2220 : pabScanline, nXSize, 1, GDT_Byte,
2221 170 : 0,0);
2222 :
2223 : /* -------------------------------------------------------------------- */
2224 : /* Write Fill Value */
2225 : /* -------------------------------------------------------------------- */
2226 :
2227 170 : cNoDataValue=(unsigned char) dfNoDataValue;
2228 : nc_put_att_uchar( fpImage,
2229 : NCDFVarID,
2230 : _FillValue,
2231 : NC_CHAR,
2232 : 1,
2233 170 : &cNoDataValue );
2234 :
2235 : /* -------------------------------------------------------------------- */
2236 : /* Write Data from Band i */
2237 : /* -------------------------------------------------------------------- */
2238 170 : start[0]=iLine;
2239 170 : start[1]=0;
2240 170 : count[0]=1;
2241 170 : count[1]=nXSize;
2242 :
2243 : /* -------------------------------------------------------------------- */
2244 : /* Put NetCDF file in data mode. */
2245 : /* -------------------------------------------------------------------- */
2246 170 : status = nc_enddef( fpImage );
2247 : status = nc_put_vara_uchar (fpImage, NCDFVarID, start,
2248 170 : count, pabScanline);
2249 :
2250 :
2251 : /* -------------------------------------------------------------------- */
2252 : /* Put NetCDF file back in define mode. */
2253 : /* -------------------------------------------------------------------- */
2254 170 : status = nc_redef( fpImage );
2255 :
2256 : }
2257 16 : CPLFree( pabScanline );
2258 : /* -------------------------------------------------------------------- */
2259 : /* Int16 */
2260 : /* -------------------------------------------------------------------- */
2261 :
2262 12 : } else if( ( eDT == GDT_UInt16 ) || ( eDT == GDT_Int16 ) ) {
2263 2 : CPLDebug( "GDAL_netCDF", "%s = GDT_Int16 ",szBandName );
2264 : status = nc_def_var( fpImage, szBandName, NC_SHORT,
2265 2 : GDALNBDIM, anBandDims, &NCDFVarID );
2266 :
2267 : pasScanline = (GInt16 *) CPLMalloc( nBands * nXSize *
2268 2 : sizeof( GInt16 ) );
2269 : /* -------------------------------------------------------------------- */
2270 : /* Write Fill Value */
2271 : /* -------------------------------------------------------------------- */
2272 2 : nsNoDataValue= (GInt16) dfNoDataValue;
2273 : nc_put_att_short( fpImage,
2274 : NCDFVarID,
2275 : _FillValue,
2276 : NC_SHORT,
2277 : 1,
2278 2 : &nsNoDataValue );
2279 22 : for( int iLine = 0; iLine < nYSize && eErr == CE_None; iLine++ ) {
2280 :
2281 : eErr = poSrcBand->RasterIO( GF_Read, 0, iLine, nXSize, 1,
2282 : pasScanline, nXSize, 1, GDT_Int16,
2283 20 : 0,0);
2284 :
2285 20 : start[0]=iLine;
2286 20 : start[1]=0;
2287 20 : count[0]=1;
2288 20 : count[1]=nXSize;
2289 :
2290 :
2291 20 : status = nc_enddef( fpImage );
2292 : status = nc_put_vara_short( fpImage, NCDFVarID, start,
2293 20 : count, pasScanline);
2294 20 : status = nc_redef( fpImage );
2295 : }
2296 2 : CPLFree( pasScanline );
2297 : /* -------------------------------------------------------------------- */
2298 : /* Int32 */
2299 : /* -------------------------------------------------------------------- */
2300 :
2301 10 : } else if( (eDT == GDT_UInt32) || (eDT == GDT_Int32) ) {
2302 2 : CPLDebug( "GDAL_netCDF", "%s = GDT_Int32 ",szBandName );
2303 : status = nc_def_var( fpImage, szBandName, NC_INT,
2304 2 : GDALNBDIM, anBandDims, &NCDFVarID );
2305 :
2306 : panScanline = (GInt32 *) CPLMalloc( nBands * nXSize *
2307 2 : sizeof( GInt32 ) );
2308 : /* -------------------------------------------------------------------- */
2309 : /* Write Fill Value */
2310 : /* -------------------------------------------------------------------- */
2311 2 : nlNoDataValue= (GInt32) dfNoDataValue;
2312 :
2313 : nc_put_att_int( fpImage,
2314 : NCDFVarID,
2315 : _FillValue,
2316 : NC_INT,
2317 : 1,
2318 2 : &nlNoDataValue );
2319 :
2320 :
2321 22 : for( int iLine = 0; iLine < nYSize && eErr == CE_None; iLine++ ) {
2322 :
2323 : eErr = poSrcBand->RasterIO( GF_Read, 0, iLine, nXSize, 1,
2324 : panScanline, nXSize, 1, GDT_Int32,
2325 20 : 0,0);
2326 :
2327 20 : start[0]=iLine;
2328 20 : start[1]=0;
2329 20 : count[0]=1;
2330 20 : count[1]=nXSize;
2331 :
2332 :
2333 20 : status = nc_enddef( fpImage );
2334 : status = nc_put_vara_int( fpImage, NCDFVarID, start,
2335 20 : count, panScanline);
2336 20 : status = nc_redef( fpImage );
2337 : }
2338 2 : CPLFree( panScanline );
2339 : /* -------------------------------------------------------------------- */
2340 : /* float */
2341 : /* -------------------------------------------------------------------- */
2342 6 : } else if( (eDT == GDT_Float32) ) {
2343 1 : CPLDebug( "GDAL_netCDF", "%s = GDT_Float32 ",szBandName );
2344 : status = nc_def_var( fpImage, szBandName, NC_FLOAT,
2345 1 : GDALNBDIM, anBandDims, &NCDFVarID );
2346 :
2347 : pafScanline = (float *) CPLMalloc( nBands * nXSize *
2348 1 : sizeof( float ) );
2349 :
2350 : /* -------------------------------------------------------------------- */
2351 : /* Write Fill Value */
2352 : /* -------------------------------------------------------------------- */
2353 1 : fNoDataValue= (float) dfNoDataValue;
2354 : nc_put_att_float( fpImage,
2355 : NCDFVarID,
2356 : _FillValue,
2357 : NC_FLOAT,
2358 : 1,
2359 1 : &fNoDataValue );
2360 :
2361 11 : for( int iLine = 0; iLine < nYSize && eErr == CE_None; iLine++ ) {
2362 :
2363 : eErr = poSrcBand->RasterIO( GF_Read, 0, iLine, nXSize, 1,
2364 : pafScanline, nXSize, 1,
2365 : GDT_Float32,
2366 10 : 0,0);
2367 :
2368 10 : start[0]=iLine;
2369 10 : start[1]=0;
2370 10 : count[0]=1;
2371 10 : count[1]=nXSize;
2372 :
2373 :
2374 10 : status = nc_enddef( fpImage );
2375 : status = nc_put_vara_float( fpImage, NCDFVarID, start,
2376 10 : count, pafScanline);
2377 10 : status = nc_redef( fpImage );
2378 : }
2379 1 : CPLFree( pafScanline );
2380 : /* -------------------------------------------------------------------- */
2381 : /* double */
2382 : /* -------------------------------------------------------------------- */
2383 5 : } else if( (eDT == GDT_Float64) ) {
2384 1 : CPLDebug( "GDAL_netCDF", "%s = GDT_Float64 ",szBandName );
2385 : status = nc_def_var( fpImage, szBandName, NC_DOUBLE,
2386 1 : GDALNBDIM, anBandDims, &NCDFVarID );
2387 :
2388 : padScanline = (double *) CPLMalloc( nBands * nXSize *
2389 1 : sizeof( double ) );
2390 :
2391 : /* -------------------------------------------------------------------- */
2392 : /* Write Fill Value */
2393 : /* -------------------------------------------------------------------- */
2394 :
2395 : nc_put_att_double( fpImage,
2396 : NCDFVarID,
2397 : _FillValue,
2398 : NC_DOUBLE,
2399 : 1,
2400 1 : &dfNoDataValue );
2401 :
2402 11 : for( int iLine = 0; iLine < nYSize && eErr == CE_None; iLine++ ) {
2403 :
2404 : eErr = poSrcBand->RasterIO( GF_Read, 0, iLine, nXSize, 1,
2405 : padScanline, nXSize, 1,
2406 : GDT_Float64,
2407 10 : 0,0);
2408 :
2409 10 : start[0]=iLine;
2410 10 : start[1]=0;
2411 10 : count[0]=1;
2412 10 : count[1]=nXSize;
2413 :
2414 :
2415 10 : status = nc_enddef( fpImage );
2416 : status = nc_put_vara_double( fpImage, NCDFVarID, start,
2417 10 : count, padScanline);
2418 10 : status = nc_redef( fpImage );
2419 : }
2420 1 : CPLFree( padScanline );
2421 : }
2422 :
2423 : /* -------------------------------------------------------------------- */
2424 : /* Write Projection for band */
2425 : /* -------------------------------------------------------------------- */
2426 26 : if( bWriteGeoTransform == TRUE ) {
2427 : /* nc_put_att_text( fpImage, NCDFVarID,
2428 : COORDINATES,
2429 : 7,
2430 : LONLAT );
2431 : */
2432 : nc_put_att_text( fpImage, NCDFVarID,
2433 : GRD_MAPPING,
2434 : strlen( pszNetcdfProjection ),
2435 26 : pszNetcdfProjection );
2436 : }
2437 :
2438 26 : sprintf( szBandName, "GDAL Band Number %d", i);
2439 26 : nBandNameLen = strlen( szBandName );
2440 : nc_put_att_text( fpImage,
2441 : NCDFVarID,
2442 : "long_name",
2443 : nBandNameLen,
2444 26 : szBandName );
2445 :
2446 26 : CopyMetadata( (void *) hBand, fpImage, NCDFVarID );
2447 :
2448 :
2449 :
2450 : }
2451 :
2452 :
2453 :
2454 : // poDstDS->SetGeoTransform( adfGeoTransform );
2455 :
2456 :
2457 : /* -------------------------------------------------------------------- */
2458 : /* Cleanup and close. */
2459 : /* -------------------------------------------------------------------- */
2460 : // CPLFree( pabScanline );
2461 :
2462 16 : nc_close( fpImage );
2463 : /* -------------------------------------------------------------------- */
2464 : /* Re-open dataset, and copy any auxilary pam information. */
2465 : /* -------------------------------------------------------------------- */
2466 16 : netCDFDataset *poDS = (netCDFDataset *) GDALOpen( pszFilename, GA_ReadOnly );
2467 :
2468 16 : if( poDS )
2469 12 : poDS->CloneInfo( poSrcDS, GCIF_PAM_DEFAULT );
2470 :
2471 16 : return poDS;
2472 :
2473 :
2474 : }
2475 :
2476 : /************************************************************************/
2477 : /* GDALRegister_netCDF() */
2478 : /************************************************************************/
2479 :
2480 338 : void GDALRegister_netCDF()
2481 :
2482 : {
2483 : GDALDriver *poDriver;
2484 :
2485 338 : if (! GDAL_CHECK_VERSION("netCDF driver"))
2486 0 : return;
2487 :
2488 338 : if( GDALGetDriverByName( "netCDF" ) == NULL )
2489 : {
2490 336 : poDriver = new GDALDriver( );
2491 :
2492 336 : poDriver->SetDescription( "netCDF" );
2493 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
2494 336 : "Network Common Data Format" );
2495 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
2496 336 : "frmt_netcdf.html" );
2497 336 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "nc" );
2498 :
2499 336 : poDriver->pfnOpen = netCDFDataset::Open;
2500 336 : poDriver->pfnCreateCopy = NCDFCreateCopy;
2501 :
2502 :
2503 336 : GetGDALDriverManager( )->RegisterDriver( poDriver );
2504 : }
2505 : }
2506 :
2507 :
2508 :
2509 : /* -------------------------------------------------------------------- */
2510 : /* Set Lambert Conformal Conic Projection */
2511 : /* -------------------------------------------------------------------- */
2512 :
2513 :
2514 :
2515 : //Albers equal area
2516 : //
2517 : //grid_mapping_name = albers_conical_equal_area
2518 : //
2519 : //Map parameters:
2520 : //
2521 : // * standard_parallel - There may be 1 or 2 values.
2522 : // * longitude_of_central_meridian
2523 : // * latitude_of_projection_origin
2524 : // * false_easting
2525 : // * false_northing
2526 : //Lambert azimuthal equal area
2527 : //
2528 : //grid_mapping_name = lambert_azimuthal_equal_area
2529 : //
2530 : //Map parameters:
2531 : //
2532 : // * longitude_of_projection_origin
2533 : // * latitude_of_projection_origin
2534 : // * false_easting
2535 : // * false_northing
2536 : //Lambert conformal
2537 : //
2538 : //grid_mapping_name = lambert_conformal_conic
2539 : //
2540 : //Map parameters:
2541 : //
2542 : // * standard_parallel - There may be 1 or 2 values.
2543 : // * longitude_of_central_meridian
2544 : // * latitude_of_projection_origin
2545 : // * false_easting
2546 : // * false_northing
2547 : //Polar stereographic
2548 : //
2549 : //grid_mapping_name = polar_stereographic
2550 : //
2551 : //Map parameters:
2552 : //
2553 : // * straight_vertical_longitude_from_pole
2554 : // * latitude_of_projection_origin - Either +90. or -90.
2555 : // * Either standard_parallel or scale_factor_at_projection_origin
2556 : // * false_easting
2557 : // * false_northing
2558 : //Rotated pole
2559 : //
2560 : //grid_mapping_name = rotated_latitude_longitude
2561 : //
2562 : //Map parameters:
2563 : //
2564 : // * grid_north_pole_latitude
2565 : // * grid_north_pole_longitude
2566 : // * north_pole_grid_longitude - This parameter is optional (default is 0.).
2567 : //Stereographic
2568 : //
2569 : //grid_mapping_name = stereographic
2570 : //
2571 : //Map parameters:
2572 : //
2573 : // * longitude_of_projection_origin
2574 : // * latitude_of_projection_origin
2575 : // * scale_factor_at_projection_origin
2576 : // * false_easting
2577 : // * false_northing
2578 : //Transverse Mercator
2579 : //
2580 : //grid_mapping_name = transverse_mercator
2581 : //
2582 : //Map parameters:
2583 : //
2584 : // * scale_factor_at_central_meridian
2585 : // * longitude_of_central_meridian
2586 : // * latitude_of_projection_origin
2587 : // * false_easting
2588 : // * false_northing
2589 : //
2590 : //
2591 : //
2592 : //Grid mapping attributes
2593 : //
2594 : //false_easting
2595 : //false_northing
2596 : //grid_mapping_name
2597 : //grid_north_pole_latitude
2598 : //grid_north_pole_longitude
2599 : //latitude_of_projection_origin
2600 : //longitude_of_central_meridian
2601 : //longitude_of_projection_origin
2602 : //north_pole_grid_longitude
2603 : //scale_factor_at_central_meridian
2604 : //scale_factor_at_projection_origin
2605 : //standard_parallel
2606 : //straight_vertical_longitude_from_pole
2607 :
2608 :
|
