1 : /******************************************************************************
2 : * $Id: netcdfdataset.cpp 23580 2011-12-15 22:24:12Z etourigny $
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 23580 2011-12-15 22:24:12Z etourigny $");
33 :
34 : #include <map> //for NCDFWriteProjAttribs()
35 :
36 : /* Internal function declarations */
37 :
38 : int NCDFIsGDALVersionGTE(const char* pszVersion, int nTarget);
39 :
40 : void NCDFAddGDALHistory( int fpImage,
41 : const char * pszFilename, const char *pszOldHist,
42 : const char * pszFunctionName );
43 :
44 : void NCDFAddHistory(int fpImage, const char *pszAddHist, const char *pszOldHist);
45 :
46 : int NCDFIsCfProjection( const char* pszProjection );
47 :
48 : void NCDFWriteProjAttribs(const OGR_SRSNode *poPROJCS,
49 : const char* pszProjection,
50 : const int fpImage, const int NCDFVarID);
51 :
52 : CPLErr NCDFSafeStrcat(char** ppszDest, char* pszSrc, size_t* nDestSize);
53 :
54 : CPLErr NCDFGetAttr( int nCdfId, int nVarId, const char *pszAttrName,
55 : double *pdfValue );
56 :
57 : CPLErr NCDFGetAttr( int nCdfId, int nVarId, const char *pszAttrName,
58 : char **pszValue );
59 :
60 : CPLErr NCDFPutAttr( int nCdfId, int nVarId,
61 : const char *pszAttrName, const char *pszValue );
62 :
63 :
64 : /************************************************************************/
65 : /* ==================================================================== */
66 : /* netCDFRasterBand */
67 : /* ==================================================================== */
68 : /************************************************************************/
69 :
70 : class netCDFRasterBand : public GDALPamRasterBand
71 : {
72 : friend class netCDFDataset;
73 :
74 : nc_type nc_datatype;
75 : int cdfid;
76 : int nZId;
77 : int nZDim;
78 : int nLevel;
79 : int nBandXPos;
80 : int nBandYPos;
81 : int *panBandZPos;
82 : int *panBandZLev;
83 : int bNoDataSet;
84 : double dfNoDataValue;
85 : double adfValidRange[2];
86 : double dfScale;
87 : double dfOffset;
88 : int bSignedData;
89 : int status;
90 :
91 : CPLErr CreateBandMetadata( int *paDimIds );
92 : template <class T> void CheckValidData ( void * pImage,
93 : int bCheckIsNan=FALSE ) ;
94 :
95 : public:
96 :
97 : netCDFRasterBand( netCDFDataset *poDS,
98 : int nZId,
99 : int nZDim,
100 : int nLevel,
101 : int *panBandZLen,
102 : int *panBandPos,
103 : int *paDimIds,
104 : int nBand );
105 : netCDFRasterBand( netCDFDataset *poDS,
106 : GDALDataType eType,
107 : int nBand,
108 : int bSigned=TRUE,
109 : char *pszBandName=NULL,
110 : char *pszLongName=NULL );
111 : ~netCDFRasterBand( );
112 :
113 : virtual double GetNoDataValue( int * );
114 : virtual CPLErr SetNoDataValue( double );
115 : virtual double GetOffset( int * );
116 : virtual CPLErr SetOffset( double );
117 : virtual double GetScale( int * );
118 : virtual CPLErr SetScale( double );
119 : virtual CPLErr IReadBlock( int, int, void * );
120 : virtual CPLErr IWriteBlock( int, int, void * );
121 :
122 : };
123 :
124 : /************************************************************************/
125 : /* netCDFRasterBand() */
126 : /************************************************************************/
127 :
128 180 : netCDFRasterBand::netCDFRasterBand( netCDFDataset *poNCDFDS,
129 : int nZId,
130 : int nZDim,
131 : int nLevel,
132 : int *panBandZLev,
133 : int *panBandDimPos,
134 : int *paDimIds,
135 180 : int nBand)
136 :
137 : {
138 180 : double dfNoData = 0.0;
139 180 : int bNoDataSet = FALSE;
140 180 : nc_type vartype=NC_NAT;
141 180 : nc_type atttype=NC_NAT;
142 : size_t attlen;
143 : char szNoValueName[NCDF_MAX_STR_LEN];
144 :
145 180 : this->poDS = poNCDFDS;
146 180 : this->panBandZPos = NULL;
147 180 : this->panBandZLev = NULL;
148 180 : this->nBand = nBand;
149 180 : this->nZId = nZId;
150 180 : this->nZDim = nZDim;
151 180 : this->nLevel = nLevel;
152 180 : this->nBandXPos = panBandDimPos[0];
153 180 : this->nBandYPos = panBandDimPos[1];
154 180 : this->bSignedData = TRUE; //default signed, except for Byte
155 180 : this->cdfid = poNCDFDS->GetCDFID();
156 180 : this->status = NC_NOERR;
157 :
158 : /* -------------------------------------------------------------------- */
159 : /* Take care of all other dimmensions */
160 : /* ------------------------------------------------------------------ */
161 180 : if( nZDim > 2 ) {
162 : this->panBandZPos =
163 24 : (int *) CPLCalloc( nZDim-1, sizeof( int ) );
164 : this->panBandZLev =
165 24 : (int *) CPLCalloc( nZDim-1, sizeof( int ) );
166 :
167 80 : for ( int i=0; i < nZDim - 2; i++ ){
168 56 : this->panBandZPos[i] = panBandDimPos[i+2];
169 56 : this->panBandZLev[i] = panBandZLev[i];
170 : }
171 : }
172 180 : bNoDataSet = FALSE;
173 180 : dfNoDataValue = 0.0;
174 :
175 180 : nRasterXSize = poDS->GetRasterXSize( );
176 180 : nRasterYSize = poDS->GetRasterYSize( );
177 180 : nBlockXSize = poDS->GetRasterXSize( );
178 180 : nBlockYSize = 1;
179 :
180 : /* -------------------------------------------------------------------- */
181 : /* Get the type of the "z" variable, our target raster array. */
182 : /* -------------------------------------------------------------------- */
183 180 : if( nc_inq_var( cdfid, nZId, NULL, &nc_datatype, NULL, NULL,
184 : NULL ) != NC_NOERR ){
185 : CPLError( CE_Failure, CPLE_AppDefined,
186 0 : "Error in nc_var_inq() on 'z'." );
187 0 : return;
188 : }
189 :
190 180 : if( nc_datatype == NC_BYTE )
191 93 : eDataType = GDT_Byte;
192 : #ifdef NETCDF_HAS_NC4
193 : /* NC_UBYTE (unsigned byte) is only available for NC4 */
194 : else if( nc_datatype == NC_UBYTE )
195 : eDataType = GDT_Byte;
196 : #endif
197 87 : else if( nc_datatype == NC_CHAR )
198 0 : eDataType = GDT_Byte;
199 87 : else if( nc_datatype == NC_SHORT )
200 9 : eDataType = GDT_Int16;
201 78 : else if( nc_datatype == NC_INT )
202 9 : eDataType = GDT_Int32;
203 69 : else if( nc_datatype == NC_FLOAT )
204 44 : eDataType = GDT_Float32;
205 25 : else if( nc_datatype == NC_DOUBLE )
206 25 : eDataType = GDT_Float64;
207 : else
208 : {
209 0 : if( nBand == 1 )
210 : CPLError( CE_Warning, CPLE_AppDefined,
211 : "Unsupported netCDF datatype (%d), treat as Float32.",
212 0 : (int) nc_datatype );
213 0 : eDataType = GDT_Float32;
214 : }
215 :
216 : /* -------------------------------------------------------------------- */
217 : /* Find out what is No Data for this variable */
218 : /* -------------------------------------------------------------------- */
219 :
220 : status = nc_inq_att( cdfid, nZId,
221 180 : _FillValue, &atttype, &attlen);
222 :
223 : /* -------------------------------------------------------------------- */
224 : /* Look for either Missing_Value or _FillValue attributes */
225 : /* -------------------------------------------------------------------- */
226 :
227 180 : if( status == NC_NOERR ) {
228 153 : strcpy(szNoValueName, _FillValue );
229 : }
230 : else {
231 : status = nc_inq_att( cdfid, nZId,
232 27 : "missing_value", &atttype, &attlen );
233 27 : if( status == NC_NOERR ) {
234 7 : strcpy( szNoValueName, "missing_value" );
235 : }
236 : }
237 :
238 180 : if( status == NC_NOERR ) {
239 160 : if ( NCDFGetAttr( cdfid, nZId, szNoValueName,
240 : &dfNoData ) == CE_None )
241 160 : bNoDataSet = TRUE;
242 : }
243 :
244 180 : nc_inq_vartype( cdfid, nZId, &vartype );
245 :
246 : /* if not found NoData, set the default one */
247 180 : if ( ! bNoDataSet ) {
248 20 : switch( vartype ) {
249 : case NC_BYTE:
250 : #ifdef NETCDF_HAS_NC4
251 : case NC_UBYTE:
252 : #endif
253 : /* don't do default fill-values for bytes, too risky */
254 0 : dfNoData = 0.0;
255 : /* should print a warning as users might not be expecting this */
256 : /* CPLError(CE_Warning, 1,"GDAL netCDF driver is setting default NoData value to 0.0 for NC_BYTE data\n"); */
257 0 : break;
258 : case NC_CHAR:
259 0 : dfNoData = NC_FILL_CHAR;
260 0 : break;
261 : case NC_SHORT:
262 0 : dfNoData = NC_FILL_SHORT;
263 0 : break;
264 : case NC_INT:
265 0 : dfNoData = NC_FILL_INT;
266 0 : break;
267 : case NC_FLOAT:
268 4 : dfNoData = NC_FILL_FLOAT;
269 4 : break;
270 : case NC_DOUBLE:
271 16 : dfNoData = NC_FILL_DOUBLE;
272 16 : break;
273 : default:
274 0 : dfNoData = 0.0;
275 : break;
276 : }
277 20 : bNoDataSet = TRUE;
278 : }
279 :
280 180 : if ( bNoDataSet )
281 180 : SetNoDataValue( dfNoData );
282 : else
283 0 : CPLDebug( "GDAL_netCDF", "did not get nodata value for variable #%d", nZId );
284 :
285 : /* -------------------------------------------------------------------- */
286 : /* Look for valid_range or valid_min/valid_max */
287 : /* -------------------------------------------------------------------- */
288 : /* set valid_range to nodata, then check for actual values */
289 180 : adfValidRange[0] = dfNoData;
290 180 : adfValidRange[1] = dfNoData;
291 : /* first look for valid_range */
292 180 : int bGotValidRange = FALSE;
293 : status = nc_inq_att( cdfid, nZId,
294 180 : "valid_range", &atttype, &attlen);
295 180 : if( (status == NC_NOERR) && (attlen == 2)) {
296 : int vrange[2];
297 : int vmin, vmax;
298 : status = nc_get_att_int( cdfid, nZId,
299 91 : "valid_range", vrange );
300 91 : if( status == NC_NOERR ) {
301 91 : bGotValidRange = TRUE;
302 91 : adfValidRange[0] = vrange[0];
303 91 : adfValidRange[1] = vrange[1];
304 : }
305 : /* if not found look for valid_min and valid_max */
306 : else {
307 : status = nc_get_att_int( cdfid, nZId,
308 0 : "valid_min", &vmin );
309 0 : if( status == NC_NOERR ) {
310 0 : adfValidRange[0] = vmin;
311 : status = nc_get_att_int( cdfid, nZId,
312 0 : "valid_max", &vmax );
313 0 : if( status == NC_NOERR ) {
314 0 : adfValidRange[1] = vmax;
315 0 : bGotValidRange = TRUE;
316 : }
317 : }
318 : }
319 : }
320 :
321 : CPLDebug( "GDAL_netCDF", "valid_range={%f,%f}",
322 180 : adfValidRange[0], adfValidRange[1] );
323 :
324 : /* -------------------------------------------------------------------- */
325 : /* Special For Byte Bands: check for signed/unsigned byte */
326 : /* -------------------------------------------------------------------- */
327 180 : if ( nc_datatype == NC_BYTE ) {
328 :
329 : /* netcdf uses signed byte by default, but GDAL uses unsigned by default */
330 : /* This may cause unexpected results, but is needed for back-compat */
331 93 : if ( poNCDFDS->bIsGdalFile )
332 91 : this->bSignedData = FALSE;
333 : else
334 2 : this->bSignedData = TRUE;
335 :
336 : /* For NC4 format NC_BYTE is signed, NC_UBYTE is unsigned */
337 93 : if ( poNCDFDS->nFormat == NCDF_FORMAT_NC4 ) {
338 0 : this->bSignedData = TRUE;
339 : }
340 : else {
341 : /* if we got valid_range, test for signed/unsigned range */
342 : /* http://www.unidata.ucar.edu/software/netcdf/docs/netcdf/Attribute-Conventions.html */
343 93 : if ( bGotValidRange == TRUE ) {
344 : /* If we got valid_range={0,255}, treat as unsigned */
345 175 : if ( (adfValidRange[0] == 0) && (adfValidRange[1] == 255) ) {
346 84 : bSignedData = FALSE;
347 : /* reset valid_range */
348 84 : adfValidRange[0] = dfNoData;
349 84 : adfValidRange[1] = dfNoData;
350 : }
351 : /* If we got valid_range={-128,127}, treat as signed */
352 7 : else if ( (adfValidRange[0] == -128) && (adfValidRange[1] == 127) ) {
353 7 : bSignedData = TRUE;
354 : /* reset valid_range */
355 7 : adfValidRange[0] = dfNoData;
356 7 : adfValidRange[1] = dfNoData;
357 : }
358 : }
359 : /* else test for _Unsigned */
360 : /* http://www.unidata.ucar.edu/software/netcdf/docs/BestPractices.html */
361 : else {
362 2 : char *pszTemp = NULL;
363 2 : if ( NCDFGetAttr( cdfid, nZId, "_Unsigned", &pszTemp )
364 :
365 : == CE_None ) {
366 2 : if ( EQUAL(pszTemp,"true"))
367 2 : bSignedData = FALSE;
368 0 : else if ( EQUAL(pszTemp,"false"))
369 0 : bSignedData = TRUE;
370 2 : CPLFree( pszTemp );
371 : }
372 : }
373 : }
374 :
375 93 : if ( bSignedData )
376 : {
377 : /* set PIXELTYPE=SIGNEDBYTE */
378 : /* See http://trac.osgeo.org/gdal/wiki/rfc14_imagestructure */
379 7 : SetMetadataItem( "PIXELTYPE", "SIGNEDBYTE", "IMAGE_STRUCTURE" );
380 7 : CPLDebug( "GDAL_netCDF", "got signed Byte" );
381 : }
382 : else
383 86 : CPLDebug( "GDAL_netCDF", "got unsigned Byte" );
384 :
385 : }
386 :
387 : /* -------------------------------------------------------------------- */
388 : /* Create Band Metadata */
389 : /* -------------------------------------------------------------------- */
390 180 : CreateBandMetadata( paDimIds );
391 :
392 : /* -------------------------------------------------------------------- */
393 : /* Attempt to fetch the scale_factor and add_offset attributes for the */
394 : /* variable and set them. If these values are not available, set */
395 : /* offset to 0 and scale to 1 */
396 : /* -------------------------------------------------------------------- */
397 180 : double dfOff = 0.0;
398 180 : double dfScale = 1.0;
399 :
400 180 : if ( nc_inq_attid ( cdfid, nZId, CF_ADD_OFFSET, NULL) == NC_NOERR ) {
401 3 : status = nc_get_att_double( cdfid, nZId, CF_ADD_OFFSET, &dfOff );
402 3 : CPLDebug( "GDAL_netCDF", "got add_offset=%.16g, status=%d", dfOff, status );
403 : }
404 180 : if ( nc_inq_attid ( cdfid, nZId,
405 : CF_SCALE_FACTOR, NULL) == NC_NOERR ) {
406 3 : status = nc_get_att_double( cdfid, nZId, CF_SCALE_FACTOR, &dfScale );
407 3 : CPLDebug( "GDAL_netCDF", "got scale_factor=%.16g, status=%d", dfScale, status );
408 : }
409 180 : SetOffset( dfOff );
410 180 : SetScale( dfScale );
411 0 : }
412 :
413 : /* constructor in create mode, assume just 2 dimensions for now */
414 136 : netCDFRasterBand::netCDFRasterBand( netCDFDataset *poNCDFDS,
415 : GDALDataType eType,
416 : int nBand,
417 : int bSigned,
418 : char *pszBandName,
419 136 : char *pszLongName )
420 :
421 : {
422 : int anBandDims[ NC_MAX_DIMS ];
423 : int status;
424 :
425 136 : double dfNoData = 0.0;
426 : char szTemp[NCDF_MAX_STR_LEN];
427 :
428 136 : this->poDS = poNCDFDS;
429 136 : this->nBand = nBand;
430 136 : this->nZId = -1;
431 136 : this->nZDim = 2;
432 136 : this->nLevel = 0;
433 136 : this->panBandZPos = NULL;
434 136 : this->panBandZLev = NULL;
435 136 : this->nBandXPos = 1;
436 136 : this->nBandYPos = 0;
437 136 : this->bSignedData = bSigned;
438 :
439 136 : this->status = NC_NOERR;
440 136 : this->cdfid = poNCDFDS->GetCDFID();
441 :
442 136 : bNoDataSet = FALSE;
443 136 : dfNoDataValue = 0.0;
444 :
445 136 : nRasterXSize = poDS->GetRasterXSize( );
446 136 : nRasterYSize = poDS->GetRasterYSize( );
447 136 : nBlockXSize = poDS->GetRasterXSize( );
448 136 : nBlockYSize = 1;
449 :
450 136 : if ( poDS->GetAccess() != GA_Update ) {
451 : CPLError( CE_Failure, CPLE_NotSupported,
452 0 : "Dataset is not in update mode, wrong netCDFRasterBand constructor" );
453 0 : return;
454 : }
455 :
456 : /* -------------------------------------------------------------------- */
457 : /* Take care of all other dimmensions */
458 : /* ------------------------------------------------------------------ */
459 : // if( nZDim > 2 ) {
460 : // this->panBandZPos =
461 : // (int *) CPLCalloc( nZDim-1, sizeof( int ) );
462 : // this->panBandZLev =
463 : // (int *) CPLCalloc( nZDim-1, sizeof( int ) );
464 :
465 : // for ( int i=0; i < nZDim - 2; i++ ){
466 : // this->panBandZPos[i] = panBandDimPos[i+2];
467 : // this->panBandZLev[i] = panBandZLev[i];
468 : // }
469 : // }
470 :
471 : /* -------------------------------------------------------------------- */
472 : /* Get the type of the "z" variable, our target raster array. */
473 : /* -------------------------------------------------------------------- */
474 136 : eDataType = eType;
475 :
476 136 : switch ( eDataType )
477 : {
478 : case GDT_Byte:
479 73 : nc_datatype = NC_BYTE;
480 : #ifdef NETCDF_HAS_NC4
481 : /* NC_UBYTE (unsigned byte) is only available for NC4 */
482 : if ( ! bSignedData && (poNCDFDS->nFormat == NCDF_FORMAT_NC4) )
483 : nc_datatype = NC_UBYTE;
484 : #endif
485 73 : break;
486 : case GDT_Int16:
487 8 : nc_datatype = NC_SHORT;
488 8 : break;
489 : case GDT_Int32:
490 8 : nc_datatype = NC_INT;
491 8 : break;
492 : case GDT_Float32:
493 13 : nc_datatype = NC_FLOAT;
494 13 : break;
495 : case GDT_Float64:
496 8 : nc_datatype = NC_DOUBLE;
497 8 : break;
498 : default:
499 26 : if( nBand == 1 )
500 : CPLError( CE_Warning, CPLE_AppDefined,
501 : "Unsupported GDAL datatype (%d), treat as NC_FLOAT.",
502 14 : (int) eDataType );
503 26 : nc_datatype = NC_FLOAT;
504 : break;
505 : }
506 :
507 : /* -------------------------------------------------------------------- */
508 : /* Define the variable */
509 : /* -------------------------------------------------------------------- */
510 136 : anBandDims[0] = poNCDFDS->nYDimID;
511 136 : anBandDims[1] = poNCDFDS->nXDimID;
512 :
513 255 : if ( !pszBandName || EQUAL(pszBandName,"") )
514 119 : sprintf( szTemp, "Band%d", nBand );
515 : else
516 17 : strcpy( szTemp, pszBandName );
517 : CPLDebug( "GDAL_netCDF", "nc_def_var(%d,%s,%d)",
518 136 : cdfid, szTemp, nc_datatype );
519 : status = nc_def_var( cdfid, szTemp, nc_datatype,
520 136 : NCDF_NBDIM, anBandDims, &nZId );
521 136 : NCDF_ERR(status);
522 :
523 256 : if ( !pszLongName || EQUAL(pszLongName,"") )
524 120 : sprintf( szTemp, "GDAL Band Number %d", nBand );
525 : else
526 16 : strcpy( szTemp, pszLongName );
527 : status = nc_put_att_text( cdfid, nZId, CF_LNG_NAME,
528 136 : strlen( szTemp ), szTemp );
529 136 : NCDF_ERR(status);
530 :
531 :
532 136 : poNCDFDS->DefVarDeflate(nZId, TRUE);
533 :
534 : /* for Byte data add signed/unsigned info */
535 136 : if ( eDataType == GDT_Byte ) {
536 :
537 73 : CPLDebug( "GDAL_netCDF", "adding valid_range attributes for Byte Band" );
538 : /* For unsigned NC_BYTE (except NC4 format) */
539 : /* add valid_range and _Unsigned ( defined in CF-1 and NUG ) */
540 73 : if ( (nc_datatype == NC_BYTE) && (poNCDFDS->nFormat != NCDF_FORMAT_NC4) ) {
541 : short int adfValidRange[2];
542 73 : if ( bSignedData ) {
543 3 : adfValidRange[0] = -128;
544 3 : adfValidRange[1] = 127;
545 : status = nc_put_att_text( cdfid,nZId,
546 3 : "_Unsigned", 5, "false" );
547 : }
548 : else {
549 70 : adfValidRange[0] = 0;
550 70 : adfValidRange[1] = 255;
551 : status = nc_put_att_text( cdfid,nZId,
552 70 : "_Unsigned", 4, "true" );
553 : }
554 : status=nc_put_att_short( cdfid,nZId, "valid_range",
555 73 : NC_SHORT, 2, adfValidRange );
556 : }
557 : /* for unsigned byte set PIXELTYPE=SIGNEDBYTE */
558 : /* See http://trac.osgeo.org/gdal/wiki/rfc14_imagestructure */
559 73 : if ( bSignedData )
560 3 : SetMetadataItem( "PIXELTYPE", "SIGNEDBYTE", "IMAGE_STRUCTURE" );
561 :
562 : }
563 :
564 : /* set default nodata */
565 136 : switch( nc_datatype ) {
566 : case NC_BYTE:
567 : #ifdef NETCDF_HAS_NC4
568 : case NC_UBYTE:
569 : #endif
570 : /* don't do default fill-values for bytes, too risky */
571 73 : dfNoData = 0.0;
572 : /* should print a warning as users might not be expecting this */
573 : /* CPLError(CE_Warning, 1,"GDAL netCDF driver is setting default NoData value to 0.0 for NC_BYTE data\n"); */
574 73 : break;
575 : case NC_CHAR:
576 0 : dfNoData = NC_FILL_CHAR;
577 0 : break;
578 : case NC_SHORT:
579 8 : dfNoData = NC_FILL_SHORT;
580 8 : break;
581 : case NC_INT:
582 8 : dfNoData = NC_FILL_INT;
583 8 : break;
584 : case NC_FLOAT:
585 39 : dfNoData = NC_FILL_FLOAT;
586 39 : break;
587 : case NC_DOUBLE:
588 8 : dfNoData = NC_FILL_DOUBLE;
589 8 : break;
590 : default:
591 0 : dfNoData = 0.0;
592 : break;
593 : }
594 :
595 136 : SetNoDataValue( dfNoData );
596 :
597 0 : }
598 :
599 : /************************************************************************/
600 : /* ~netCDFRasterBand() */
601 : /************************************************************************/
602 :
603 316 : netCDFRasterBand::~netCDFRasterBand()
604 : {
605 316 : FlushCache();
606 316 : if( panBandZPos )
607 24 : CPLFree( panBandZPos );
608 316 : if( panBandZLev )
609 24 : CPLFree( panBandZLev );
610 316 : }
611 :
612 : /************************************************************************/
613 : /* GetOffset() */
614 : /************************************************************************/
615 21 : double netCDFRasterBand::GetOffset( int *pbSuccess )
616 : {
617 21 : if( pbSuccess != NULL )
618 21 : *pbSuccess = TRUE;
619 :
620 21 : return dfOffset;
621 : }
622 :
623 : /************************************************************************/
624 : /* SetOffset() */
625 : /************************************************************************/
626 180 : CPLErr netCDFRasterBand::SetOffset( double dfNewOffset )
627 : {
628 180 : dfOffset = dfNewOffset;
629 :
630 : /* write value if in update mode */
631 180 : if ( poDS->GetAccess() == GA_Update ) {
632 :
633 : /* make sure we are in define mode */
634 0 : ( ( netCDFDataset * ) poDS )->SetDefineMode( TRUE );
635 :
636 : status = nc_put_att_double( cdfid, nZId, CF_ADD_OFFSET,
637 0 : NC_DOUBLE, 1, &dfOffset );
638 :
639 0 : NCDF_ERR(status);
640 0 : if ( status == NC_NOERR )
641 0 : return CE_None;
642 : else
643 0 : return CE_Failure;
644 :
645 : }
646 :
647 180 : return CE_None;
648 : }
649 :
650 : /************************************************************************/
651 : /* GetScale() */
652 : /************************************************************************/
653 21 : double netCDFRasterBand::GetScale( int *pbSuccess )
654 : {
655 21 : if( pbSuccess != NULL )
656 21 : *pbSuccess = TRUE;
657 21 : return dfScale;
658 : }
659 :
660 : /************************************************************************/
661 : /* SetScale() */
662 : /************************************************************************/
663 180 : CPLErr netCDFRasterBand::SetScale( double dfNewScale )
664 : {
665 180 : dfScale = dfNewScale;
666 :
667 : /* write value if in update mode */
668 180 : if ( poDS->GetAccess() == GA_Update ) {
669 :
670 : /* make sure we are in define mode */
671 0 : ( ( netCDFDataset * ) poDS )->SetDefineMode( TRUE );
672 :
673 : status = nc_put_att_double( cdfid, nZId, CF_SCALE_FACTOR,
674 0 : NC_DOUBLE, 1, &dfScale );
675 :
676 0 : NCDF_ERR(status);
677 0 : if ( status == NC_NOERR )
678 0 : return CE_None;
679 : else
680 0 : return CE_Failure;
681 :
682 : }
683 :
684 180 : return CE_None;
685 : }
686 :
687 : /************************************************************************/
688 : /* GetNoDataValue() */
689 : /************************************************************************/
690 :
691 121 : double netCDFRasterBand::GetNoDataValue( int * pbSuccess )
692 :
693 : {
694 121 : if( pbSuccess )
695 105 : *pbSuccess = bNoDataSet;
696 :
697 121 : if( bNoDataSet )
698 121 : return dfNoDataValue;
699 : else
700 0 : return GDALPamRasterBand::GetNoDataValue( pbSuccess );
701 : }
702 :
703 : /************************************************************************/
704 : /* SetNoDataValue() */
705 : /************************************************************************/
706 :
707 389 : CPLErr netCDFRasterBand::SetNoDataValue( double dfNoData )
708 :
709 : {
710 389 : bNoDataSet = TRUE;
711 389 : dfNoDataValue = dfNoData;
712 :
713 : /* write value if in update mode */
714 389 : if ( poDS->GetAccess() == GA_Update ) {
715 :
716 : /* make sure we are in define mode */
717 209 : ( ( netCDFDataset * ) poDS )->SetDefineMode( TRUE );
718 :
719 209 : if ( eDataType == GDT_Byte) {
720 127 : if ( bSignedData ) {
721 5 : signed char cNoDataValue = (signed char) dfNoData;
722 : status = nc_put_att_schar( cdfid, nZId, _FillValue,
723 5 : nc_datatype, 1, &cNoDataValue );
724 : }
725 : else {
726 122 : unsigned char ucNoDataValue = (unsigned char) dfNoDataValue;
727 : status = nc_put_att_uchar( cdfid, nZId, _FillValue,
728 122 : nc_datatype, 1, &ucNoDataValue );
729 : }
730 : }
731 82 : else if ( eDataType == GDT_Int16 ) {
732 11 : short int nsNoDataValue = (short int) dfNoDataValue;
733 : status = nc_put_att_short( cdfid, nZId, _FillValue,
734 11 : nc_datatype, 1, &nsNoDataValue );
735 : }
736 71 : else if ( eDataType == GDT_Int32) {
737 11 : int nNoDataValue = (int) dfNoDataValue;
738 : status = nc_put_att_int( cdfid, nZId, _FillValue,
739 11 : nc_datatype, 1, &nNoDataValue );
740 : }
741 60 : else if ( eDataType == GDT_Float32) {
742 21 : float fNoDataValue = (float) dfNoDataValue;
743 : status = nc_put_att_float( cdfid, nZId, _FillValue,
744 21 : nc_datatype, 1, &fNoDataValue );
745 : }
746 : else
747 : status = nc_put_att_double( cdfid, nZId, _FillValue,
748 39 : nc_datatype, 1, &dfNoDataValue );
749 :
750 209 : NCDF_ERR(status);
751 209 : if ( status == NC_NOERR )
752 209 : return CE_None;
753 : else
754 0 : return CE_Failure;
755 :
756 : }
757 180 : return CE_None;
758 : }
759 :
760 : /************************************************************************/
761 : /* CreateBandMetadata() */
762 : /************************************************************************/
763 :
764 180 : CPLErr netCDFRasterBand::CreateBandMetadata( int *paDimIds )
765 :
766 : {
767 : char szVarName[NC_MAX_NAME];
768 : char szMetaName[NC_MAX_NAME];
769 : char szMetaTemp[NCDF_MAX_STR_LEN];
770 180 : char *pszMetaValue = NULL;
771 : char szTemp[NC_MAX_NAME];
772 : const char *pszValue;
773 :
774 : int nd;
775 : int i,j;
776 180 : int Sum = 1;
777 180 : int Taken = 0;
778 180 : int result = 0;
779 : int status;
780 180 : int nVarID = -1;
781 : int nDims;
782 : size_t start[1];
783 : size_t count[1];
784 180 : nc_type nVarType = NC_NAT;
785 180 : int nAtt=0;
786 :
787 180 : netCDFDataset *poDS = (netCDFDataset *) this->poDS;
788 :
789 : /* -------------------------------------------------------------------- */
790 : /* Compute all dimensions from Band number and save in Metadata */
791 : /* -------------------------------------------------------------------- */
792 180 : nc_inq_varname( cdfid, nZId, szVarName );
793 180 : nc_inq_varndims( cdfid, nZId, &nd );
794 : /* -------------------------------------------------------------------- */
795 : /* Compute multidimention band position */
796 : /* */
797 : /* BandPosition = (Total - sum(PastBandLevels) - 1)/sum(remainingLevels)*/
798 : /* if Data[2,3,4,x,y] */
799 : /* */
800 : /* BandPos0 = (nBand ) / (3*4) */
801 : /* BandPos1 = (nBand - BandPos0*(3*4) ) / (4) */
802 : /* BandPos2 = (nBand - BandPos0*(3*4) ) % (4) */
803 : /* -------------------------------------------------------------------- */
804 :
805 180 : sprintf( szMetaName,"NETCDF_VARNAME");
806 180 : sprintf( szMetaTemp,"%s",szVarName);
807 180 : SetMetadataItem( szMetaName, szMetaTemp );
808 180 : if( nd == 3 ) {
809 8 : Sum *= panBandZLev[0];
810 : }
811 :
812 236 : for( i=0; i < nd-2 ; i++ ) {
813 56 : if( i != nd - 2 -1 ) {
814 32 : Sum = 1;
815 80 : for( j=i+1; j < nd-2; j++ ) {
816 48 : Sum *= panBandZLev[j];
817 : }
818 32 : result = (int) ( ( nLevel-Taken) / Sum );
819 : }
820 : else {
821 24 : result = (int) ( ( nLevel-Taken) % Sum );
822 : }
823 :
824 : strcpy(szVarName,
825 56 : poDS->papszDimName[paDimIds[panBandZPos[i]]] );
826 :
827 56 : sprintf( szMetaName,"NETCDF_DIMENSION_%s", szVarName );
828 :
829 56 : status=nc_inq_varid( cdfid, szVarName, &nVarID );
830 :
831 : /* -------------------------------------------------------------------- */
832 : /* Try to uppercase the first letter of the variable */
833 : /* -------------------------------------------------------------------- */
834 :
835 56 : if( status != NC_NOERR ) {
836 48 : szVarName[0]=(char) toupper(szVarName[0]);
837 48 : status=nc_inq_varid( cdfid, szVarName, &nVarID );
838 : }
839 :
840 56 : status = nc_inq_vartype( cdfid, nVarID, &nVarType );
841 :
842 56 : nDims = 0;
843 56 : status = nc_inq_varndims( cdfid, nVarID, &nDims );
844 :
845 56 : if( nDims == 1 ) {
846 8 : count[0]=1;
847 8 : start[0]=result;
848 8 : switch( nVarType ) {
849 : case NC_SHORT:
850 : short sData;
851 : status = nc_get_vara_short( cdfid, nVarID,
852 : start,
853 0 : count, &sData );
854 0 : sprintf( szMetaTemp,"%d", sData );
855 0 : break;
856 : case NC_INT:
857 : int nData;
858 : status = nc_get_vara_int( cdfid, nVarID,
859 : start,
860 0 : count, &nData );
861 0 : sprintf( szMetaTemp,"%d", nData );
862 0 : break;
863 : case NC_FLOAT:
864 : float fData;
865 : status = nc_get_vara_float( cdfid, nVarID,
866 : start,
867 0 : count, &fData );
868 0 : sprintf( szMetaTemp,"%.8g", fData );
869 0 : break;
870 : case NC_DOUBLE:
871 : double dfData;
872 : status = nc_get_vara_double( cdfid, nVarID,
873 : start,
874 8 : count, &dfData);
875 8 : sprintf( szMetaTemp,"%.16g", dfData );
876 8 : break;
877 : default:
878 : CPLDebug( "GDAL_netCDF", "invalid dim %s, type=%d",
879 0 : szMetaTemp, nVarType);
880 : break;
881 : }
882 : }
883 : else
884 48 : sprintf( szMetaTemp,"%d", result+1);
885 :
886 : CPLDebug( "GDAL_netCDF", "setting dimension metadata %s=%s",
887 56 : szMetaName, szMetaTemp );
888 :
889 56 : SetMetadataItem( szMetaName, szMetaTemp );
890 :
891 : /* -------------------------------------------------------------------- */
892 : /* Fetch dimension units */
893 : /* -------------------------------------------------------------------- */
894 :
895 56 : strcpy( szTemp, szVarName );
896 56 : strcat( szTemp, "#units" );
897 56 : pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp);
898 56 : if( pszValue != NULL ) {
899 8 : if( EQUAL( pszValue, "T") ) {
900 0 : strcpy( szTemp, szVarName );
901 0 : strcat( szTemp, "#original_units" );
902 0 : pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp);
903 0 : strcpy( szTemp, "NETCDF_");
904 0 : strcat( szTemp, szVarName );
905 0 : strcat( szTemp, "_original_units" );
906 0 : SetMetadataItem( szTemp, pszValue );
907 : }
908 : else {
909 8 : strcpy( szTemp, "NETCDF_");
910 8 : strcat( szTemp, szVarName );
911 8 : strcat( szTemp, "_units" );
912 8 : SetMetadataItem( szTemp, pszValue );
913 : }
914 : }
915 56 : Taken += result * Sum;
916 : }
917 :
918 : /* -------------------------------------------------------------------- */
919 : /* Get all other metadata */
920 : /* -------------------------------------------------------------------- */
921 180 : nc_inq_varnatts( cdfid, nZId, &nAtt );
922 :
923 995 : for( i=0; i < nAtt ; i++ ) {
924 :
925 815 : status = nc_inq_attname( cdfid, nZId, i, szTemp);
926 : // if(strcmp(szTemp,_FillValue) ==0) continue;
927 815 : sprintf( szMetaName,"%s",szTemp);
928 :
929 815 : if ( NCDFGetAttr( cdfid, nZId, szMetaName, &pszMetaValue)
930 : == CE_None ) {
931 815 : SetMetadataItem( szMetaName, pszMetaValue );
932 : }
933 : else {
934 0 : CPLDebug( "GDAL_netCDF", "invalid Band metadata %s", szMetaName );
935 : }
936 :
937 815 : if ( pszMetaValue ) {
938 815 : CPLFree( pszMetaValue );
939 815 : pszMetaValue = NULL;
940 : }
941 :
942 : }
943 :
944 180 : return CE_None;
945 : }
946 :
947 : /************************************************************************/
948 : /* CheckValidData() */
949 : /************************************************************************/
950 : template <class T>
951 5804 : void netCDFRasterBand::CheckValidData ( void * pImage, int bCheckIsNan )
952 : {
953 : int i;
954 5804 : CPLAssert( pImage != NULL );
955 :
956 : /* check if needed or requested */
957 5804 : if ( (adfValidRange[0] != dfNoDataValue) ||
958 : (adfValidRange[1] != dfNoDataValue) ||
959 : bCheckIsNan ) {
960 :
961 44038 : for( i=0; i<nBlockXSize; i++ ) {
962 : /* check for nodata and nan */
963 43193 : if ( CPLIsEqual( (double) ((T *)pImage)[i], dfNoDataValue ) )
964 0 : continue;
965 43193 : if( bCheckIsNan && CPLIsNan( ( (T *) pImage)[i] ) ) {
966 0 : ( (T *)pImage )[i] = (T)dfNoDataValue;
967 0 : continue;
968 : }
969 : /* check for valid_range */
970 43193 : if ( ( ( adfValidRange[0] != dfNoDataValue ) &&
971 : ( ((T *)pImage)[i] < (T)adfValidRange[0] ) )
972 : ||
973 : ( ( adfValidRange[1] != dfNoDataValue ) &&
974 : ( ((T *)pImage)[i] > (T)adfValidRange[1] ) ) ) {
975 0 : ( (T *)pImage )[i] = (T)dfNoDataValue;
976 : }
977 : }
978 :
979 : }
980 :
981 5804 : }
982 :
983 : /************************************************************************/
984 : /* IReadBlock() */
985 : /************************************************************************/
986 :
987 5804 : CPLErr netCDFRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
988 : void * pImage )
989 :
990 : {
991 : size_t start[ MAX_NC_DIMS ];
992 : size_t edge[ MAX_NC_DIMS ];
993 : char pszName[ NCDF_MAX_STR_LEN ];
994 : int i,j;
995 5804 : int Sum=-1;
996 5804 : int Taken=-1;
997 : int nd;
998 :
999 5804 : *pszName='\0';
1000 5804 : memset( start, 0, sizeof( start ) );
1001 5804 : memset( edge, 0, sizeof( edge ) );
1002 5804 : nc_inq_varndims ( cdfid, nZId, &nd );
1003 :
1004 : /* -------------------------------------------------------------------- */
1005 : /* Locate X, Y and Z position in the array */
1006 : /* -------------------------------------------------------------------- */
1007 :
1008 5804 : start[nBandXPos] = 0; // x dim can move arround in array
1009 : // check y order
1010 5804 : if( ( ( netCDFDataset *) poDS )->bBottomUp ) {
1011 184 : start[nBandYPos] = nRasterYSize - 1 - nBlockYOff;
1012 : } else {
1013 5620 : start[nBandYPos] = nBlockYOff; // y
1014 : }
1015 :
1016 5804 : edge[nBandXPos] = nBlockXSize;
1017 5804 : edge[nBandYPos] = 1;
1018 :
1019 5804 : if( nd == 3 ) {
1020 80 : start[panBandZPos[0]] = nLevel; // z
1021 80 : edge [panBandZPos[0]] = 1;
1022 : }
1023 :
1024 : /* -------------------------------------------------------------------- */
1025 : /* Compute multidimention band position */
1026 : /* */
1027 : /* BandPosition = (Total - sum(PastBandLevels) - 1)/sum(remainingLevels)*/
1028 : /* if Data[2,3,4,x,y] */
1029 : /* */
1030 : /* BandPos0 = (nBand ) / (3*4) */
1031 : /* BandPos1 = (nBand - (3*4) ) / (4) */
1032 : /* BandPos2 = (nBand - (3*4) ) % (4) */
1033 : /* -------------------------------------------------------------------- */
1034 5804 : if (nd > 3)
1035 : {
1036 0 : Taken = 0;
1037 0 : for( i=0; i < nd-2 ; i++ )
1038 : {
1039 0 : if( i != nd - 2 -1 ) {
1040 0 : Sum = 1;
1041 0 : for( j=i+1; j < nd-2; j++ ) {
1042 0 : Sum *= panBandZLev[j];
1043 : }
1044 0 : start[panBandZPos[i]] = (int) ( ( nLevel-Taken) / Sum );
1045 0 : edge[panBandZPos[i]] = 1;
1046 : } else {
1047 0 : start[panBandZPos[i]] = (int) ( ( nLevel-Taken) % Sum );
1048 0 : edge[panBandZPos[i]] = 1;
1049 : }
1050 0 : Taken += start[panBandZPos[i]] * Sum;
1051 : }
1052 : }
1053 :
1054 5804 : if( eDataType == GDT_Byte )
1055 : {
1056 4839 : if (this->bSignedData)
1057 : {
1058 : status = nc_get_vara_schar( cdfid, nZId, start, edge,
1059 60 : (signed char *) pImage );
1060 60 : if ( status == NC_NOERR ) CheckValidData<signed char>( pImage );
1061 : }
1062 : else {
1063 : status = nc_get_vara_uchar( cdfid, nZId, start, edge,
1064 4779 : (unsigned char *) pImage );
1065 4779 : if ( status == NC_NOERR ) CheckValidData<unsigned char>( pImage );
1066 : }
1067 : }
1068 965 : else if( eDataType == GDT_Int16 )
1069 : {
1070 : status = nc_get_vara_short( cdfid, nZId, start, edge,
1071 60 : (short int *) pImage );
1072 60 : if ( status == NC_NOERR ) CheckValidData<short int>( pImage );
1073 : }
1074 905 : else if( eDataType == GDT_Int32 )
1075 : {
1076 : if( sizeof(long) == 4 )
1077 : {
1078 : status = nc_get_vara_long( cdfid, nZId, start, edge,
1079 : (long *) pImage );
1080 : if ( status == NC_NOERR ) CheckValidData<long>( pImage );
1081 : }
1082 : else
1083 : {
1084 : status = nc_get_vara_int( cdfid, nZId, start, edge,
1085 60 : (int *) pImage );
1086 60 : if ( status == NC_NOERR ) CheckValidData<int>( pImage );
1087 : }
1088 : }
1089 845 : else if( eDataType == GDT_Float32 )
1090 : {
1091 : status = nc_get_vara_float( cdfid, nZId, start, edge,
1092 785 : (float *) pImage );
1093 785 : if ( status == NC_NOERR ) CheckValidData<float>( pImage, TRUE );
1094 : }
1095 60 : else if( eDataType == GDT_Float64 )
1096 : {
1097 : status = nc_get_vara_double( cdfid, nZId, start, edge,
1098 60 : (double *) pImage );
1099 60 : if ( status == NC_NOERR ) CheckValidData<double>( pImage, TRUE );
1100 : }
1101 : else
1102 0 : status = NC_EBADTYPE;
1103 :
1104 5804 : if( status != NC_NOERR )
1105 : {
1106 : CPLError( CE_Failure, CPLE_AppDefined,
1107 : "netCDF scanline fetch failed: %s",
1108 0 : nc_strerror( status ) );
1109 0 : return CE_Failure;
1110 : }
1111 : else
1112 5804 : return CE_None;
1113 : }
1114 :
1115 : /************************************************************************/
1116 : /* IWriteBlock() */
1117 : /************************************************************************/
1118 :
1119 3478 : CPLErr netCDFRasterBand::IWriteBlock( int nBlockXOff, int nBlockYOff,
1120 : void * pImage )
1121 :
1122 : {
1123 3478 : int status=-1;//, status2=-1;
1124 : // int cdfid = ( ( netCDFDataset * ) poDS )->cdfid;
1125 : size_t start[ NCDF_NBDIM ];
1126 : size_t edge[ MAX_NC_DIMS ];
1127 : char pszName[ NCDF_MAX_STR_LEN ];
1128 : int i,j;
1129 3478 : int Sum=-1;
1130 3478 : int Taken=-1;
1131 : int nd;
1132 :
1133 3478 : *pszName='\0';
1134 3478 : memset( start, 0, sizeof( start ) );
1135 3478 : memset( edge, 0, sizeof( edge ) );
1136 : // nc_inq_varndims ( cdfid, nZId, &nd );
1137 3478 : nd=2;
1138 :
1139 : /* -------------------------------------------------------------------- */
1140 : /* Locate X, Y and Z position in the array */
1141 : /* -------------------------------------------------------------------- */
1142 :
1143 3478 : start[nBandXPos] = 0; // x dim can move arround in array
1144 : // check y order
1145 3478 : if( ( ( netCDFDataset *) poDS )->bBottomUp ) {
1146 0 : start[nBandYPos] = nRasterYSize - 1 - nBlockYOff;
1147 : } else {
1148 3478 : start[nBandYPos] = nBlockYOff; // y
1149 : }
1150 :
1151 3478 : edge[nBandXPos] = nBlockXSize;
1152 3478 : edge[nBandYPos] = 1;
1153 :
1154 3478 : if( nd == 3 ) {
1155 0 : start[panBandZPos[0]] = nLevel; // z
1156 0 : edge [panBandZPos[0]] = 1;
1157 : }
1158 :
1159 : /* -------------------------------------------------------------------- */
1160 : /* Compute multidimention band position */
1161 : /* */
1162 : /* BandPosition = (Total - sum(PastBandLevels) - 1)/sum(remainingLevels)*/
1163 : /* if Data[2,3,4,x,y] */
1164 : /* */
1165 : /* BandPos0 = (nBand ) / (3*4) */
1166 : /* BandPos1 = (nBand - (3*4) ) / (4) */
1167 : /* BandPos2 = (nBand - (3*4) ) % (4) */
1168 : /* -------------------------------------------------------------------- */
1169 3478 : if (nd > 3)
1170 : {
1171 0 : Taken = 0;
1172 0 : for( i=0; i < nd-2 ; i++ )
1173 : {
1174 0 : if( i != nd - 2 -1 ) {
1175 0 : Sum = 1;
1176 0 : for( j=i+1; j < nd-2; j++ ) {
1177 0 : Sum *= panBandZLev[j];
1178 : }
1179 0 : start[panBandZPos[i]] = (int) ( ( nLevel-Taken) / Sum );
1180 0 : edge[panBandZPos[i]] = 1;
1181 : } else {
1182 0 : start[panBandZPos[i]] = (int) ( ( nLevel-Taken) % Sum );
1183 0 : edge[panBandZPos[i]] = 1;
1184 : }
1185 0 : Taken += start[panBandZPos[i]] * Sum;
1186 : }
1187 : }
1188 :
1189 : /* make sure we are in data mode */
1190 3478 : ( ( netCDFDataset * ) poDS )->SetDefineMode( FALSE );
1191 :
1192 : /* copy data according to type */
1193 3478 : if( eDataType == GDT_Byte ) {
1194 3067 : if ( this->bSignedData )
1195 : status = nc_put_vara_schar( cdfid, nZId, start, edge,
1196 40 : (signed char*) pImage);
1197 : else
1198 : status = nc_put_vara_uchar( cdfid, nZId, start, edge,
1199 3027 : (unsigned char*) pImage);
1200 : }
1201 471 : else if( ( eDataType == GDT_UInt16 ) || ( eDataType == GDT_Int16 ) ) {
1202 : status = nc_put_vara_short( cdfid, nZId, start, edge,
1203 60 : (short int *) pImage);
1204 : }
1205 351 : else if( eDataType == GDT_Int32 ) {
1206 : status = nc_put_vara_int( cdfid, nZId, start, edge,
1207 50 : (int *) pImage);
1208 : }
1209 301 : else if( eDataType == GDT_Float32 ) {
1210 : status = nc_put_vara_float( cdfid, nZId, start, edge,
1211 250 : (float *) pImage);
1212 : }
1213 51 : else if( eDataType == GDT_Float64 ) {
1214 : status = nc_put_vara_double( cdfid, nZId, start, edge,
1215 50 : (double *) pImage);
1216 : }
1217 : else {
1218 : CPLError( CE_Failure, CPLE_NotSupported,
1219 : "The NetCDF driver does not support GDAL data type %d",
1220 1 : eDataType );
1221 1 : status = NC_EBADTYPE;
1222 : }
1223 3478 : NCDF_ERR(status);
1224 :
1225 3478 : if( status != NC_NOERR )
1226 : {
1227 : CPLError( CE_Failure, CPLE_AppDefined,
1228 : "netCDF scanline write failed: %s",
1229 1 : nc_strerror( status ) );
1230 1 : return CE_Failure;
1231 : }
1232 : else
1233 3477 : return CE_None;
1234 :
1235 : }
1236 :
1237 : /************************************************************************/
1238 : /* ==================================================================== */
1239 : /* netCDFDataset */
1240 : /* ==================================================================== */
1241 : /************************************************************************/
1242 :
1243 : /************************************************************************/
1244 : /* netCDFDataset() */
1245 : /************************************************************************/
1246 :
1247 284 : netCDFDataset::netCDFDataset()
1248 :
1249 : {
1250 : /* basic dataset vars */
1251 284 : cdfid = -1;
1252 284 : papszSubDatasets = NULL;
1253 284 : papszMetadata = NULL;
1254 284 : bBottomUp = FALSE;
1255 284 : nFormat = NCDF_FORMAT_NONE;
1256 284 : bIsGdalFile = FALSE;
1257 284 : bIsGdalCfFile = FALSE;
1258 :
1259 : /* projection/GT */
1260 284 : adfGeoTransform[0] = 0.0;
1261 284 : adfGeoTransform[1] = 1.0;
1262 284 : adfGeoTransform[2] = 0.0;
1263 284 : adfGeoTransform[3] = 0.0;
1264 284 : adfGeoTransform[4] = 0.0;
1265 284 : adfGeoTransform[5] = 1.0;
1266 284 : pszProjection = NULL;
1267 284 : nXDimID = -1;
1268 284 : nYDimID = -1;
1269 284 : bIsProjected = FALSE;
1270 284 : bIsGeographic = FALSE; /* can be not projected, and also not geographic */
1271 :
1272 : /* state vars */
1273 284 : status = NC_NOERR;
1274 284 : bDefineMode = TRUE;
1275 284 : bSetProjection = FALSE;
1276 284 : bSetGeoTransform = FALSE;
1277 284 : bAddedProjectionVars = FALSE;
1278 :
1279 : /* create vars */
1280 284 : papszCreationOptions = NULL;
1281 284 : nCompress = NCDF_COMPRESS_NONE;
1282 284 : nZLevel = NCDF_DEFLATE_LEVEL;
1283 284 : nCreateMode = NC_CLOBBER;
1284 284 : bSignedData = TRUE;
1285 284 : }
1286 :
1287 :
1288 : /************************************************************************/
1289 : /* ~netCDFDataset() */
1290 : /************************************************************************/
1291 :
1292 284 : netCDFDataset::~netCDFDataset()
1293 :
1294 : {
1295 : /* make sure projection is written if GeoTransform OR Projection are missing */
1296 284 : if( (GetAccess() == GA_Update) && (! bAddedProjectionVars) ) {
1297 15 : if ( bSetProjection && ! bSetGeoTransform )
1298 0 : AddProjectionVars();
1299 15 : else if ( bSetGeoTransform && ! bSetProjection )
1300 1 : AddProjectionVars();
1301 : // CPLError( CE_Warning, CPLE_AppDefined,
1302 : // "netCDFDataset::~netCDFDataset() Projection was not defined, projection will be missing" );
1303 : }
1304 284 : FlushCache();
1305 :
1306 284 : CSLDestroy( papszMetadata );
1307 284 : CSLDestroy( papszSubDatasets );
1308 284 : CSLDestroy( papszCreationOptions );
1309 :
1310 284 : if( pszProjection )
1311 223 : CPLFree( pszProjection );
1312 :
1313 284 : if( cdfid )
1314 284 : nc_close( cdfid );
1315 284 : }
1316 :
1317 : /************************************************************************/
1318 : /* SetDefineMode() */
1319 : /************************************************************************/
1320 4013 : int netCDFDataset::SetDefineMode( int bNewDefineMode )
1321 : {
1322 : /* do nothing, already in same mode */
1323 4013 : if ( bDefineMode == bNewDefineMode ) return CE_None;
1324 :
1325 : CPLDebug( "GDAL_netCDF", "SetDefineMode(%d) old=%d",
1326 91 : bNewDefineMode, bDefineMode );
1327 :
1328 91 : bDefineMode = bNewDefineMode;
1329 :
1330 91 : if ( bDefineMode == TRUE )
1331 0 : status = nc_redef( cdfid );
1332 : else
1333 91 : status = nc_enddef( cdfid );
1334 :
1335 91 : NCDF_ERR(status);
1336 91 : return status;
1337 : }
1338 :
1339 : /************************************************************************/
1340 : /* GetMetadata() */
1341 : /************************************************************************/
1342 18 : char **netCDFDataset::GetMetadata( const char *pszDomain )
1343 : {
1344 18 : if( pszDomain != NULL && EQUALN( pszDomain, "SUBDATASETS", 11 ) )
1345 0 : return papszSubDatasets;
1346 : else
1347 18 : return GDALDataset::GetMetadata( pszDomain );
1348 : }
1349 :
1350 : /************************************************************************/
1351 : /* GetProjectionRef() */
1352 : /************************************************************************/
1353 :
1354 67 : const char * netCDFDataset::GetProjectionRef()
1355 : {
1356 67 : if( bSetProjection )
1357 67 : return pszProjection;
1358 : else
1359 0 : return GDALPamDataset::GetProjectionRef();
1360 : }
1361 :
1362 : /************************************************************************/
1363 : /* FetchCopyParm() */
1364 : /************************************************************************/
1365 :
1366 991 : double netCDFDataset::FetchCopyParm( const char *pszGridMappingValue,
1367 : const char *pszParm, double dfDefault )
1368 :
1369 : {
1370 : char szTemp[ MAX_NC_NAME ];
1371 : const char *pszValue;
1372 :
1373 991 : strcpy(szTemp,pszGridMappingValue);
1374 991 : strcat( szTemp, "#" );
1375 991 : strcat( szTemp, pszParm );
1376 991 : pszValue = CSLFetchNameValue(papszMetadata, szTemp);
1377 :
1378 991 : if( pszValue )
1379 : {
1380 728 : return CPLAtofM(pszValue);
1381 : }
1382 : else
1383 263 : return dfDefault;
1384 : }
1385 :
1386 : /************************************************************************/
1387 : /* FetchStandardParallels() */
1388 : /************************************************************************/
1389 :
1390 43 : char** netCDFDataset::FetchStandardParallels( const char *pszGridMappingValue )
1391 : {
1392 : char szTemp[ MAX_NC_NAME ];
1393 : const char *pszValue;
1394 43 : char **papszValues = NULL;
1395 : //cf-1.0 tags
1396 43 : strcpy( szTemp,pszGridMappingValue );
1397 43 : strcat( szTemp, "#" );
1398 43 : strcat( szTemp, CF_PP_STD_PARALLEL );
1399 43 : pszValue = CSLFetchNameValue( papszMetadata, szTemp );
1400 43 : if( pszValue != NULL ) {
1401 : // papszValues = CSLTokenizeString2( pszValue, ",", CSLT_STRIPLEADSPACES |
1402 : // CSLT_STRIPENDSPACES );
1403 : /* format has changed to { std1, std2 }, must remove { and } */
1404 37 : strcpy( szTemp,pszValue );
1405 37 : int last_char = strlen(pszValue) - 1;
1406 37 : if ( szTemp[0] == '{' ) szTemp[0] = ' ';
1407 37 : if ( szTemp[last_char] == '}' ) szTemp[last_char] = ' ';
1408 : papszValues = CSLTokenizeString2( szTemp, ",", CSLT_STRIPLEADSPACES |
1409 37 : CSLT_STRIPENDSPACES );
1410 : }
1411 : //try gdal tags
1412 : else
1413 : {
1414 6 : strcpy( szTemp, pszGridMappingValue );
1415 6 : strcat( szTemp, "#" );
1416 6 : strcat( szTemp, CF_PP_STD_PARALLEL_1 );
1417 :
1418 6 : pszValue = CSLFetchNameValue( papszMetadata, szTemp );
1419 :
1420 6 : if ( pszValue != NULL )
1421 0 : papszValues = CSLAddString( papszValues, pszValue );
1422 :
1423 6 : strcpy( szTemp,pszGridMappingValue );
1424 6 : strcat( szTemp, "#" );
1425 6 : strcat( szTemp, CF_PP_STD_PARALLEL_2 );
1426 :
1427 6 : pszValue = CSLFetchNameValue( papszMetadata, szTemp );
1428 :
1429 6 : if( pszValue != NULL )
1430 0 : papszValues = CSLAddString( papszValues, pszValue );
1431 : }
1432 :
1433 43 : return papszValues;
1434 : }
1435 :
1436 : /************************************************************************/
1437 : /* SetProjectionFromVar() */
1438 : /************************************************************************/
1439 166 : void netCDFDataset::SetProjectionFromVar( int var )
1440 : {
1441 : size_t start[2], edge[2];
1442 : unsigned int i;
1443 : const char *pszValue;
1444 : int nVarProjectionID;
1445 : char szVarName[ MAX_NC_NAME ];
1446 : char szTemp[ MAX_NC_NAME ];
1447 : char szGridMappingName[ MAX_NC_NAME ];
1448 : char szGridMappingValue[ MAX_NC_NAME ];
1449 :
1450 166 : double dfStdP1=0.0;
1451 166 : double dfStdP2=0.0;
1452 : double dfCenterLat;
1453 : double dfCenterLon;
1454 : double dfScale;
1455 : double dfFalseEasting;
1456 : double dfFalseNorthing;
1457 : double dfCentralMeridian;
1458 : double dfEarthRadius;
1459 : double dfInverseFlattening;
1460 : double dfLonPrimeMeridian;
1461 : double dfSemiMajorAxis;
1462 : double dfSemiMinorAxis;
1463 :
1464 166 : int bGotGeogCS = FALSE;
1465 166 : int bGotCfSRS = FALSE;
1466 166 : int bGotGdalSRS = FALSE;
1467 166 : int bLookForWellKnownGCS = FALSE; //this could be a Config Option
1468 :
1469 : /* These values from CF metadata */
1470 166 : OGRSpatialReference oSRS;
1471 166 : int nVarDimXID = -1;
1472 166 : int nVarDimYID = -1;
1473 : double *pdfXCoord;
1474 : double *pdfYCoord;
1475 : char szDimNameX[ MAX_NC_NAME ];
1476 : char szDimNameY[ MAX_NC_NAME ];
1477 : int nSpacingBegin;
1478 : int nSpacingMiddle;
1479 : int nSpacingLast;
1480 166 : size_t xdim = nRasterXSize;
1481 166 : size_t ydim = nRasterYSize;
1482 :
1483 166 : const char *pszUnits = NULL;
1484 :
1485 : /* These values from GDAL metadata */
1486 166 : const char *pszWKT = NULL;
1487 166 : const char *pszGeoTransform = NULL;
1488 166 : char **papszGeoTransform = NULL;
1489 :
1490 166 : netCDFDataset * poDS = this; /* perhaps this should be removed for clarity */
1491 :
1492 : /* temp variables to use in SetGeoTransform() and SetProjection() */
1493 : double adfTempGeoTransform[6];
1494 : char *pszTempProjection;
1495 :
1496 166 : int bGotGeoTransform = FALSE;
1497 :
1498 166 : CPLDebug( "GDAL_netCDF", "\n=====\nSetProjectionFromVar( %d )\n", var );
1499 :
1500 : /* -------------------------------------------------------------------- */
1501 : /* Get x/y range information. */
1502 : /* -------------------------------------------------------------------- */
1503 :
1504 166 : adfTempGeoTransform[0] = 0.0;
1505 166 : adfTempGeoTransform[1] = 1.0;
1506 166 : adfTempGeoTransform[2] = 0.0;
1507 166 : adfTempGeoTransform[3] = 0.0;
1508 166 : adfTempGeoTransform[4] = 0.0;
1509 166 : adfTempGeoTransform[5] = 1.0;
1510 166 : pszTempProjection = NULL;
1511 :
1512 : /* -------------------------------------------------------------------- */
1513 : /* Look for grid_mapping metadata */
1514 : /* -------------------------------------------------------------------- */
1515 :
1516 166 : strcpy( szGridMappingValue, "" );
1517 166 : strcpy( szGridMappingName, "" );
1518 :
1519 166 : nc_inq_varname( cdfid, var, szVarName );
1520 166 : strcpy(szTemp,szVarName);
1521 166 : strcat(szTemp,"#");
1522 166 : strcat(szTemp,CF_GRD_MAPPING);
1523 :
1524 166 : pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp);
1525 166 : if( pszValue ) {
1526 116 : strcpy(szGridMappingName,szTemp);
1527 116 : strcpy(szGridMappingValue,pszValue);
1528 : }
1529 :
1530 166 : if( !EQUAL( szGridMappingValue, "" ) ) {
1531 :
1532 : /* Read grid_mapping metadata */
1533 116 : nc_inq_varid( cdfid, szGridMappingValue, &nVarProjectionID );
1534 116 : poDS->ReadAttributes( cdfid, nVarProjectionID );
1535 :
1536 : /* -------------------------------------------------------------------- */
1537 : /* Look for GDAL spatial_ref and GeoTransform within grid_mapping */
1538 : /* -------------------------------------------------------------------- */
1539 116 : CPLDebug( "GDAL_netCDF", "got grid_mapping %s", szGridMappingValue );
1540 116 : strcpy( szTemp,szGridMappingValue);
1541 116 : strcat( szTemp, "#" );
1542 116 : strcat( szTemp, NCDF_SPATIAL_REF);
1543 :
1544 116 : pszWKT = CSLFetchNameValue(poDS->papszMetadata, szTemp);
1545 :
1546 116 : if( pszWKT != NULL ) {
1547 111 : strcpy( szTemp,szGridMappingValue);
1548 111 : strcat( szTemp, "#" );
1549 111 : strcat( szTemp, NCDF_GEOTRANSFORM);
1550 111 : pszGeoTransform = CSLFetchNameValue(poDS->papszMetadata, szTemp);
1551 : }
1552 : }
1553 :
1554 : /* -------------------------------------------------------------------- */
1555 : /* Get information about the file. */
1556 : /* -------------------------------------------------------------------- */
1557 : /* Was this file created by the GDAL netcdf driver? */
1558 : /* Was this file created by the newer (CF-conformant) driver? */
1559 : /* -------------------------------------------------------------------- */
1560 : /* 1) If GDAL netcdf metadata is set, and version >= 1.9, */
1561 : /* it was created with the new driver */
1562 : /* 2) Else, if spatial_ref and GeoTransform are present in the */
1563 : /* grid_mapping variable, it was created by the old driver */
1564 : /* -------------------------------------------------------------------- */
1565 166 : pszValue = CSLFetchNameValue(poDS->papszMetadata, "NC_GLOBAL#GDAL");
1566 :
1567 166 : if( pszValue && NCDFIsGDALVersionGTE(pszValue, 1900)) {
1568 147 : bIsGdalFile = TRUE;
1569 147 : bIsGdalCfFile = TRUE;
1570 : }
1571 19 : else if( pszWKT != NULL && pszGeoTransform != NULL ) {
1572 0 : bIsGdalFile = TRUE;
1573 0 : bIsGdalCfFile = FALSE;
1574 : }
1575 :
1576 : /* -------------------------------------------------------------------- */
1577 : /* Set default bottom-up default value */
1578 : /* Config option GDAL_NETCDF_BOTTOMUP and Y axis dimension variable*/
1579 : /* override the default */
1580 : /* -------------------------------------------------------------------- */
1581 166 : pszValue = CPLGetConfigOption( "GDAL_NETCDF_BOTTOMUP", NULL );
1582 166 : if ( pszValue ) {
1583 0 : poDS->bBottomUp = CSLTestBoolean( pszValue ) != FALSE;
1584 : }
1585 : else {
1586 : // if ( bIsGdalFile && ! bIsGdalCfFile )
1587 : /* new driver is not not bottom by default any more */
1588 166 : if ( bIsGdalFile )
1589 147 : poDS->bBottomUp = FALSE;
1590 : else
1591 19 : poDS->bBottomUp = TRUE;
1592 : }
1593 : CPLDebug( "GDAL_netCDF",
1594 : "bIsGdalFile=%d bIsGdalCfFile=%d bBottomUp=%d",
1595 166 : bIsGdalFile, bIsGdalCfFile, bBottomUp );
1596 :
1597 : /* -------------------------------------------------------------------- */
1598 : /* Look for dimension: lon */
1599 : /* -------------------------------------------------------------------- */
1600 :
1601 166 : memset( szDimNameX, '\0', sizeof( char ) * MAX_NC_NAME );
1602 166 : memset( szDimNameY, '\0', sizeof( char ) * MAX_NC_NAME );
1603 :
1604 416 : for( i = 0; (i < strlen( poDS->papszDimName[ poDS->nXDimID ] ) &&
1605 : i < 3 ); i++ ) {
1606 250 : szDimNameX[i]=(char)tolower( ( poDS->papszDimName[poDS->nXDimID] )[i] );
1607 : }
1608 166 : szDimNameX[3] = '\0';
1609 416 : for( i = 0; (i < strlen( poDS->papszDimName[ poDS->nYDimID ] ) &&
1610 : i < 3 ); i++ ) {
1611 250 : szDimNameY[i]=(char)tolower( ( poDS->papszDimName[poDS->nYDimID] )[i] );
1612 : }
1613 166 : szDimNameY[3] = '\0';
1614 :
1615 : /* -------------------------------------------------------------------- */
1616 : /* Read grid_mapping information and set projections */
1617 : /* -------------------------------------------------------------------- */
1618 :
1619 166 : if( !( EQUAL(szGridMappingName,"" ) ) ) {
1620 :
1621 116 : strcpy( szTemp, szGridMappingValue );
1622 116 : strcat( szTemp, "#" );
1623 116 : strcat( szTemp, CF_GRD_MAPPING_NAME );
1624 116 : pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp);
1625 :
1626 116 : if( pszValue != NULL ) {
1627 :
1628 : /* -------------------------------------------------------------------- */
1629 : /* Check for datum/spheroid information */
1630 : /* -------------------------------------------------------------------- */
1631 : dfEarthRadius =
1632 : poDS->FetchCopyParm( szGridMappingValue,
1633 : CF_PP_EARTH_RADIUS,
1634 116 : -1.0 );
1635 :
1636 : dfLonPrimeMeridian =
1637 : poDS->FetchCopyParm( szGridMappingValue,
1638 : CF_PP_LONG_PRIME_MERIDIAN,
1639 116 : 0.0 );
1640 :
1641 : dfInverseFlattening =
1642 : poDS->FetchCopyParm( szGridMappingValue,
1643 : CF_PP_INVERSE_FLATTENING,
1644 116 : -1.0 );
1645 :
1646 : dfSemiMajorAxis =
1647 : poDS->FetchCopyParm( szGridMappingValue,
1648 : CF_PP_SEMI_MAJOR_AXIS,
1649 116 : -1.0 );
1650 :
1651 : dfSemiMinorAxis =
1652 : poDS->FetchCopyParm( szGridMappingValue,
1653 : CF_PP_SEMI_MINOR_AXIS,
1654 116 : -1.0 );
1655 :
1656 : //see if semi-major exists if radius doesn't
1657 116 : if( dfEarthRadius < 0.0 )
1658 116 : dfEarthRadius = dfSemiMajorAxis;
1659 :
1660 : //if still no radius, check old tag
1661 116 : if( dfEarthRadius < 0.0 )
1662 : dfEarthRadius = poDS->FetchCopyParm( szGridMappingValue,
1663 : CF_PP_EARTH_RADIUS_OLD,
1664 5 : -1.0 );
1665 :
1666 : //has radius value
1667 116 : if( dfEarthRadius > 0.0 ) {
1668 : //check for inv_flat tag
1669 114 : if( dfInverseFlattening < 0.0 ) {
1670 : //no inv_flat tag, check for semi_minor
1671 2 : if( dfSemiMinorAxis < 0.0 ) {
1672 : //no way to get inv_flat, use sphere
1673 : oSRS.SetGeogCS( "unknown",
1674 : NULL,
1675 : "Sphere",
1676 2 : dfEarthRadius, 0.0 );
1677 2 : bGotGeogCS = TRUE;
1678 : }
1679 : else {
1680 0 : if( dfSemiMajorAxis < 0.0 )
1681 0 : dfSemiMajorAxis = dfEarthRadius;
1682 : //set inv_flat using semi_minor/major
1683 : dfInverseFlattening =
1684 0 : 1.0 / ( dfSemiMajorAxis - dfSemiMinorAxis ) / dfSemiMajorAxis;
1685 : oSRS.SetGeogCS( "unknown",
1686 : NULL,
1687 : "Spheroid",
1688 0 : dfEarthRadius, dfInverseFlattening );
1689 0 : bGotGeogCS = TRUE;
1690 : }
1691 : }
1692 : else {
1693 : oSRS.SetGeogCS( "unknown",
1694 : NULL,
1695 : "Spheroid",
1696 112 : dfEarthRadius, dfInverseFlattening );
1697 112 : bGotGeogCS = TRUE;
1698 : }
1699 :
1700 114 : if ( bGotGeogCS )
1701 114 : CPLDebug( "GDAL_netCDF", "got spheroid from CF: (%f , %f)", dfEarthRadius, dfInverseFlattening );
1702 :
1703 : }
1704 : //no radius, set as wgs84 as default?
1705 : else {
1706 : // This would be too indiscrimant. But we should set
1707 : // it if we know the data is geographic.
1708 : //oSRS.SetWellKnownGeogCS( "WGS84" );
1709 : }
1710 :
1711 : /* -------------------------------------------------------------------- */
1712 : /* Transverse Mercator */
1713 : /* -------------------------------------------------------------------- */
1714 :
1715 116 : if( EQUAL( pszValue, CF_PT_TM ) ) {
1716 :
1717 : dfScale =
1718 : poDS->FetchCopyParm( szGridMappingValue,
1719 26 : CF_PP_SCALE_FACTOR_MERIDIAN, 1.0 );
1720 :
1721 : dfCenterLon =
1722 : poDS->FetchCopyParm( szGridMappingValue,
1723 26 : CF_PP_LONG_CENTRAL_MERIDIAN, 0.0 );
1724 :
1725 : dfCenterLat =
1726 : poDS->FetchCopyParm( szGridMappingValue,
1727 26 : CF_PP_LAT_PROJ_ORIGIN, 0.0 );
1728 :
1729 : dfFalseEasting =
1730 : poDS->FetchCopyParm( szGridMappingValue,
1731 26 : CF_PP_FALSE_EASTING, 0.0 );
1732 :
1733 : dfFalseNorthing =
1734 : poDS->FetchCopyParm( szGridMappingValue,
1735 26 : CF_PP_FALSE_NORTHING, 0.0 );
1736 :
1737 26 : bGotCfSRS = TRUE;
1738 : oSRS.SetTM( dfCenterLat,
1739 : dfCenterLon,
1740 : dfScale,
1741 : dfFalseEasting,
1742 26 : dfFalseNorthing );
1743 :
1744 26 : if( !bGotGeogCS )
1745 0 : oSRS.SetWellKnownGeogCS( "WGS84" );
1746 : }
1747 :
1748 : /* -------------------------------------------------------------------- */
1749 : /* Albers Equal Area */
1750 : /* -------------------------------------------------------------------- */
1751 :
1752 116 : if( EQUAL( pszValue, CF_PT_AEA ) ) {
1753 7 : char **papszStdParallels = NULL;
1754 :
1755 : dfCenterLon =
1756 : poDS->FetchCopyParm( szGridMappingValue,
1757 7 : CF_PP_LONG_CENTRAL_MERIDIAN, 0.0 );
1758 :
1759 : dfCenterLat =
1760 : poDS->FetchCopyParm( szGridMappingValue,
1761 7 : CF_PP_LAT_PROJ_ORIGIN, 0.0 );
1762 :
1763 : dfFalseEasting =
1764 : poDS->FetchCopyParm( szGridMappingValue,
1765 7 : CF_PP_FALSE_EASTING, 0.0 );
1766 :
1767 : dfFalseNorthing =
1768 : poDS->FetchCopyParm( szGridMappingValue,
1769 7 : CF_PP_FALSE_NORTHING, 0.0 );
1770 :
1771 : papszStdParallels =
1772 7 : FetchStandardParallels( szGridMappingValue );
1773 :
1774 7 : if( papszStdParallels != NULL ) {
1775 :
1776 7 : if ( CSLCount( papszStdParallels ) == 1 ) {
1777 : /* TODO CF-1 standard says it allows AEA to be encoded with only 1 standard parallel */
1778 : /* how should this actually map to a 2StdP OGC WKT version? */
1779 : CPLError( CE_Warning, CPLE_NotSupported,
1780 0 : "NetCDF driver import of AEA-1SP is not tested, using identical std. parallels\n" );
1781 0 : dfStdP1 = CPLAtofM( papszStdParallels[0] );
1782 0 : dfStdP2 = dfStdP1;
1783 :
1784 : }
1785 :
1786 7 : else if( CSLCount( papszStdParallels ) == 2 ) {
1787 7 : dfStdP1 = CPLAtofM( papszStdParallels[0] );
1788 7 : dfStdP2 = CPLAtofM( papszStdParallels[1] );
1789 : }
1790 : }
1791 : //old default
1792 : else {
1793 : dfStdP1 =
1794 : poDS->FetchCopyParm( szGridMappingValue,
1795 0 : CF_PP_STD_PARALLEL_1, 0.0 );
1796 :
1797 : dfStdP2 =
1798 : poDS->FetchCopyParm( szGridMappingValue,
1799 0 : CF_PP_STD_PARALLEL_2, 0.0 );
1800 : }
1801 :
1802 : dfCenterLat =
1803 : poDS->FetchCopyParm( szGridMappingValue,
1804 7 : CF_PP_LAT_PROJ_ORIGIN, 0.0 );
1805 :
1806 7 : bGotCfSRS = TRUE;
1807 : oSRS.SetACEA( dfStdP1, dfStdP2, dfCenterLat, dfCenterLon,
1808 7 : dfFalseEasting, dfFalseNorthing );
1809 :
1810 7 : if( !bGotGeogCS )
1811 0 : oSRS.SetWellKnownGeogCS( "WGS84" );
1812 :
1813 7 : CSLDestroy( papszStdParallels );
1814 : }
1815 :
1816 : /* -------------------------------------------------------------------- */
1817 : /* Cylindrical Equal Area */
1818 : /* -------------------------------------------------------------------- */
1819 :
1820 115 : else if( EQUAL( pszValue, CF_PT_CEA ) || EQUAL( pszValue, CF_PT_LCEA ) ) {
1821 6 : char **papszStdParallels = NULL;
1822 :
1823 : papszStdParallels =
1824 6 : FetchStandardParallels( szGridMappingValue );
1825 :
1826 6 : if( papszStdParallels != NULL ) {
1827 6 : dfStdP1 = CPLAtofM( papszStdParallels[0] );
1828 : }
1829 : else {
1830 : //TODO: add support for 'scale_factor_at_projection_origin' variant to standard parallel
1831 : //Probably then need to calc a std parallel equivalent
1832 : CPLError( CE_Failure, CPLE_NotSupported,
1833 : "NetCDF driver does not support import of CF-1 LCEA "
1834 0 : "'scale_factor_at_projection_origin' variant yet.\n" );
1835 : }
1836 :
1837 : dfCentralMeridian =
1838 : poDS->FetchCopyParm( szGridMappingValue,
1839 6 : CF_PP_LONG_CENTRAL_MERIDIAN, 0.0 );
1840 :
1841 : dfFalseEasting =
1842 : poDS->FetchCopyParm( szGridMappingValue,
1843 6 : CF_PP_FALSE_EASTING, 0.0 );
1844 :
1845 : dfFalseNorthing =
1846 : poDS->FetchCopyParm( szGridMappingValue,
1847 6 : CF_PP_FALSE_NORTHING, 0.0 );
1848 :
1849 6 : bGotCfSRS = TRUE;
1850 : oSRS.SetCEA( dfStdP1, dfCentralMeridian,
1851 6 : dfFalseEasting, dfFalseNorthing );
1852 :
1853 6 : if( !bGotGeogCS )
1854 0 : oSRS.SetWellKnownGeogCS( "WGS84" );
1855 :
1856 : }
1857 :
1858 : /* -------------------------------------------------------------------- */
1859 : /* lambert_azimuthal_equal_area */
1860 : /* -------------------------------------------------------------------- */
1861 103 : else if( EQUAL( pszValue, CF_PT_LAEA ) ) {
1862 : dfCenterLon =
1863 : poDS->FetchCopyParm( szGridMappingValue,
1864 6 : CF_PP_LON_PROJ_ORIGIN, 0.0 );
1865 :
1866 : dfCenterLat =
1867 : poDS->FetchCopyParm( szGridMappingValue,
1868 6 : CF_PP_LAT_PROJ_ORIGIN, 0.0 );
1869 :
1870 : dfFalseEasting =
1871 : poDS->FetchCopyParm( szGridMappingValue,
1872 6 : CF_PP_FALSE_EASTING, 0.0 );
1873 :
1874 : dfFalseNorthing =
1875 : poDS->FetchCopyParm( szGridMappingValue,
1876 6 : CF_PP_FALSE_NORTHING, 0.0 );
1877 :
1878 6 : oSRS.SetProjCS( "LAEA (WGS84) " );
1879 :
1880 6 : bGotCfSRS = TRUE;
1881 : oSRS.SetLAEA( dfCenterLat, dfCenterLon,
1882 6 : dfFalseEasting, dfFalseNorthing );
1883 :
1884 6 : if( !bGotGeogCS )
1885 0 : oSRS.SetWellKnownGeogCS( "WGS84" );
1886 :
1887 : }
1888 :
1889 : /* -------------------------------------------------------------------- */
1890 : /* Azimuthal Equidistant */
1891 : /* -------------------------------------------------------------------- */
1892 97 : else if( EQUAL( pszValue, CF_PT_AE ) ) {
1893 : dfCenterLon =
1894 : poDS->FetchCopyParm( szGridMappingValue,
1895 6 : CF_PP_LON_PROJ_ORIGIN, 0.0 );
1896 :
1897 : dfCenterLat =
1898 : poDS->FetchCopyParm( szGridMappingValue,
1899 6 : CF_PP_LAT_PROJ_ORIGIN, 0.0 );
1900 :
1901 : dfFalseEasting =
1902 : poDS->FetchCopyParm( szGridMappingValue,
1903 6 : CF_PP_FALSE_EASTING, 0.0 );
1904 :
1905 : dfFalseNorthing =
1906 : poDS->FetchCopyParm( szGridMappingValue,
1907 6 : CF_PP_FALSE_NORTHING, 0.0 );
1908 :
1909 6 : bGotCfSRS = TRUE;
1910 : oSRS.SetAE( dfCenterLat, dfCenterLon,
1911 6 : dfFalseEasting, dfFalseNorthing );
1912 :
1913 6 : if( !bGotGeogCS )
1914 0 : oSRS.SetWellKnownGeogCS( "WGS84" );
1915 :
1916 : }
1917 :
1918 : /* -------------------------------------------------------------------- */
1919 : /* Lambert conformal conic */
1920 : /* -------------------------------------------------------------------- */
1921 91 : else if( EQUAL( pszValue, CF_PT_LCC ) ) {
1922 :
1923 12 : char **papszStdParallels = NULL;
1924 :
1925 : dfCenterLon =
1926 : poDS->FetchCopyParm( szGridMappingValue,
1927 12 : CF_PP_LONG_CENTRAL_MERIDIAN, 0.0 );
1928 :
1929 : dfCenterLat =
1930 : poDS->FetchCopyParm( szGridMappingValue,
1931 12 : CF_PP_LAT_PROJ_ORIGIN, 0.0 );
1932 :
1933 : dfFalseEasting =
1934 : poDS->FetchCopyParm( szGridMappingValue,
1935 12 : CF_PP_FALSE_EASTING, 0.0 );
1936 :
1937 : dfFalseNorthing =
1938 : poDS->FetchCopyParm( szGridMappingValue,
1939 12 : CF_PP_FALSE_NORTHING, 0.0 );
1940 :
1941 : papszStdParallels =
1942 12 : FetchStandardParallels( szGridMappingValue );
1943 :
1944 : /* 2SP variant */
1945 12 : if( CSLCount( papszStdParallels ) == 2 ) {
1946 11 : dfStdP1 = CPLAtofM( papszStdParallels[0] );
1947 11 : dfStdP2 = CPLAtofM( papszStdParallels[1] );
1948 : oSRS.SetLCC( dfStdP1, dfStdP2, dfCenterLat, dfCenterLon,
1949 11 : dfFalseEasting, dfFalseNorthing );
1950 : }
1951 : /* 1SP variant (with standard_parallel or center lon) */
1952 : /* See comments in netcdfdataset.h for this projection. */
1953 : else {
1954 :
1955 : dfScale =
1956 : poDS->FetchCopyParm( szGridMappingValue,
1957 1 : CF_PP_SCALE_FACTOR_ORIGIN, -1.0 );
1958 :
1959 : /* CF definition, without scale factor */
1960 1 : if( CPLIsEqual(dfScale, -1.0) ) {
1961 :
1962 : /* with standard_parallel */
1963 1 : if( CSLCount( papszStdParallels ) == 1 )
1964 1 : dfStdP1 = CPLAtofM( papszStdParallels[0] );
1965 : /* with center lon instead */
1966 : else
1967 0 : dfStdP1 = dfCenterLat;
1968 1 : dfStdP2 = dfStdP1;
1969 :
1970 : /* test if we should actually compute scale factor */
1971 1 : if ( ! CPLIsEqual( dfStdP1, dfCenterLat ) ) {
1972 : CPLError( CE_Warning, CPLE_NotSupported,
1973 : "NetCDF driver import of LCC-1SP with standard_parallel1 != latitude_of_projection_origin\n"
1974 0 : "(which forces a computation of scale_factor) is experimental (bug #3324)\n" );
1975 : /* use Snyder eq. 15-4 to compute dfScale from dfStdP1 and dfCenterLat */
1976 : /* only tested for dfStdP1=dfCenterLat and (25,26), needs more data for testing */
1977 : /* other option: use the 2SP variant - how to compute new standard parallels? */
1978 : #define NCDF_PI 3.14159265358979323846
1979 : dfScale = ( cos(dfStdP1) * pow( tan(NCDF_PI/4 + dfStdP1/2), sin(dfStdP1) ) ) /
1980 0 : ( cos(dfCenterLat) * pow( tan(NCDF_PI/4 + dfCenterLat/2), sin(dfCenterLat) ) );
1981 : }
1982 : /* default is 1.0 */
1983 : else
1984 1 : dfScale = 1.0;
1985 :
1986 : oSRS.SetLCC1SP( dfCenterLat, dfCenterLon, dfScale,
1987 1 : dfFalseEasting, dfFalseNorthing );
1988 : /* store dfStdP1 so we can output it to CF later */
1989 1 : oSRS.SetNormProjParm( SRS_PP_STANDARD_PARALLEL_1, dfStdP1 );
1990 : }
1991 : /* OGC/PROJ.4 definition with scale factor */
1992 : else {
1993 : oSRS.SetLCC1SP( dfCenterLat, dfCenterLon, dfScale,
1994 0 : dfFalseEasting, dfFalseNorthing );
1995 : }
1996 : }
1997 :
1998 :
1999 12 : bGotCfSRS = TRUE;
2000 12 : if( !bGotGeogCS )
2001 2 : oSRS.SetWellKnownGeogCS( "WGS84" );
2002 :
2003 12 : CSLDestroy( papszStdParallels );
2004 : }
2005 :
2006 : /* -------------------------------------------------------------------- */
2007 : /* Is this Latitude/Longitude Grid explicitly */
2008 : /* -------------------------------------------------------------------- */
2009 :
2010 79 : else if ( EQUAL ( pszValue, CF_PT_LATITUDE_LONGITUDE ) ) {
2011 23 : bGotCfSRS = TRUE;
2012 23 : if( !bGotGeogCS )
2013 0 : oSRS.SetWellKnownGeogCS( "WGS84" );
2014 : }
2015 : /* -------------------------------------------------------------------- */
2016 : /* Mercator */
2017 : /* -------------------------------------------------------------------- */
2018 :
2019 56 : else if ( EQUAL ( pszValue, CF_PT_MERCATOR ) ) {
2020 12 : char **papszStdParallels = NULL;
2021 :
2022 : /* If there is a standard_parallel, know it is Mercator 2SP */
2023 : papszStdParallels =
2024 12 : FetchStandardParallels( szGridMappingValue );
2025 :
2026 12 : if (NULL != papszStdParallels) {
2027 : /* CF-1 Mercator 2SP always has lat centered at equator */
2028 6 : dfStdP1 = CPLAtofM( papszStdParallels[0] );
2029 :
2030 6 : dfCenterLat = 0.0;
2031 :
2032 : dfCenterLon =
2033 : poDS->FetchCopyParm( szGridMappingValue,
2034 6 : CF_PP_LON_PROJ_ORIGIN, 0.0 );
2035 :
2036 : dfFalseEasting =
2037 : poDS->FetchCopyParm( szGridMappingValue,
2038 6 : CF_PP_FALSE_EASTING, 0.0 );
2039 :
2040 : dfFalseNorthing =
2041 : poDS->FetchCopyParm( szGridMappingValue,
2042 6 : CF_PP_FALSE_NORTHING, 0.0 );
2043 :
2044 : oSRS.SetMercator2SP( dfStdP1, dfCenterLat, dfCenterLon,
2045 6 : dfFalseEasting, dfFalseNorthing );
2046 : }
2047 : else {
2048 : dfCenterLon =
2049 : poDS->FetchCopyParm( szGridMappingValue,
2050 6 : CF_PP_LON_PROJ_ORIGIN, 0.0 );
2051 :
2052 : dfCenterLat =
2053 : poDS->FetchCopyParm( szGridMappingValue,
2054 6 : CF_PP_LAT_PROJ_ORIGIN, 0.0 );
2055 :
2056 : dfScale =
2057 : poDS->FetchCopyParm( szGridMappingValue,
2058 : CF_PP_SCALE_FACTOR_ORIGIN,
2059 6 : 1.0 );
2060 :
2061 : dfFalseEasting =
2062 : poDS->FetchCopyParm( szGridMappingValue,
2063 6 : CF_PP_FALSE_EASTING, 0.0 );
2064 :
2065 : dfFalseNorthing =
2066 : poDS->FetchCopyParm( szGridMappingValue,
2067 6 : CF_PP_FALSE_NORTHING, 0.0 );
2068 :
2069 : oSRS.SetMercator( dfCenterLat, dfCenterLon, dfScale,
2070 6 : dfFalseEasting, dfFalseNorthing );
2071 : }
2072 :
2073 12 : bGotCfSRS = TRUE;
2074 :
2075 12 : if( !bGotGeogCS )
2076 0 : oSRS.SetWellKnownGeogCS( "WGS84" );
2077 : }
2078 :
2079 : /* -------------------------------------------------------------------- */
2080 : /* Orthographic */
2081 : /* -------------------------------------------------------------------- */
2082 :
2083 :
2084 44 : else if ( EQUAL ( pszValue, CF_PT_ORTHOGRAPHIC ) ) {
2085 : dfCenterLon =
2086 : poDS->FetchCopyParm( szGridMappingValue,
2087 6 : CF_PP_LON_PROJ_ORIGIN, 0.0 );
2088 :
2089 : dfCenterLat =
2090 : poDS->FetchCopyParm( szGridMappingValue,
2091 6 : CF_PP_LAT_PROJ_ORIGIN, 0.0 );
2092 :
2093 : dfFalseEasting =
2094 : poDS->FetchCopyParm( szGridMappingValue,
2095 6 : CF_PP_FALSE_EASTING, 0.0 );
2096 :
2097 : dfFalseNorthing =
2098 : poDS->FetchCopyParm( szGridMappingValue,
2099 6 : CF_PP_FALSE_NORTHING, 0.0 );
2100 :
2101 6 : bGotCfSRS = TRUE;
2102 :
2103 : oSRS.SetOrthographic( dfCenterLat, dfCenterLon,
2104 6 : dfFalseEasting, dfFalseNorthing );
2105 :
2106 6 : if( !bGotGeogCS )
2107 0 : oSRS.SetWellKnownGeogCS( "WGS84" );
2108 : }
2109 :
2110 : /* -------------------------------------------------------------------- */
2111 : /* Polar Stereographic */
2112 : /* -------------------------------------------------------------------- */
2113 :
2114 38 : else if ( EQUAL ( pszValue, CF_PT_POLAR_STEREO ) ) {
2115 : /* Note: reversing our current mapping on export, from
2116 : 'latitude_of_origin' in OGC WKT to 'standard_parallel' in CF-1
2117 : See comments in netcdfdataset.h for this projection. */
2118 6 : char **papszStdParallels = NULL;
2119 :
2120 : papszStdParallels =
2121 6 : FetchStandardParallels( szGridMappingValue );
2122 :
2123 6 : if (NULL != papszStdParallels) {
2124 6 : dfCenterLat = CPLAtofM( papszStdParallels[0] );
2125 : }
2126 : else {
2127 : //TODO: not sure how to handle if CF-1 doesn't include a std_parallel:
2128 0 : dfCenterLat = 0.0; //just to avoid warning at compilation
2129 : CPLError( CE_Failure, CPLE_NotSupported,
2130 : "The NetCDF driver does not yet to support import of CF-1 Polar stereographic "
2131 0 : "without a std_parallel attribute.\n" );
2132 : }
2133 : dfScale =
2134 : poDS->FetchCopyParm( szGridMappingValue,
2135 : CF_PP_SCALE_FACTOR_ORIGIN,
2136 6 : 1.0 );
2137 :
2138 : dfCenterLon =
2139 : poDS->FetchCopyParm( szGridMappingValue,
2140 6 : CF_PP_VERT_LONG_FROM_POLE, 0.0 );
2141 :
2142 : dfFalseEasting =
2143 : poDS->FetchCopyParm( szGridMappingValue,
2144 6 : CF_PP_FALSE_EASTING, 0.0 );
2145 :
2146 : dfFalseNorthing =
2147 : poDS->FetchCopyParm( szGridMappingValue,
2148 6 : CF_PP_FALSE_NORTHING, 0.0 );
2149 :
2150 6 : bGotCfSRS = TRUE;
2151 : oSRS.SetPS( dfCenterLat, dfCenterLon, dfScale,
2152 6 : dfFalseEasting, dfFalseNorthing );
2153 :
2154 6 : if( !bGotGeogCS )
2155 0 : oSRS.SetWellKnownGeogCS( "WGS84" );
2156 : }
2157 :
2158 : /* -------------------------------------------------------------------- */
2159 : /* Stereographic */
2160 : /* -------------------------------------------------------------------- */
2161 :
2162 32 : else if ( EQUAL ( pszValue, CF_PT_STEREO ) ) {
2163 :
2164 : dfCenterLon =
2165 : poDS->FetchCopyParm( szGridMappingValue,
2166 6 : CF_PP_LON_PROJ_ORIGIN, 0.0 );
2167 :
2168 : dfCenterLat =
2169 : poDS->FetchCopyParm( szGridMappingValue,
2170 6 : CF_PP_LAT_PROJ_ORIGIN, 0.0 );
2171 :
2172 : dfScale =
2173 : poDS->FetchCopyParm( szGridMappingValue,
2174 : CF_PP_SCALE_FACTOR_ORIGIN,
2175 6 : 1.0 );
2176 :
2177 : dfFalseEasting =
2178 : poDS->FetchCopyParm( szGridMappingValue,
2179 6 : CF_PP_FALSE_EASTING, 0.0 );
2180 :
2181 : dfFalseNorthing =
2182 : poDS->FetchCopyParm( szGridMappingValue,
2183 6 : CF_PP_FALSE_NORTHING, 0.0 );
2184 :
2185 6 : bGotCfSRS = TRUE;
2186 : oSRS.SetStereographic( dfCenterLat, dfCenterLon, dfScale,
2187 6 : dfFalseEasting, dfFalseNorthing );
2188 :
2189 6 : if( !bGotGeogCS )
2190 0 : oSRS.SetWellKnownGeogCS( "WGS84" );
2191 : }
2192 :
2193 : /* -------------------------------------------------------------------- */
2194 : /* Is this Latitude/Longitude Grid, default */
2195 : /* -------------------------------------------------------------------- */
2196 :
2197 0 : } else if( EQUAL( szDimNameX,"lon" ) ) {
2198 0 : oSRS.SetWellKnownGeogCS( "WGS84" );
2199 :
2200 : } else {
2201 : // This would be too indiscrimant. But we should set
2202 : // it if we know the data is geographic.
2203 : //oSRS.SetWellKnownGeogCS( "WGS84" );
2204 : }
2205 : }
2206 : /* -------------------------------------------------------------------- */
2207 : /* Read projection coordinates */
2208 : /* -------------------------------------------------------------------- */
2209 :
2210 166 : nc_inq_varid( cdfid, poDS->papszDimName[nXDimID], &nVarDimXID );
2211 166 : nc_inq_varid( cdfid, poDS->papszDimName[nYDimID], &nVarDimYID );
2212 :
2213 166 : if( ( nVarDimXID != -1 ) && ( nVarDimYID != -1 ) ) {
2214 140 : pdfXCoord = (double *) CPLCalloc( xdim, sizeof(double) );
2215 140 : pdfYCoord = (double *) CPLCalloc( ydim, sizeof(double) );
2216 :
2217 140 : start[0] = 0;
2218 140 : edge[0] = xdim;
2219 : status = nc_get_vara_double( cdfid, nVarDimXID,
2220 140 : start, edge, pdfXCoord);
2221 :
2222 140 : edge[0] = ydim;
2223 : status = nc_get_vara_double( cdfid, nVarDimYID,
2224 140 : start, edge, pdfYCoord);
2225 :
2226 : /* -------------------------------------------------------------------- */
2227 : /* Check for bottom-up from the Y-axis order */
2228 : /* see bugs #4284 and #4251 */
2229 : /* -------------------------------------------------------------------- */
2230 :
2231 140 : if ( pdfYCoord[0] > pdfYCoord[1] )
2232 124 : poDS->bBottomUp = FALSE;
2233 : else
2234 16 : poDS->bBottomUp = TRUE;
2235 :
2236 140 : CPLDebug( "GDAL_netCDF", "set bBottomUp = %d from Y axis", poDS->bBottomUp );
2237 :
2238 :
2239 : /* -------------------------------------------------------------------- */
2240 : /* Is pixel spacing is uniform accross the map? */
2241 : /* -------------------------------------------------------------------- */
2242 :
2243 : /* -------------------------------------------------------------------- */
2244 : /* Check Longitude */
2245 : /* -------------------------------------------------------------------- */
2246 :
2247 140 : nSpacingBegin = (int) poDS->rint((pdfXCoord[1]-pdfXCoord[0]) * 1000);
2248 :
2249 140 : nSpacingMiddle = (int) poDS->rint((pdfXCoord[xdim / 2] -
2250 140 : pdfXCoord[(xdim / 2) + 1]) * 1000);
2251 :
2252 140 : nSpacingLast = (int) poDS->rint((pdfXCoord[xdim - 2] -
2253 140 : pdfXCoord[xdim-1]) * 1000);
2254 :
2255 279 : if( ( abs( nSpacingBegin ) == abs( nSpacingLast ) ) &&
2256 : ( abs( nSpacingBegin ) == abs( nSpacingMiddle ) ) &&
2257 : ( abs( nSpacingMiddle ) == abs( nSpacingLast ) ) ) {
2258 :
2259 : /* -------------------------------------------------------------------- */
2260 : /* Longitude is equally spaced, check latitude */
2261 : /* -------------------------------------------------------------------- */
2262 139 : nSpacingBegin = (int) poDS->rint((pdfYCoord[1]-pdfYCoord[0]) *
2263 139 : 1000);
2264 :
2265 139 : nSpacingMiddle = (int) poDS->rint((pdfYCoord[ydim / 2] -
2266 139 : pdfYCoord[(ydim / 2) + 1]) *
2267 139 : 1000);
2268 :
2269 139 : nSpacingLast = (int) poDS->rint((pdfYCoord[ydim - 2] -
2270 139 : pdfYCoord[ydim-1]) *
2271 139 : 1000);
2272 :
2273 :
2274 : /* -------------------------------------------------------------------- */
2275 : /* For Latitude we allow an error of 0.1 degrees for gaussian */
2276 : /* gridding */
2277 : /* -------------------------------------------------------------------- */
2278 :
2279 278 : if((( abs( abs(nSpacingBegin) - abs(nSpacingLast) ) ) < 100 ) &&
2280 : (( abs( abs(nSpacingBegin) - abs(nSpacingMiddle) ) ) < 100 ) &&
2281 : (( abs( abs(nSpacingMiddle) - abs(nSpacingLast) ) ) < 100) ) {
2282 :
2283 139 : if( ( abs( nSpacingBegin ) != abs( nSpacingLast ) ) ||
2284 : ( abs( nSpacingBegin ) != abs( nSpacingMiddle ) ) ||
2285 : ( abs( nSpacingMiddle ) != abs( nSpacingLast ) ) ) {
2286 :
2287 0 : CPLError(CE_Warning, 1,"Latitude grid not spaced evenly.\nSeting projection for grid spacing is within 0.1 degrees threshold.\n");
2288 :
2289 : }
2290 : /* -------------------------------------------------------------------- */
2291 : /* We have gridded data s we can set the Gereferencing info. */
2292 : /* -------------------------------------------------------------------- */
2293 :
2294 : /* -------------------------------------------------------------------- */
2295 : /* Enable GeoTransform */
2296 : /* -------------------------------------------------------------------- */
2297 : /* ----------------------------------------------------------*/
2298 : /* In the following "actual_range" and "node_offset" */
2299 : /* are attributes used by netCDF files created by GMT. */
2300 : /* If we find them we know how to proceed. Else, use */
2301 : /* the original algorithm. */
2302 : /* --------------------------------------------------------- */
2303 : double dummy[2], xMinMax[2], yMinMax[2];
2304 139 : int node_offset = 0;
2305 :
2306 139 : bGotGeoTransform = TRUE;
2307 :
2308 139 : nc_get_att_int (cdfid, NC_GLOBAL, "node_offset", &node_offset);
2309 :
2310 139 : if (!nc_get_att_double (cdfid, nVarDimXID, "actual_range", dummy)) {
2311 5 : xMinMax[0] = dummy[0];
2312 5 : xMinMax[1] = dummy[1];
2313 : }
2314 : else {
2315 134 : xMinMax[0] = pdfXCoord[0];
2316 134 : xMinMax[1] = pdfXCoord[xdim-1];
2317 134 : node_offset = 0;
2318 : }
2319 :
2320 139 : if (!nc_get_att_double (cdfid, nVarDimYID, "actual_range", dummy)) {
2321 5 : yMinMax[0] = dummy[0];
2322 5 : yMinMax[1] = dummy[1];
2323 : }
2324 : else {
2325 134 : yMinMax[0] = pdfYCoord[0];
2326 134 : yMinMax[1] = pdfYCoord[ydim-1];
2327 134 : node_offset = 0;
2328 : }
2329 :
2330 : /* for CF-1 conventions, assume bottom first */
2331 : /* was here, now is higher up in the code */
2332 : // if( ( EQUAL( szDimNameY, "lat" ) || EQUAL( szDimNameY, "y" ) )
2333 : // && pdfYCoord[0] < pdfYCoord[1] )
2334 : // poDS->bBottomUp = TRUE;
2335 :
2336 : /* Check for reverse order of y-coordinate */
2337 139 : if ( yMinMax[0] > yMinMax[1] ) {
2338 123 : dummy[0] = yMinMax[1];
2339 123 : dummy[1] = yMinMax[0];
2340 123 : yMinMax[0] = dummy[0];
2341 123 : yMinMax[1] = dummy[1];
2342 : }
2343 :
2344 : /* ----------------------------------------------------------*/
2345 : /* Many netcdf files are weather files distributed */
2346 : /* in km for the x/y resolution. This isn't perfect, */
2347 : /* but geotransforms can be terribly off if this isn't */
2348 : /* checked and accounted for. Maybe one more level of */
2349 : /* checking (grid_mapping_value#GRIB_param_Dx, or */
2350 : /* x#grid_spacing), but those are not cf tags. */
2351 : /* Have to change metadata value if change Create() to */
2352 : /* write cf tags */
2353 : /* ----------------------------------------------------------*/
2354 :
2355 : //check units for x and y, expand to other values
2356 : //and conversions.
2357 139 : if( oSRS.IsProjected( ) ) {
2358 93 : strcpy( szTemp, "x" );
2359 93 : strcat( szTemp, "#units" );
2360 : pszValue = CSLFetchNameValue( poDS->papszMetadata,
2361 93 : szTemp );
2362 93 : if( pszValue != NULL ) {
2363 93 : pszUnits = pszValue;
2364 93 : if( EQUAL( pszValue, "km" ) ) {
2365 5 : xMinMax[0] = xMinMax[0] * 1000;
2366 5 : xMinMax[1] = xMinMax[1] * 1000;
2367 : }
2368 : }
2369 93 : strcpy( szTemp, "y" );
2370 93 : strcat( szTemp, "#units" );
2371 : pszValue = CSLFetchNameValue( poDS->papszMetadata,
2372 93 : szTemp );
2373 93 : if( pszValue != NULL ) {
2374 : /* TODO: see how to deal with diff. values */
2375 : // if ( ! EQUAL( pszValue, szUnits ) )
2376 : // strcpy( szUnits, "\0" );
2377 93 : if( EQUAL( pszValue, "km" ) ) {
2378 5 : yMinMax[0] = yMinMax[0] * 1000;
2379 5 : yMinMax[1] = yMinMax[1] * 1000;
2380 : }
2381 : }
2382 : }
2383 :
2384 139 : adfTempGeoTransform[0] = xMinMax[0];
2385 139 : adfTempGeoTransform[2] = 0;
2386 139 : adfTempGeoTransform[3] = yMinMax[1];
2387 139 : adfTempGeoTransform[4] = 0;
2388 139 : adfTempGeoTransform[1] = ( xMinMax[1] - xMinMax[0] ) /
2389 139 : ( poDS->nRasterXSize + (node_offset - 1) );
2390 139 : adfTempGeoTransform[5] = ( yMinMax[0] - yMinMax[1] ) /
2391 139 : ( poDS->nRasterYSize + (node_offset - 1) );
2392 :
2393 : /* -------------------------------------------------------------------- */
2394 : /* Compute the center of the pixel */
2395 : /* -------------------------------------------------------------------- */
2396 139 : if ( !node_offset ) { // Otherwise its already the pixel center
2397 139 : adfTempGeoTransform[0] -= (adfTempGeoTransform[1] / 2);
2398 139 : adfTempGeoTransform[3] -= (adfTempGeoTransform[5] / 2);
2399 : }
2400 : }// end if (Latitude is equally spaced, within 0.1 degrees)
2401 : else {
2402 : CPLDebug( "GDAL_netCDF",
2403 0 : "Latitude is not equally spaced." );
2404 : }
2405 : }// end if (Longitude is equally spaced)
2406 : else {
2407 : CPLDebug( "GDAL_netCDF",
2408 1 : "Longitude is not equally spaced." );
2409 : }
2410 :
2411 140 : CPLFree( pdfXCoord );
2412 140 : CPLFree( pdfYCoord );
2413 : }// end if (has dims)
2414 :
2415 : CPLDebug( "GDAL_netCDF",
2416 : "bGotGeogCS=%d bGotCfSRS=%d bGotGeoTransform=%d",
2417 166 : bGotGeogCS, bGotCfSRS, bGotGeoTransform );
2418 :
2419 : /* -------------------------------------------------------------------- */
2420 : /* Set Projection if we got a geotransform */
2421 : /* -------------------------------------------------------------------- */
2422 166 : if ( bGotGeoTransform ) {
2423 : /* Set SRS Units */
2424 : /* TODO: check for other units */
2425 139 : if ( pszUnits != NULL && ! EQUAL(pszUnits,"") ) {
2426 93 : if ( EQUAL(pszUnits,"m") ) {
2427 88 : oSRS.SetLinearUnits( CF_UNITS_M, 1.0 );
2428 88 : oSRS.SetAuthority( "PROJCS|UNIT", "EPSG", 9001 );
2429 : }
2430 5 : else if ( EQUAL(pszUnits,"km") ) {
2431 5 : oSRS.SetLinearUnits( CF_UNITS_M, 1000.0 );
2432 5 : oSRS.SetAuthority( "PROJCS|UNIT", "EPSG", 9001 );
2433 : }
2434 0 : else if ( EQUALN(pszUnits,"degrees",7) ) {
2435 0 : oSRS.SetAngularUnits( CF_UNITS_D, CPLAtof(SRS_UA_DEGREE_CONV) );
2436 0 : oSRS.SetAuthority( "GEOGCS|UNIT", "EPSG", 9108 );
2437 : }
2438 : // else
2439 : // oSRS.SetLinearUnits(pszUnits, 1.0);
2440 : }
2441 :
2442 : /* Set GeoTransform */
2443 139 : SetGeoTransform( adfTempGeoTransform );
2444 : /* Set Projection */
2445 139 : oSRS.exportToWkt( &(pszTempProjection) );
2446 139 : SetProjection( pszTempProjection );
2447 139 : CPLFree( pszTempProjection );
2448 139 : CPLDebug( "GDAL_netCDF", "set WKT from CF [%s]\n", pszProjection );
2449 : }
2450 27 : else if ( bGotGeogCS || bGotCfSRS ) {
2451 0 : CPLError(CE_Warning, 1,"got SRS but no geotransform from CF!");
2452 : }
2453 : /* -------------------------------------------------------------------- */
2454 : /* Process custom GDAL values (spatial_ref, GeoTransform) */
2455 : /* -------------------------------------------------------------------- */
2456 166 : if( !EQUAL( szGridMappingValue, "" ) ) {
2457 :
2458 116 : if( pszWKT != NULL ) {
2459 :
2460 : /* -------------------------------------------------------------------- */
2461 : /* Compare CRS obtained from CF attributes and GDAL WKT */
2462 : /* If possible use the more complete GDAL WKT */
2463 : /* -------------------------------------------------------------------- */
2464 : /* Set the CRS to the one written by GDAL */
2465 111 : if ( ! bGotCfSRS || poDS->pszProjection == NULL || ! bIsGdalCfFile ) {
2466 0 : bGotGdalSRS = TRUE;
2467 0 : SetProjection( pszWKT );
2468 0 : CPLDebug( "GDAL_netCDF", "set WKT from GDAL [%s]\n", pszProjection );
2469 : }
2470 : else { /* use the SRS from GDAL if it doesn't conflict with the one from CF */
2471 111 : char *pszProjectionGDAL = (char*) pszWKT ;
2472 111 : OGRSpatialReference oSRSGDAL;
2473 111 : oSRSGDAL.importFromWkt( &pszProjectionGDAL );
2474 : /* set datum to unknown or else datums will not match, see bug #4281 */
2475 111 : if ( oSRSGDAL.GetAttrNode( "DATUM" ) )
2476 111 : oSRSGDAL.GetAttrNode( "DATUM" )->GetChild(0)->SetValue( "unknown" );
2477 : /* need this for setprojection autotest */
2478 111 : if ( oSRSGDAL.GetAttrNode( "PROJCS" ) )
2479 88 : oSRSGDAL.GetAttrNode( "PROJCS" )->GetChild(0)->SetValue( "unnamed" );
2480 111 : if ( oSRSGDAL.GetAttrNode( "GEOGCS" ) )
2481 111 : oSRSGDAL.GetAttrNode( "GEOGCS" )->GetChild(0)->SetValue( "unknown" );
2482 111 : oSRSGDAL.GetRoot()->StripNodes( "UNIT" );
2483 111 : if ( oSRS.IsSame(&oSRSGDAL) ) {
2484 : // printf("ARE SAME, using GDAL WKT\n");
2485 111 : bGotGdalSRS = TRUE;
2486 111 : SetProjection( pszWKT );
2487 111 : CPLDebug( "GDAL_netCDF", "set WKT from GDAL [%s]\n", pszProjection );
2488 : }
2489 : else {
2490 : CPLDebug( "GDAL_netCDF",
2491 : "got WKT from GDAL \n[%s]\nbut not using it because conflicts with CF\n[%s]\n",
2492 0 : pszWKT, poDS->pszProjection );
2493 111 : }
2494 : }
2495 : /* -------------------------------------------------------------------- */
2496 : /* Look for GeoTransform Array, if not found previously */
2497 : /* -------------------------------------------------------------------- */
2498 111 : if ( !bGotGeoTransform ) {
2499 :
2500 : /* TODO read the GT values and detect for conflict with CF */
2501 : /* this could resolve the GT precision loss issue */
2502 :
2503 0 : if( pszGeoTransform != NULL ) {
2504 : papszGeoTransform = CSLTokenizeString2( pszGeoTransform,
2505 : " ",
2506 0 : CSLT_HONOURSTRINGS );
2507 0 : bGotGeoTransform = TRUE;
2508 :
2509 0 : adfTempGeoTransform[0] = atof( papszGeoTransform[0] );
2510 0 : adfTempGeoTransform[1] = atof( papszGeoTransform[1] );
2511 0 : adfTempGeoTransform[2] = atof( papszGeoTransform[2] );
2512 0 : adfTempGeoTransform[3] = atof( papszGeoTransform[3] );
2513 0 : adfTempGeoTransform[4] = atof( papszGeoTransform[4] );
2514 0 : adfTempGeoTransform[5] = atof( papszGeoTransform[5] );
2515 : /* -------------------------------------------------------------------- */
2516 : /* Look for corner array values */
2517 : /* -------------------------------------------------------------------- */
2518 : } else {
2519 0 : double dfNN=0.0, dfSN=0.0, dfEE=0.0, dfWE=0.0;
2520 0 : int bGotNN=FALSE, bGotSN=FALSE, bGotEE=FALSE, bGotWE=FALSE;
2521 : // CPLDebug( "GDAL_netCDF", "looking for geotransform corners\n" );
2522 0 : strcpy(szTemp,szGridMappingValue);
2523 0 : strcat( szTemp, "#" );
2524 0 : strcat( szTemp, "Northernmost_Northing");
2525 0 : pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp);
2526 :
2527 0 : if( pszValue != NULL ) {
2528 0 : dfNN = atof( pszValue );
2529 0 : bGotNN = TRUE;
2530 : }
2531 0 : strcpy(szTemp,szGridMappingValue);
2532 0 : strcat( szTemp, "#" );
2533 0 : strcat( szTemp, "Southernmost_Northing");
2534 0 : pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp);
2535 :
2536 0 : if( pszValue != NULL ) {
2537 0 : dfSN = atof( pszValue );
2538 0 : bGotSN = TRUE;
2539 : }
2540 :
2541 0 : strcpy(szTemp,szGridMappingValue);
2542 0 : strcat( szTemp, "#" );
2543 0 : strcat( szTemp, "Easternmost_Easting");
2544 0 : pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp);
2545 :
2546 0 : if( pszValue != NULL ) {
2547 0 : dfEE = atof( pszValue );
2548 0 : bGotEE = TRUE;
2549 : }
2550 :
2551 0 : strcpy(szTemp,szGridMappingValue);
2552 0 : strcat( szTemp, "#" );
2553 0 : strcat( szTemp, "Westernmost_Easting");
2554 0 : pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp);
2555 :
2556 0 : if( pszValue != NULL ) {
2557 0 : dfWE = atof( pszValue );
2558 0 : bGotWE = TRUE;
2559 : }
2560 :
2561 : /* Only set the GeoTransform if we got all the values */
2562 0 : if ( bGotNN && bGotSN && bGotEE && bGotWE ) {
2563 0 : adfTempGeoTransform[0] = dfWE;
2564 : adfTempGeoTransform[1] = (dfEE - dfWE) /
2565 0 : ( poDS->GetRasterXSize() - 1 );
2566 0 : adfTempGeoTransform[2] = 0.0;
2567 0 : adfTempGeoTransform[3] = dfNN;
2568 0 : adfTempGeoTransform[4] = 0.0;
2569 : adfTempGeoTransform[5] = (dfSN - dfNN) /
2570 0 : ( poDS->GetRasterYSize() - 1 );
2571 : /* -------------------------------------------------------------------- */
2572 : /* Compute the center of the pixel */
2573 : /* -------------------------------------------------------------------- */
2574 : adfTempGeoTransform[0] = dfWE
2575 0 : - (adfTempGeoTransform[1] / 2);
2576 :
2577 : adfTempGeoTransform[3] = dfNN
2578 0 : - (adfTempGeoTransform[5] / 2);
2579 :
2580 :
2581 0 : bGotGeoTransform = TRUE;
2582 : }
2583 : } // (pszGeoTransform != NULL)
2584 0 : CSLDestroy( papszGeoTransform );
2585 :
2586 : /* Issue a warning if we did not get a geotransform from GDAL */
2587 0 : if ( !bGotGeoTransform ) {
2588 0 : CPLError(CE_Warning, 1,"got SRS but not geotransform from GDAL!");
2589 : }
2590 : } // (!bGotGeoTransform)
2591 : }
2592 : }
2593 :
2594 166 : if ( bGotGeoTransform ) {
2595 : CPLDebug( "GDAL_netCDF", "Got GeoTransform:"
2596 : " %.16g, %.16g, %.16g"
2597 : " %.16g, %.16g, %.16g",
2598 : adfGeoTransform[0], adfGeoTransform[1],
2599 : adfGeoTransform[2], adfGeoTransform[3],
2600 139 : adfGeoTransform[4], adfGeoTransform[5] );
2601 : }
2602 :
2603 : /* -------------------------------------------------------------------- */
2604 : /* Search for Well-known GeogCS if got only CF WKT */
2605 : /* Disabled for now, as a named datum also include control points */
2606 : /* (see mailing list and bug#4281 */
2607 : /* For example, WGS84 vs. GDA94 (EPSG:3577) - AEA in netcdf_cf.py */
2608 : /* -------------------------------------------------------------------- */
2609 : /* disabled for now, but could be set in a config option */
2610 166 : bLookForWellKnownGCS = FALSE;
2611 166 : if ( bLookForWellKnownGCS && bGotCfSRS && ! bGotGdalSRS ) {
2612 : /* ET - could use a more exhaustive method by scanning all EPSG codes in data/gcs.csv */
2613 : /* as proposed by Even in the gdal-dev mailing list "help for comparing two WKT" */
2614 : /* this code could be contributed to a new function */
2615 : /* OGRSpatialReference * OGRSpatialReference::FindMatchingGeogCS( const OGRSpatialReference *poOther ) */
2616 0 : CPLDebug( "GDAL_netCDF", "Searching for Well-known GeogCS" );
2617 0 : const char *pszWKGCSList[] = { "WGS84", "WGS72", "NAD27", "NAD83" };
2618 0 : char *pszWKGCS = NULL;
2619 0 : oSRS.exportToPrettyWkt( &pszWKGCS );
2620 0 : for( size_t i=0; i<sizeof(pszWKGCSList)/8; i++ ) {
2621 0 : pszWKGCS = CPLStrdup( pszWKGCSList[i] );
2622 0 : OGRSpatialReference oSRSTmp;
2623 0 : oSRSTmp.SetWellKnownGeogCS( pszWKGCSList[i] );
2624 : /* set datum to unknown, bug #4281 */
2625 0 : if ( oSRSTmp.GetAttrNode( "DATUM" ) )
2626 0 : oSRSTmp.GetAttrNode( "DATUM" )->GetChild(0)->SetValue( "unknown" );
2627 : /* could use OGRSpatialReference::StripCTParms() but let's keep TOWGS84 */
2628 0 : oSRSTmp.GetRoot()->StripNodes( "AXIS" );
2629 0 : oSRSTmp.GetRoot()->StripNodes( "AUTHORITY" );
2630 0 : oSRSTmp.GetRoot()->StripNodes( "EXTENSION" );
2631 :
2632 0 : oSRSTmp.exportToPrettyWkt( &pszWKGCS );
2633 0 : if ( oSRS.IsSameGeogCS(&oSRSTmp) ) {
2634 0 : oSRS.SetWellKnownGeogCS( pszWKGCSList[i] );
2635 0 : oSRS.exportToWkt( &(pszTempProjection) );
2636 0 : SetProjection( pszTempProjection );
2637 0 : CPLFree( pszTempProjection );
2638 : }
2639 : }
2640 166 : }
2641 :
2642 :
2643 166 : }
2644 :
2645 : /************************************************************************/
2646 : /* SetProjection() */
2647 : /************************************************************************/
2648 334 : CPLErr netCDFDataset::SetProjection( const char * pszNewProjection )
2649 : {
2650 : /* TODO look if proj. already defined, like in geotiff */
2651 334 : if( pszNewProjection == NULL )
2652 : {
2653 0 : CPLError( CE_Failure, CPLE_AppDefined, "NULL projection." );
2654 0 : return CE_Failure;
2655 : }
2656 :
2657 334 : if( bSetProjection && (GetAccess() == GA_Update) )
2658 : {
2659 : CPLError( CE_Warning, CPLE_AppDefined,
2660 : "netCDFDataset::SetProjection() should only be called once "
2661 : "in update mode!\npszNewProjection=\n%s",
2662 0 : pszNewProjection );
2663 : }
2664 :
2665 334 : CPLDebug( "GDAL_netCDF", "SetProjection, WKT = %s", pszNewProjection );
2666 :
2667 334 : if( !EQUALN(pszNewProjection,"GEOGCS",6)
2668 : && !EQUALN(pszNewProjection,"PROJCS",6)
2669 : && !EQUAL(pszNewProjection,"") )
2670 : {
2671 : CPLError( CE_Failure, CPLE_AppDefined,
2672 : "Only OGC WKT GEOGCS and PROJCS Projections supported for writing to NetCDF.\n"
2673 : "%s not supported.",
2674 0 : pszNewProjection );
2675 :
2676 0 : return CE_Failure;
2677 : }
2678 :
2679 334 : CPLFree( pszProjection );
2680 334 : pszProjection = CPLStrdup( pszNewProjection );
2681 :
2682 334 : if( GetAccess() == GA_Update )
2683 : {
2684 84 : if ( bSetGeoTransform && ! bSetProjection ) {
2685 27 : bSetProjection = TRUE;
2686 27 : return AddProjectionVars();
2687 : }
2688 : }
2689 :
2690 307 : bSetProjection = TRUE;
2691 :
2692 307 : return CE_None;
2693 :
2694 : }
2695 :
2696 : /************************************************************************/
2697 : /* SetGeoTransform() */
2698 : /************************************************************************/
2699 :
2700 224 : CPLErr netCDFDataset::SetGeoTransform ( double * padfTransform )
2701 : {
2702 224 : memcpy( adfGeoTransform, padfTransform, sizeof(double)*6 );
2703 : // bGeoTransformValid = TRUE;
2704 : // bGeoTIFFInfoChanged = TRUE;
2705 :
2706 : CPLDebug( "GDAL_netCDF",
2707 : "SetGeoTransform(%f,%f,%f,%f,%f,%f)",
2708 : padfTransform[0],padfTransform[1],padfTransform[2],
2709 224 : padfTransform[3],padfTransform[4],padfTransform[5]);
2710 :
2711 224 : if( GetAccess() == GA_Update )
2712 : {
2713 85 : if ( bSetProjection && ! bSetGeoTransform ) {
2714 0 : bSetGeoTransform = TRUE;
2715 0 : return AddProjectionVars();
2716 : }
2717 : }
2718 :
2719 224 : bSetGeoTransform = TRUE;
2720 :
2721 224 : return CE_None;
2722 :
2723 : }
2724 :
2725 : /************************************************************************/
2726 : /* AddProjectionVars() */
2727 : /************************************************************************/
2728 :
2729 85 : CPLErr netCDFDataset::AddProjectionVars( GDALProgressFunc pfnProgress,
2730 : void * pProgressData )
2731 : {
2732 85 : OGRSpatialReference oSRS;
2733 85 : int NCDFVarID = -1;
2734 85 : double dfTemp = 0.0;
2735 85 : const char *pszValue = NULL;
2736 : char szTemp[ NCDF_MAX_STR_LEN ];
2737 :
2738 : char szGeoTransform[ NCDF_MAX_STR_LEN ];
2739 85 : *szGeoTransform = '\0';
2740 85 : char *pszWKT = NULL;
2741 :
2742 85 : int bWriteGridMapping = FALSE;
2743 85 : int bWriteLonLat = FALSE;
2744 85 : int bWriteGDALTags = FALSE;
2745 85 : int bWriteGeoTransform = FALSE;
2746 :
2747 85 : nc_type eLonLatType = NC_NAT;
2748 : // int nLonSize = nRasterXSize;
2749 : // int nLatSize = nRasterYSize;
2750 85 : int nVarLonID=-1, nVarLatID=-1;
2751 85 : int nVarXID=-1, nVarYID=-1;
2752 :
2753 : char szNetcdfProjection[ NC_MAX_NAME ];
2754 85 : szNetcdfProjection[0]='\0';
2755 :
2756 85 : bAddedProjectionVars = TRUE;
2757 :
2758 85 : pszWKT = (char *) pszProjection;
2759 85 : oSRS.importFromWkt( &pszWKT );
2760 :
2761 85 : if( oSRS.IsProjected() )
2762 41 : bIsProjected = TRUE;
2763 44 : else if( oSRS.IsGeographic() )
2764 39 : bIsGeographic = TRUE;
2765 :
2766 : CPLDebug( "GDAL_netCDF", "SetProjection, WKT now = [%s]\nprojected: %d geographic: %d",
2767 85 : pszProjection,bIsProjected,bIsGeographic );
2768 :
2769 85 : if ( ! bSetGeoTransform )
2770 : CPLError( CE_Warning, CPLE_AppDefined,
2771 0 : "netCDFDataset::DefineDims() called, but GeoTransform has not yet been defined!" );
2772 :
2773 85 : if ( ! bSetProjection )
2774 : CPLError( CE_Warning, CPLE_AppDefined,
2775 1 : "netCDFDataset::DefineDims() called, but Projection has not yet been defined!" );
2776 :
2777 : /* go through ProcessOptions again as defaults depend on projection type */
2778 : // ProcessCreationOptions( );
2779 :
2780 : /* process options */
2781 85 : if( bIsProjected )
2782 : {
2783 41 : int bIsCfProjection = NCDFIsCfProjection( oSRS.GetAttrValue( "PROJECTION" ) );
2784 41 : bWriteGridMapping = TRUE;
2785 41 : bWriteGDALTags = CSLFetchBoolean( papszCreationOptions, "WRITE_GDAL_TAGS", TRUE );
2786 : /* force WRITE_GDAL_TAGS if is not a CF projection */
2787 41 : if ( ! bWriteGDALTags && ! bIsCfProjection )
2788 0 : bWriteGDALTags = TRUE;
2789 41 : if ( bWriteGDALTags )
2790 41 : bWriteGeoTransform = TRUE;
2791 :
2792 41 : pszValue = CSLFetchNameValueDef(papszCreationOptions,"WRITE_LONLAT", "NO");
2793 41 : if ( EQUAL( pszValue, "IF_NEEDED" ) ) {
2794 0 : if ( bIsCfProjection )
2795 0 : bWriteLonLat = FALSE;
2796 : else
2797 0 : bWriteLonLat = TRUE;
2798 : }
2799 41 : else bWriteLonLat = CSLTestBoolean( pszValue );
2800 : // if ( bWriteLonLat == TRUE ) {
2801 : // nLonSize = nRasterXSize * nRasterYSize;
2802 : // nLatSize = nRasterXSize * nRasterYSize;
2803 : // }
2804 41 : eLonLatType = NC_FLOAT;
2805 41 : pszValue = CSLFetchNameValueDef(papszCreationOptions,"TYPE_LONLAT", "FLOAT");
2806 41 : if ( EQUAL(pszValue, "DOUBLE" ) )
2807 0 : eLonLatType = NC_DOUBLE;
2808 : }
2809 : else
2810 : {
2811 : /* files without a Datum will not have a grid_mapping variable and geographic information */
2812 44 : if ( bIsGeographic ) bWriteGridMapping = TRUE;
2813 5 : else bWriteGridMapping = FALSE;
2814 44 : bWriteGDALTags = CSLFetchBoolean( papszCreationOptions, "WRITE_GDAL_TAGS", bWriteGridMapping );
2815 44 : if ( bWriteGDALTags )
2816 39 : bWriteGeoTransform = TRUE;
2817 :
2818 44 : pszValue = CSLFetchNameValueDef(papszCreationOptions,"WRITE_LONLAT", "YES");
2819 44 : if ( EQUAL( pszValue, "IF_NEEDED" ) )
2820 0 : bWriteLonLat = TRUE;
2821 44 : else bWriteLonLat = CSLTestBoolean( pszValue );
2822 : /* Don't write lon/lat if no source geotransform */
2823 44 : if ( ! bSetGeoTransform )
2824 0 : bWriteLonLat = FALSE;
2825 : /* If we don't write lon/lat, set dimnames to X/Y and write gdal tags*/
2826 44 : if ( ! bWriteLonLat ) {
2827 : CPLError( CE_Warning, CPLE_AppDefined,
2828 0 : "creating geographic file without lon/lat values!");
2829 0 : if ( bSetGeoTransform ) {
2830 0 : bWriteGDALTags = TRUE; //not desireable if no geotransform
2831 0 : bWriteGeoTransform = TRUE;
2832 : }
2833 : // pszLonDimName = NCDF_DIMNAME_X;
2834 : // pszLatDimName = NCDF_DIMNAME_Y;
2835 : // bBottomUp = FALSE;
2836 : }
2837 :
2838 : // nLonSize = nRasterXSize;
2839 : // nLatSize = nRasterYSize;
2840 :
2841 44 : eLonLatType = NC_DOUBLE;
2842 44 : pszValue = CSLFetchNameValueDef(papszCreationOptions,"TYPE_LONLAT", "DOUBLE");
2843 44 : if ( EQUAL(pszValue, "FLOAT" ) )
2844 0 : eLonLatType = NC_FLOAT;
2845 : }
2846 :
2847 : /* make sure we write grid_mapping if we need to write GDAL tags */
2848 85 : if ( bWriteGDALTags ) bWriteGridMapping = TRUE;
2849 :
2850 : CPLDebug( "GDAL_netCDF",
2851 : "bIsProjected=%d bIsGeographic=%d bWriteGridMapping=%d bWriteGDALTags=%d bWriteLonLat=%d bBottomUp=%d",
2852 85 : bIsProjected,bIsGeographic,bWriteGridMapping,bWriteGDALTags,bWriteLonLat,bBottomUp );
2853 :
2854 :
2855 : /* -------------------------------------------------------------------- */
2856 : /* Define dimension names */
2857 : /* -------------------------------------------------------------------- */
2858 : /* make sure we are in define mode */
2859 85 : SetDefineMode( TRUE );
2860 :
2861 :
2862 : /* -------------------------------------------------------------------- */
2863 : /* Rename dimensions if lon/lat */
2864 : /* -------------------------------------------------------------------- */
2865 85 : if( ! bIsProjected )
2866 : {
2867 : /* rename dims to lat/lon */
2868 44 : papszDimName.Clear(); //if we add other dims one day, this has to change
2869 44 : papszDimName.AddString( NCDF_DIMNAME_LAT );
2870 44 : papszDimName.AddString( NCDF_DIMNAME_LON );
2871 :
2872 44 : status = nc_rename_dim(cdfid, nYDimID, NCDF_DIMNAME_LAT );
2873 44 : NCDF_ERR(status);
2874 44 : status = nc_rename_dim(cdfid, nXDimID, NCDF_DIMNAME_LON );
2875 44 : NCDF_ERR(status);
2876 : // anBandDims[0] = nLatDimID;
2877 : // anBandDims[1] = nLonDimID;
2878 : }
2879 :
2880 : /* -------------------------------------------------------------------- */
2881 : /* Write projection attributes */
2882 : /* -------------------------------------------------------------------- */
2883 85 : if( bIsProjected )
2884 : {
2885 : /* -------------------------------------------------------------------- */
2886 : /* Write CF-1.5 compliant Projected attributes */
2887 : /* -------------------------------------------------------------------- */
2888 :
2889 41 : const OGR_SRSNode *poPROJCS = oSRS.GetAttrNode( "PROJCS" );
2890 : const char *pszProjName;
2891 41 : pszProjName = oSRS.GetAttrValue( "PROJECTION" );
2892 :
2893 : /* Basic Projection info (grid_mapping and datum) */
2894 823 : for( int i=0; poNetcdfSRS_PT[i].WKT_SRS != NULL; i++ ) {
2895 823 : if( EQUAL( poNetcdfSRS_PT[i].WKT_SRS, pszProjName ) ) {
2896 : CPLDebug( "GDAL_netCDF", "GDAL PROJECTION = %s , NCDF PROJECTION = %s",
2897 : poNetcdfSRS_PT[i].WKT_SRS,
2898 41 : poNetcdfSRS_PT[i].CF_SRS);
2899 41 : strcpy( szNetcdfProjection, poNetcdfSRS_PT[i].CF_SRS );
2900 : CPLDebug( "GDAL_netCDF", "nc_def_var(%d,%s,%d)",
2901 41 : cdfid, poNetcdfSRS_PT[i].CF_SRS, NC_CHAR );
2902 : status = nc_def_var( cdfid,
2903 : poNetcdfSRS_PT[i].CF_SRS,
2904 : NC_CHAR,
2905 41 : 0, NULL, &NCDFVarID );
2906 41 : NCDF_ERR(status);
2907 41 : break;
2908 : }
2909 : }
2910 : nc_put_att_text( cdfid, NCDFVarID, CF_GRD_MAPPING_NAME,
2911 : strlen( szNetcdfProjection ),
2912 41 : szNetcdfProjection );
2913 :
2914 : /* Various projection attributes */
2915 : // PDS: keep in synch with SetProjection function
2916 41 : NCDFWriteProjAttribs(poPROJCS, pszProjName, cdfid, NCDFVarID);
2917 :
2918 : }
2919 : else
2920 : {
2921 : /* -------------------------------------------------------------------- */
2922 : /* Write CF-1.5 compliant Geographics attributes */
2923 : /* Note: WKT information will not be preserved (e.g. WGS84) */
2924 : /* -------------------------------------------------------------------- */
2925 :
2926 44 : if( bWriteGridMapping == TRUE )
2927 : {
2928 39 : strcpy( szNetcdfProjection, "crs" );
2929 : CPLDebug( "GDAL_netCDF", "nc_def_var(%d,%s,%d)",
2930 39 : cdfid, szNetcdfProjection, NC_CHAR );
2931 : nc_def_var( cdfid, szNetcdfProjection, NC_CHAR,
2932 39 : 0, NULL, &NCDFVarID );
2933 39 : NCDF_ERR(status);
2934 : nc_put_att_text( cdfid, NCDFVarID, CF_GRD_MAPPING_NAME,
2935 : strlen(CF_PT_LATITUDE_LONGITUDE),
2936 39 : CF_PT_LATITUDE_LONGITUDE );
2937 : }
2938 :
2939 : }
2940 :
2941 : /* -------------------------------------------------------------------- */
2942 : /* Write CF-1.5 compliant common attributes */
2943 : /* -------------------------------------------------------------------- */
2944 :
2945 : /* DATUM information */
2946 85 : dfTemp = oSRS.GetPrimeMeridian();
2947 : nc_put_att_double( cdfid, NCDFVarID, CF_PP_LONG_PRIME_MERIDIAN,
2948 85 : NC_DOUBLE, 1, &dfTemp );
2949 85 : dfTemp = oSRS.GetSemiMajor();
2950 : nc_put_att_double( cdfid, NCDFVarID, CF_PP_SEMI_MAJOR_AXIS,
2951 85 : NC_DOUBLE, 1, &dfTemp );
2952 85 : dfTemp = oSRS.GetInvFlattening();
2953 : nc_put_att_double( cdfid, NCDFVarID, CF_PP_INVERSE_FLATTENING,
2954 85 : NC_DOUBLE, 1, &dfTemp );
2955 :
2956 :
2957 : /* Optional GDAL custom projection tags */
2958 85 : if ( bWriteGDALTags == TRUE ) {
2959 :
2960 80 : *szGeoTransform = '\0';
2961 560 : for( int i=0; i<6; i++ ) {
2962 : sprintf( szTemp, "%.16g ",
2963 480 : adfGeoTransform[i] );
2964 480 : strcat( szGeoTransform, szTemp );
2965 : }
2966 80 : CPLDebug( "GDAL_netCDF", "szGeoTranform = %s", szGeoTransform );
2967 :
2968 : // if ( strlen(pszProj4Defn) > 0 ) {
2969 : // nc_put_att_text( cdfid, NCDFVarID, "proj4",
2970 : // strlen( pszProj4Defn ), pszProj4Defn );
2971 : // }
2972 : nc_put_att_text( cdfid, NCDFVarID, NCDF_SPATIAL_REF,
2973 80 : strlen( pszProjection ), pszProjection );
2974 : /* for now write the geotransform for back-compat or else
2975 : the old (1.8.1) driver overrides the CF geotransform with
2976 : empty values from dfNN, dfSN, dfEE, dfWE; */
2977 : /* TODO: fix this in 1.8 branch */
2978 80 : if ( bWriteGeoTransform && bSetGeoTransform ) {
2979 : nc_put_att_text( cdfid, NCDFVarID, NCDF_GEOTRANSFORM,
2980 : strlen( szGeoTransform ),
2981 80 : szGeoTransform );
2982 : }
2983 : }
2984 :
2985 : /* -------------------------------------------------------------------- */
2986 : /* Write Projection var in Bands */
2987 : /* -------------------------------------------------------------------- */
2988 85 : if( bWriteGridMapping == TRUE ) {
2989 199 : for( int i=1; i <= nBands; i++ ) {
2990 : netCDFRasterBand *poSrcBand =
2991 119 : (netCDFRasterBand *) GetRasterBand( i );
2992 : status = nc_put_att_text( cdfid, poSrcBand->nZId,
2993 : CF_GRD_MAPPING,
2994 : strlen( szNetcdfProjection ),
2995 119 : szNetcdfProjection );
2996 119 : NCDF_ERR(status);
2997 119 : if ( bWriteLonLat == TRUE ) {
2998 : status = nc_put_att_text( cdfid, poSrcBand->nZId,
2999 : CF_COORDINATES,
3000 78 : strlen( NCDF_LONLAT ), NCDF_LONLAT );
3001 78 : NCDF_ERR(status);
3002 : }
3003 : }
3004 : }
3005 :
3006 85 : pfnProgress( 0.10, NULL, pProgressData );
3007 :
3008 : /* -------------------------------------------------------------------- */
3009 : /* Write CF Projection vars */
3010 : /* -------------------------------------------------------------------- */
3011 :
3012 : /* -------------------------------------------------------------------- */
3013 : /* Write X/Y attributes */
3014 : /* -------------------------------------------------------------------- */
3015 85 : if( bIsProjected )
3016 : {
3017 : /* X */
3018 : int anXDims[1];
3019 41 : anXDims[0] = nXDimID;
3020 : CPLDebug( "GDAL_netCDF", "nc_def_var(%d,%s,%d)",
3021 41 : cdfid, NCDF_DIMNAME_X, NC_DOUBLE );
3022 : status = nc_def_var( cdfid, NCDF_DIMNAME_X, NC_DOUBLE,
3023 41 : 1, anXDims, &NCDFVarID );
3024 41 : NCDF_ERR(status);
3025 41 : nVarXID=NCDFVarID;
3026 : nc_put_att_text( cdfid, NCDFVarID, CF_STD_NAME,
3027 : strlen(CF_PROJ_X_COORD),
3028 41 : CF_PROJ_X_COORD );
3029 : nc_put_att_text( cdfid, NCDFVarID, CF_LNG_NAME,
3030 : strlen("x coordinate of projection"),
3031 41 : "x coordinate of projection" );
3032 : /* TODO verify this is ok */
3033 41 : nc_put_att_text( cdfid, NCDFVarID, CF_UNITS, 1, "m" );
3034 :
3035 : /* Y */
3036 : int anYDims[1];
3037 41 : anYDims[0] = nYDimID;
3038 : CPLDebug( "GDAL_netCDF", "nc_def_var(%d,%s,%d)",
3039 41 : cdfid, NCDF_DIMNAME_Y, NC_DOUBLE );
3040 : status = nc_def_var( cdfid, NCDF_DIMNAME_Y, NC_DOUBLE,
3041 41 : 1, anYDims, &NCDFVarID );
3042 41 : NCDF_ERR(status);
3043 41 : nVarYID=NCDFVarID;
3044 : nc_put_att_text( cdfid, NCDFVarID, CF_STD_NAME,
3045 : strlen(CF_PROJ_Y_COORD),
3046 41 : CF_PROJ_Y_COORD );
3047 : nc_put_att_text( cdfid, NCDFVarID, CF_LNG_NAME,
3048 : strlen("y coordinate of projection"),
3049 41 : "y coordinate of projection" );
3050 41 : nc_put_att_text( cdfid, NCDFVarID, CF_UNITS, 1, "m" );
3051 : }
3052 :
3053 : /* -------------------------------------------------------------------- */
3054 : /* Write lat/lon attributes if needed */
3055 : /* -------------------------------------------------------------------- */
3056 85 : if ( bWriteLonLat == TRUE ) {
3057 :
3058 : /* latitude attributes */
3059 44 : if ( bIsProjected ) {
3060 : int anLatDims[2];
3061 0 : anLatDims[0] = nYDimID;
3062 0 : anLatDims[1] = nXDimID;
3063 : CPLDebug( "GDAL_netCDF", "nc_def_var(%d,%s,%d)",
3064 0 : cdfid, NCDF_DIMNAME_LAT, eLonLatType );
3065 : status = nc_def_var( cdfid, NCDF_DIMNAME_LAT, eLonLatType,
3066 0 : 2, anLatDims, &NCDFVarID );
3067 0 : NCDF_ERR(status);
3068 : /* compress lon/lat to save space */
3069 0 : DefVarDeflate( NCDFVarID );
3070 : }
3071 : else {
3072 : int anLatDims[1];
3073 : // anLatDims[0] = nLatDimID;
3074 44 : anLatDims[0] = nYDimID;
3075 : CPLDebug( "GDAL_netCDF", "nc_def_var(%d,%s,%d)",
3076 44 : cdfid, NCDF_DIMNAME_LAT, eLonLatType );
3077 : status = nc_def_var( cdfid, NCDF_DIMNAME_LAT, eLonLatType,
3078 44 : 1, anLatDims, &NCDFVarID );
3079 44 : NCDF_ERR(status);
3080 : }
3081 44 : nVarLatID = NCDFVarID;
3082 : nc_put_att_text( cdfid, NCDFVarID, CF_STD_NAME,
3083 44 : 8,"latitude" );
3084 : nc_put_att_text( cdfid, NCDFVarID, CF_LNG_NAME,
3085 44 : 8, "latitude" );
3086 : nc_put_att_text( cdfid, NCDFVarID, CF_UNITS,
3087 44 : 13, "degrees_north" );
3088 :
3089 : /* longitude attributes */
3090 44 : if ( bIsProjected ) {
3091 : int anLonDims[2];
3092 0 : anLonDims[0] = nYDimID;
3093 0 : anLonDims[1] = nXDimID;
3094 : CPLDebug( "GDAL_netCDF", "nc_def_var(%d,%s,%d)",
3095 0 : cdfid, NCDF_DIMNAME_LON, eLonLatType );
3096 : status = nc_def_var( cdfid, NCDF_DIMNAME_LON, eLonLatType,
3097 0 : 2, anLonDims, &NCDFVarID );
3098 0 : NCDF_ERR(status);
3099 : /* compress lon/lat to save space */
3100 0 : DefVarDeflate( NCDFVarID );
3101 : }
3102 : else {
3103 : int anLonDims[1];
3104 : // anLonDims[0] = nLonDimID;
3105 44 : anLonDims[0] = nXDimID;
3106 : CPLDebug( "GDAL_netCDF", "nc_def_var(%d,%s,%d)",
3107 44 : cdfid, NCDF_DIMNAME_LON, eLonLatType );
3108 : status = nc_def_var( cdfid, NCDF_DIMNAME_LON, eLonLatType,
3109 44 : 1, anLonDims, &NCDFVarID );
3110 44 : NCDF_ERR(status);
3111 : }
3112 44 : nVarLonID = NCDFVarID;
3113 : nc_put_att_text( cdfid, NCDFVarID, CF_STD_NAME,
3114 44 : 9, "longitude" );
3115 : nc_put_att_text( cdfid, NCDFVarID, CF_LNG_NAME,
3116 44 : 9, "longitude" );
3117 : nc_put_att_text( cdfid, NCDFVarID, CF_UNITS,
3118 44 : 12, "degrees_east" );
3119 : }
3120 :
3121 85 : pfnProgress( 0.50, NULL, pProgressData );
3122 :
3123 : /* -------------------------------------------------------------------- */
3124 : /* Get projection values */
3125 : /* -------------------------------------------------------------------- */
3126 :
3127 : double dfX0, dfDX, dfY0, dfDY;
3128 85 : dfX0=0.0, dfDX=0.0, dfY0=0.0, dfDY=0.0;
3129 85 : double *padLonVal = NULL;
3130 85 : double *padLatVal = NULL; /* should use float for projected, save space */
3131 :
3132 85 : if( bIsProjected )
3133 : {
3134 : // const char *pszProjection;
3135 41 : OGRSpatialReference oSRS;
3136 41 : OGRSpatialReference *poLatLonCRS = NULL;
3137 41 : OGRCoordinateTransformation *poTransform = NULL;
3138 :
3139 41 : char *pszWKT = (char *) pszProjection;
3140 41 : oSRS.importFromWkt( &pszWKT );
3141 :
3142 41 : double *padYVal = NULL;
3143 41 : double *padXVal = NULL;
3144 : size_t startX[1];
3145 : size_t countX[1];
3146 : size_t startY[1];
3147 : size_t countY[1];
3148 :
3149 41 : CPLDebug("GDAL_netCDF", "Getting (X,Y) values" );
3150 :
3151 41 : padXVal = (double *) CPLMalloc( nRasterXSize * sizeof( double ) );
3152 41 : padYVal = (double *) CPLMalloc( nRasterYSize * sizeof( double ) );
3153 :
3154 : /* -------------------------------------------------------------------- */
3155 : /* Get Y values */
3156 : /* -------------------------------------------------------------------- */
3157 41 : if ( ! bBottomUp )
3158 41 : dfY0 = adfGeoTransform[3];
3159 : else /* invert latitude values */
3160 0 : dfY0 = adfGeoTransform[3] + ( adfGeoTransform[5] * nRasterYSize );
3161 41 : dfDY = adfGeoTransform[5];
3162 :
3163 3015 : for( int j=0; j<nRasterYSize; j++ ) {
3164 : /* The data point is centered inside the pixel */
3165 2974 : if ( ! bBottomUp )
3166 2974 : padYVal[j] = dfY0 + (j+0.5)*dfDY ;
3167 : else /* invert latitude values */
3168 0 : padYVal[j] = dfY0 - (j+0.5)*dfDY ;
3169 : }
3170 41 : startX[0] = 0;
3171 41 : countX[0] = nRasterXSize;
3172 :
3173 : /* -------------------------------------------------------------------- */
3174 : /* Get X values */
3175 : /* -------------------------------------------------------------------- */
3176 41 : dfX0 = adfGeoTransform[0];
3177 41 : dfDX = adfGeoTransform[1];
3178 :
3179 6585 : for( int i=0; i<nRasterXSize; i++ ) {
3180 : /* The data point is centered inside the pixel */
3181 6544 : padXVal[i] = dfX0 + (i+0.5)*dfDX ;
3182 : }
3183 41 : startY[0] = 0;
3184 41 : countY[0] = nRasterYSize;
3185 :
3186 : /* -------------------------------------------------------------------- */
3187 : /* Write X/Y values */
3188 : /* -------------------------------------------------------------------- */
3189 : /* make sure we are in data mode */
3190 41 : SetDefineMode( FALSE );
3191 :
3192 41 : CPLDebug("GDAL_netCDF", "Writing X values" );
3193 : status = nc_put_vara_double( cdfid, nVarXID, startX,
3194 41 : countX, padXVal);
3195 41 : NCDF_ERR(status);
3196 :
3197 41 : CPLDebug("GDAL_netCDF", "Writing Y values" );
3198 : status = nc_put_vara_double( cdfid, nVarYID, startY,
3199 41 : countY, padYVal);
3200 41 : NCDF_ERR(status);
3201 :
3202 41 : pfnProgress( 0.20, NULL, pProgressData );
3203 :
3204 : /* -------------------------------------------------------------------- */
3205 : /* Transform (X,Y) values to (lon,lat) */
3206 : /* -------------------------------------------------------------------- */
3207 :
3208 : /* Get OGR transform */
3209 41 : if ( bWriteLonLat == TRUE ) {
3210 0 : poLatLonCRS = oSRS.CloneGeogCS();
3211 0 : if ( poLatLonCRS != NULL )
3212 0 : poTransform = OGRCreateCoordinateTransformation( &oSRS, poLatLonCRS );
3213 : /* if no OGR transform, then don't write CF lon/lat */
3214 0 : if( poTransform == NULL ) {
3215 : CPLError( CE_Failure, CPLE_AppDefined,
3216 0 : "Unable to get Coordinate Transform" );
3217 0 : bWriteLonLat = FALSE;
3218 : }
3219 : }
3220 :
3221 41 : if ( bWriteLonLat == TRUE ) {
3222 :
3223 0 : CPLDebug("GDAL_netCDF", "Transforming (X,Y)->(lon,lat)" );
3224 :
3225 0 : padLatVal = (double *) CPLMalloc( nRasterXSize * sizeof( double ) );
3226 0 : padLonVal = (double *) CPLMalloc( nRasterXSize * sizeof( double ) );
3227 : size_t start[2], count[2];
3228 0 : start[1] = 0; //X
3229 0 : count[1] = nRasterXSize;
3230 0 : start[0] = 0; //Y
3231 0 : count[0] = 1;
3232 :
3233 0 : int bOK = TRUE;
3234 : double dfTemp;
3235 : int i,j;
3236 0 : dfTemp = 0.2;
3237 :
3238 0 : for( j = 0; (j < nRasterYSize) && bOK && (status == NC_NOERR); j++ ) {
3239 :
3240 0 : start[0] = j;
3241 :
3242 : /* fill values to transform */
3243 0 : for( i=0; i<nRasterXSize; i++ ) {
3244 0 : padLatVal[i] = padYVal[j];
3245 0 : padLonVal[i] = padXVal[i];
3246 : }
3247 :
3248 : /* do the transform */
3249 : bOK = poTransform->Transform( nRasterXSize,
3250 0 : padLonVal, padLatVal, NULL );
3251 : /* write data */
3252 0 : if ( bOK ) {
3253 : status = nc_put_vara_double( cdfid, nVarLatID, start,
3254 0 : count, padLatVal);
3255 0 : NCDF_ERR(status);
3256 : status = nc_put_vara_double( cdfid, nVarLonID, start,
3257 0 : count, padLonVal);
3258 0 : NCDF_ERR(status);
3259 : }
3260 : else
3261 : CPLError( CE_Failure, CPLE_AppDefined,
3262 0 : "Unable to Transform (X,Y) to (lon,lat).\n" );
3263 :
3264 0 : if ( j % (nRasterYSize/10) == 0 ) {
3265 0 : dfTemp += 0.08;
3266 0 : pfnProgress( dfTemp , NULL, pProgressData );
3267 : }
3268 : }
3269 :
3270 : }
3271 :
3272 : /* Free the srs and transform objects */
3273 41 : if ( poLatLonCRS != NULL ) CPLFree( poLatLonCRS );
3274 41 : if ( poTransform != NULL ) CPLFree( poTransform );
3275 :
3276 : /* Free data */
3277 41 : CPLFree( padXVal );
3278 41 : CPLFree( padYVal );
3279 41 : CPLFree( padLonVal );
3280 41 : CPLFree( padLatVal);
3281 :
3282 : } // projected
3283 :
3284 : else { /* If not Projected assume Geographic to catch grids without Datum */
3285 :
3286 : /* -------------------------------------------------------------------- */
3287 : /* Get latitude values */
3288 : /* -------------------------------------------------------------------- */
3289 44 : if ( ! bBottomUp )
3290 44 : dfY0 = adfGeoTransform[3];
3291 : else /* invert latitude values */
3292 0 : dfY0 = adfGeoTransform[3] + ( adfGeoTransform[5] * nRasterYSize );
3293 44 : dfDY = adfGeoTransform[5];
3294 :
3295 44 : padLatVal = (double *) CPLMalloc( nRasterYSize * sizeof( double ) );
3296 2977 : for( int i=0; i<nRasterYSize; i++ ) {
3297 : /* The data point is centered inside the pixel */
3298 2933 : if ( ! bBottomUp )
3299 2933 : padLatVal[i] = dfY0 + (i+0.5)*dfDY ;
3300 : else /* invert latitude values */
3301 0 : padLatVal[i] = dfY0 - (i+0.5)*dfDY ;
3302 : }
3303 :
3304 : size_t startLat[1];
3305 : size_t countLat[1];
3306 44 : startLat[0] = 0;
3307 44 : countLat[0] = nRasterYSize;
3308 :
3309 : /* -------------------------------------------------------------------- */
3310 : /* Get longitude values */
3311 : /* -------------------------------------------------------------------- */
3312 44 : dfX0 = adfGeoTransform[0];
3313 44 : dfDX = adfGeoTransform[1];
3314 :
3315 44 : padLonVal = (double *) CPLMalloc( nRasterXSize * sizeof( double ) );
3316 2977 : for( int i=0; i<nRasterXSize; i++ ) {
3317 : /* The data point is centered inside the pixel */
3318 2933 : padLonVal[i] = dfX0 + (i+0.5)*dfDX ;
3319 : }
3320 :
3321 : size_t startLon[1];
3322 : size_t countLon[1];
3323 44 : startLon[0] = 0;
3324 44 : countLon[0] = nRasterXSize;
3325 :
3326 : /* -------------------------------------------------------------------- */
3327 : /* Write latitude and longitude values */
3328 : /* -------------------------------------------------------------------- */
3329 : /* latitude values */
3330 44 : CPLDebug("GDAL_netCDF", "Writing lat values" );
3331 :
3332 : /* make sure we are in data mode */
3333 44 : SetDefineMode( FALSE );
3334 :
3335 : status = nc_put_vara_double( cdfid, nVarLatID, startLat,
3336 44 : countLat, padLatVal);
3337 44 : NCDF_ERR(status);
3338 :
3339 : /* free values */
3340 44 : CPLFree( padLatVal );
3341 :
3342 : /* longitude values */
3343 44 : CPLDebug("GDAL_netCDF", "Writing lon values" );
3344 : status = nc_put_vara_double( cdfid, nVarLonID, startLon,
3345 44 : countLon, padLonVal);
3346 44 : NCDF_ERR(status);
3347 :
3348 : /* free values */
3349 44 : CPLFree( padLonVal );
3350 :
3351 : }// not projected
3352 :
3353 :
3354 85 : pfnProgress( 1.00, NULL, pProgressData );
3355 :
3356 85 : return CE_None;
3357 : }
3358 :
3359 : /************************************************************************/
3360 : /* GetGeoTransform() */
3361 : /************************************************************************/
3362 :
3363 61 : CPLErr netCDFDataset::GetGeoTransform( double * padfTransform )
3364 :
3365 : {
3366 61 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
3367 61 : if( bSetGeoTransform )
3368 61 : return CE_None;
3369 : else
3370 0 : return GDALPamDataset::GetGeoTransform( padfTransform );
3371 : }
3372 :
3373 : /************************************************************************/
3374 : /* rint() */
3375 : /************************************************************************/
3376 :
3377 837 : double netCDFDataset::rint( double dfX)
3378 : {
3379 837 : if( dfX > 0 ) {
3380 402 : int nX = (int) (dfX+0.5);
3381 402 : if( nX % 2 ) {
3382 10 : double dfDiff = dfX - (double)nX;
3383 10 : if( dfDiff == -0.5 )
3384 0 : return double( nX-1 );
3385 : }
3386 402 : return double( nX );
3387 : } else {
3388 435 : int nX= (int) (dfX-0.5);
3389 435 : if( nX % 2 ) {
3390 21 : double dfDiff = dfX - (double)nX;
3391 21 : if( dfDiff == 0.5 )
3392 0 : return double(nX+1);
3393 : }
3394 435 : return double(nX);
3395 : }
3396 : }
3397 :
3398 : /************************************************************************/
3399 : /* ReadAttributes() */
3400 : /************************************************************************/
3401 757 : CPLErr netCDFDataset::ReadAttributes( int cdfid, int var)
3402 :
3403 : {
3404 : char szAttrName[ NC_MAX_NAME ];
3405 : char szVarName [ NC_MAX_NAME ];
3406 : char szMetaName[ NC_MAX_NAME * 2 ];
3407 757 : char *pszMetaTemp = NULL;
3408 : int nbAttr;
3409 :
3410 757 : nc_inq_varnatts( cdfid, var, &nbAttr );
3411 757 : if( var == NC_GLOBAL ) {
3412 185 : strcpy( szVarName,"NC_GLOBAL" );
3413 : }
3414 : else {
3415 572 : nc_inq_varname( cdfid, var, szVarName );
3416 : }
3417 :
3418 4431 : for( int l=0; l < nbAttr; l++) {
3419 :
3420 3674 : nc_inq_attname( cdfid, var, l, szAttrName);
3421 3674 : sprintf( szMetaName, "%s#%s", szVarName, szAttrName );
3422 :
3423 3674 : if ( NCDFGetAttr( cdfid, var, szAttrName, &pszMetaTemp )
3424 : == CE_None ) {
3425 : papszMetadata = CSLSetNameValue(papszMetadata,
3426 : szMetaName,
3427 3674 : pszMetaTemp);
3428 3674 : CPLFree(pszMetaTemp);
3429 3674 : pszMetaTemp = NULL;
3430 : }
3431 : else {
3432 0 : CPLDebug( "GDAL_netCDF", "invalid global metadata %s", szMetaName );
3433 : }
3434 :
3435 : }
3436 :
3437 757 : return CE_None;
3438 :
3439 : }
3440 :
3441 :
3442 : /************************************************************************/
3443 : /* netCDFDataset::CreateSubDatasetList() */
3444 : /************************************************************************/
3445 18 : void netCDFDataset::CreateSubDatasetList( )
3446 : {
3447 :
3448 : char szDim[ MAX_NC_NAME ];
3449 : char szTemp[ MAX_NC_NAME ];
3450 : char szType[ MAX_NC_NAME ];
3451 : char szName[ MAX_NC_NAME ];
3452 : char szVarStdName[ MAX_NC_NAME ];
3453 : int nDims;
3454 : int nVar;
3455 : int nVarCount;
3456 : int i;
3457 : nc_type nVarType;
3458 : int *ponDimIds;
3459 : size_t nDimLen;
3460 : int nSub;
3461 : nc_type nAttype;
3462 : size_t nAttlen;
3463 :
3464 : netCDFDataset *poDS;
3465 18 : poDS = this;
3466 :
3467 18 : nSub=1;
3468 18 : nc_inq_nvars ( cdfid, &nVarCount );
3469 130 : for ( nVar = 0; nVar < nVarCount; nVar++ ) {
3470 :
3471 112 : nc_inq_varndims ( cdfid, nVar, &nDims );
3472 112 : if( nDims >= 2 ) {
3473 58 : ponDimIds = (int *) CPLCalloc( nDims, sizeof( int ) );
3474 58 : nc_inq_vardimid ( cdfid, nVar, ponDimIds );
3475 :
3476 : /* -------------------------------------------------------------------- */
3477 : /* Create Sub dataset list */
3478 : /* -------------------------------------------------------------------- */
3479 58 : szDim[0]='\0';
3480 174 : for( i = 0; i < nDims; i++ ) {
3481 116 : nc_inq_dimlen ( cdfid, ponDimIds[i], &nDimLen );
3482 116 : sprintf(szTemp, "%d", (int) nDimLen);
3483 116 : strcat(szTemp, "x" );
3484 116 : strcat(szDim, szTemp);
3485 : }
3486 :
3487 58 : nc_inq_vartype( cdfid, nVar, &nVarType );
3488 : /* -------------------------------------------------------------------- */
3489 : /* Get rid of the last "x" character */
3490 : /* -------------------------------------------------------------------- */
3491 58 : szDim[strlen(szDim) - 1] = '\0';
3492 58 : switch( nVarType ) {
3493 :
3494 : case NC_BYTE:
3495 : #ifdef NETCDF_HAS_NC4
3496 : case NC_UBYTE:
3497 : #endif
3498 : case NC_CHAR:
3499 28 : strcpy(szType, "8-bit character");
3500 28 : break;
3501 :
3502 : case NC_SHORT:
3503 3 : strcpy(szType, "8-bit integer");
3504 3 : break;
3505 : case NC_INT:
3506 3 : strcpy(szType, "16-bit integer");
3507 3 : break;
3508 : case NC_FLOAT:
3509 21 : strcpy(szType, "32-bit floating-point");
3510 21 : break;
3511 : case NC_DOUBLE:
3512 3 : strcpy(szType, "64-bit floating-point");
3513 : break;
3514 :
3515 : default:
3516 : break;
3517 : }
3518 58 : nc_inq_varname( cdfid, nVar, szName);
3519 58 : nc_inq_att( cdfid, nVar, CF_STD_NAME, &nAttype, &nAttlen);
3520 58 : if( nc_get_att_text ( cdfid, nVar, CF_STD_NAME,
3521 : szVarStdName ) == NC_NOERR ) {
3522 0 : szVarStdName[nAttlen] = '\0';
3523 : }
3524 : else {
3525 58 : strcpy( szVarStdName, szName );
3526 : }
3527 :
3528 58 : sprintf( szTemp, "SUBDATASET_%d_NAME", nSub) ;
3529 :
3530 : poDS->papszSubDatasets =
3531 : CSLSetNameValue( poDS->papszSubDatasets, szTemp,
3532 : CPLSPrintf( "NETCDF:\"%s\":%s",
3533 : poDS->osFilename.c_str(),
3534 58 : szName) ) ;
3535 :
3536 58 : sprintf( szTemp, "SUBDATASET_%d_DESC", nSub++ );
3537 :
3538 : poDS->papszSubDatasets =
3539 : CSLSetNameValue( poDS->papszSubDatasets, szTemp,
3540 : CPLSPrintf( "[%s] %s (%s)",
3541 : szDim,
3542 : szVarStdName,
3543 58 : szType ) );
3544 58 : CPLFree(ponDimIds);
3545 : }
3546 : }
3547 :
3548 18 : }
3549 :
3550 : /************************************************************************/
3551 : /* IdentifyFormat() */
3552 : /************************************************************************/
3553 :
3554 12247 : int netCDFDataset::IdentifyFormat( GDALOpenInfo * poOpenInfo, bool bCheckExt = TRUE )
3555 :
3556 : {
3557 : /* -------------------------------------------------------------------- */
3558 : /* Does this appear to be a netcdf file? If so, which format? */
3559 : /* http://www.unidata.ucar.edu/software/netcdf/docs/faq.html#fv1_5 */
3560 : /* -------------------------------------------------------------------- */
3561 :
3562 12247 : if( EQUALN(poOpenInfo->pszFilename,"NETCDF:",7) )
3563 0 : return NCDF_FORMAT_UNKNOWN;
3564 12247 : if ( poOpenInfo->nHeaderBytes < 4 )
3565 10860 : return NCDF_FORMAT_NONE;
3566 1387 : if ( EQUALN((char*)poOpenInfo->pabyHeader,"CDF\001",4) )
3567 184 : return NCDF_FORMAT_NC;
3568 1203 : else if ( EQUALN((char*)poOpenInfo->pabyHeader,"CDF\002",4) )
3569 1 : return NCDF_FORMAT_NC2;
3570 1202 : else if ( EQUALN((char*)poOpenInfo->pabyHeader,"\211HDF\r\n\032\n",8) ) {
3571 : /* Requires netCDF-4/HDF5 support in libnetcdf (not just libnetcdf-v4).
3572 : If HDF5 is not supported in GDAL, this driver will try to open the file
3573 : Else, make sure this driver does not try to open HDF5 files
3574 : If user really wants to open with this driver, use NETCDF:file.h5 format.
3575 : This check should be relaxed, but there is no clear way to make a difference.
3576 : */
3577 :
3578 : /* Check for HDF5 support in GDAL */
3579 : #ifdef HAVE_HDF5
3580 5 : if ( bCheckExt ) { /* Check by default */
3581 5 : const char* pszExtension = CPLGetExtension( poOpenInfo->pszFilename );
3582 5 : if ( ! ( EQUAL( pszExtension, "nc") || EQUAL( pszExtension, "cdf")
3583 : || EQUAL( pszExtension, "nc2") || EQUAL( pszExtension, "nc4") ) )
3584 5 : return NCDF_FORMAT_HDF5;
3585 : }
3586 : #endif
3587 :
3588 : /* Check for netcdf-4 support in libnetcdf */
3589 : #ifdef NETCDF_HAS_NC4
3590 : return NCDF_FORMAT_NC4;
3591 : #else
3592 0 : return NCDF_FORMAT_HDF5;
3593 : #endif
3594 :
3595 : }
3596 1197 : else if ( EQUALN((char*)poOpenInfo->pabyHeader,"\016\003\023\001",4) ) {
3597 : /* Requires HDF4 support in libnetcdf, but if HF4 is supported by GDAL don't try to open. */
3598 : /* If user really wants to open with this driver, use NETCDF:file.hdf syntax. */
3599 :
3600 : /* Check for HDF4 support in GDAL */
3601 : #ifdef HAVE_HDF4
3602 257 : if ( bCheckExt ) { /* Check by default */
3603 : /* Always treat as HDF4 file */
3604 257 : return NCDF_FORMAT_HDF4;
3605 : }
3606 : #endif
3607 :
3608 : /* Check for HDF4 support in libnetcdf */
3609 : #ifdef NETCDF_HAS_HDF4
3610 : return NCDF_FORMAT_NC4;
3611 : #else
3612 0 : return NCDF_FORMAT_HDF4;
3613 : #endif
3614 : }
3615 :
3616 940 : return NCDF_FORMAT_NONE;
3617 : }
3618 :
3619 : /************************************************************************/
3620 : /* Identify() */
3621 : /************************************************************************/
3622 :
3623 9246 : int netCDFDataset::Identify( GDALOpenInfo * poOpenInfo )
3624 :
3625 : {
3626 9246 : if( EQUALN(poOpenInfo->pszFilename,"NETCDF:",7) ) {
3627 0 : return TRUE;
3628 : }
3629 9246 : int nTmpFormat = IdentifyFormat( poOpenInfo );
3630 9246 : if( NCDF_FORMAT_NC == nTmpFormat ||
3631 : NCDF_FORMAT_NC2 == nTmpFormat ||
3632 : NCDF_FORMAT_NC4 == nTmpFormat ||
3633 : NCDF_FORMAT_NC4C == nTmpFormat )
3634 0 : return TRUE;
3635 : else
3636 9246 : return FALSE;
3637 : }
3638 :
3639 : /************************************************************************/
3640 : /* Open() */
3641 : /************************************************************************/
3642 :
3643 3001 : GDALDataset *netCDFDataset::Open( GDALOpenInfo * poOpenInfo )
3644 :
3645 : {
3646 : int j;
3647 : unsigned int k;
3648 : int nd;
3649 : int cdfid, dim_count, var, var_count;
3650 3001 : int i = 0;
3651 : size_t lev_count;
3652 3001 : size_t nTotLevCount = 1;
3653 3001 : int nDim = 2;
3654 : int status;
3655 : int nDimID;
3656 : char attname[NC_MAX_NAME];
3657 : int ndims, nvars, ngatts, unlimdimid;
3658 3001 : int nCount=0;
3659 3001 : int nVarID=-1;
3660 :
3661 3001 : int nTmpFormat=NCDF_FORMAT_NONE;
3662 : int *panBandDimPos; // X, Y, Z postion in array
3663 : int *panBandZLev;
3664 : int *paDimIds;
3665 : size_t xdim, ydim;
3666 : char szTemp[NC_MAX_NAME];
3667 :
3668 3001 : CPLString osSubdatasetName;
3669 : int bTreatAsSubdataset;
3670 :
3671 : /* -------------------------------------------------------------------- */
3672 : /* Does this appear to be a netcdf file? */
3673 : /* -------------------------------------------------------------------- */
3674 3001 : if( ! EQUALN(poOpenInfo->pszFilename,"NETCDF:",7) ) {
3675 2996 : nTmpFormat = IdentifyFormat( poOpenInfo );
3676 : /* Note: not calling Identify() directly, because we want the file type */
3677 : /* Only support NCDF_FORMAT* formats */
3678 2996 : if( ! ( NCDF_FORMAT_NC == nTmpFormat ||
3679 : NCDF_FORMAT_NC2 == nTmpFormat ||
3680 : NCDF_FORMAT_NC4 == nTmpFormat ||
3681 : NCDF_FORMAT_NC4C == nTmpFormat ) )
3682 2816 : return NULL;
3683 : }
3684 :
3685 : netCDFDataset *poDS;
3686 185 : poDS = new netCDFDataset();
3687 :
3688 : /* -------------------------------------------------------------------- */
3689 : /* Disable PAM, at least temporarily. See bug #4244 */
3690 : /* -------------------------------------------------------------------- */
3691 185 : poDS->nPamFlags |= GPF_DISABLED;
3692 :
3693 :
3694 185 : poDS->SetDescription( poOpenInfo->pszFilename );
3695 :
3696 : /* -------------------------------------------------------------------- */
3697 : /* Check if filename start with NETCDF: tag */
3698 : /* -------------------------------------------------------------------- */
3699 185 : if( EQUALN( poOpenInfo->pszFilename,"NETCDF:",7) )
3700 : {
3701 : char **papszName =
3702 : CSLTokenizeString2( poOpenInfo->pszFilename,
3703 5 : ":", CSLT_HONOURSTRINGS|CSLT_PRESERVEESCAPES );
3704 :
3705 : /* -------------------------------------------------------------------- */
3706 : /* Check for drive name in windows NETCDF:"D:\... */
3707 : /* -------------------------------------------------------------------- */
3708 5 : if ( CSLCount(papszName) == 4 &&
3709 0 : strlen(papszName[1]) == 1 &&
3710 0 : (papszName[2][0] == '/' || papszName[2][0] == '\\') )
3711 : {
3712 0 : poDS->osFilename = papszName[1];
3713 0 : poDS->osFilename += ':';
3714 0 : poDS->osFilename += papszName[2];
3715 0 : osSubdatasetName = papszName[3];
3716 0 : bTreatAsSubdataset = TRUE;
3717 0 : CSLDestroy( papszName );
3718 : }
3719 5 : else if( CSLCount(papszName) == 3 )
3720 : {
3721 5 : poDS->osFilename = papszName[1];
3722 5 : osSubdatasetName = papszName[2];
3723 5 : bTreatAsSubdataset = TRUE;
3724 5 : CSLDestroy( papszName );
3725 : }
3726 0 : else if( CSLCount(papszName) == 2 )
3727 : {
3728 0 : poDS->osFilename = papszName[1];
3729 0 : osSubdatasetName = "";
3730 0 : bTreatAsSubdataset = FALSE;
3731 0 : CSLDestroy( papszName );
3732 : }
3733 : else
3734 : {
3735 0 : CSLDestroy( papszName );
3736 0 : delete poDS;
3737 : CPLError( CE_Failure, CPLE_AppDefined,
3738 0 : "Failed to parse NETCDF: prefix string into expected 2, 3 or 4 fields." );
3739 0 : return NULL;
3740 : }
3741 : /* Identify Format from real file, with bCheckExt=FALSE */
3742 5 : GDALOpenInfo* poOpenInfo2 = new GDALOpenInfo(poDS->osFilename.c_str(), GA_ReadOnly );
3743 10 : poDS->nFormat = IdentifyFormat( poOpenInfo2, FALSE );
3744 5 : delete poOpenInfo2;
3745 5 : if( NCDF_FORMAT_NONE == poDS->nFormat ||
3746 : NCDF_FORMAT_UNKNOWN == poDS->nFormat ) {
3747 0 : delete poDS;
3748 0 : return NULL;
3749 : }
3750 : }
3751 : else
3752 : {
3753 180 : poDS->osFilename = poOpenInfo->pszFilename;
3754 180 : bTreatAsSubdataset = FALSE;
3755 180 : poDS->nFormat = nTmpFormat;
3756 : }
3757 :
3758 : /* -------------------------------------------------------------------- */
3759 : /* Try opening the dataset. */
3760 : /* -------------------------------------------------------------------- */
3761 185 : CPLDebug( "GDAL_netCDF", "\n=====\ncalling nc_open( %s )\n", poDS->osFilename.c_str() );
3762 185 : if( nc_open( poDS->osFilename, NC_NOWRITE, &cdfid ) != NC_NOERR ) {
3763 0 : delete poDS;
3764 0 : return NULL;
3765 : }
3766 :
3767 : /* -------------------------------------------------------------------- */
3768 : /* Is this a real netCDF file? */
3769 : /* -------------------------------------------------------------------- */
3770 185 : status = nc_inq(cdfid, &ndims, &nvars, &ngatts, &unlimdimid);
3771 185 : if( status != NC_NOERR ) {
3772 0 : delete poDS;
3773 0 : return NULL;
3774 : }
3775 :
3776 : /* -------------------------------------------------------------------- */
3777 : /* Get file type from netcdf */
3778 : /* -------------------------------------------------------------------- */
3779 185 : status = nc_inq_format (cdfid, &nTmpFormat);
3780 185 : if ( status != NC_NOERR ) {
3781 0 : NCDF_ERR(status);
3782 : }
3783 : else {
3784 : CPLDebug( "GDAL_netCDF",
3785 : "driver detected file type=%d, libnetcdf detected type=%d",
3786 185 : poDS->nFormat, nTmpFormat );
3787 185 : if ( nTmpFormat != poDS->nFormat ) {
3788 : /* warn if file detection conflicts with that from libnetcdf */
3789 : /* except for NC4C, which we have no way of detecting initially */
3790 0 : if ( nTmpFormat != NCDF_FORMAT_NC4C ) {
3791 : CPLError( CE_Warning, CPLE_AppDefined,
3792 : "NetCDF driver detected file type=%d, but libnetcdf detected type=%d",
3793 0 : poDS->nFormat, nTmpFormat );
3794 : }
3795 : CPLDebug( "GDAL_netCDF", "seting file type to %d, was %d",
3796 0 : nTmpFormat, poDS->nFormat );
3797 0 : poDS->nFormat = nTmpFormat;
3798 : }
3799 : }
3800 :
3801 : /* -------------------------------------------------------------------- */
3802 : /* Confirm the requested access is supported. */
3803 : /* -------------------------------------------------------------------- */
3804 185 : if( poOpenInfo->eAccess == GA_Update )
3805 : {
3806 : CPLError( CE_Failure, CPLE_NotSupported,
3807 : "The NETCDF driver does not support update access to existing"
3808 0 : " datasets.\n" );
3809 0 : nc_close( cdfid );
3810 0 : delete poDS;
3811 0 : return NULL;
3812 : }
3813 :
3814 : /* -------------------------------------------------------------------- */
3815 : /* Does the request variable exist? */
3816 : /* -------------------------------------------------------------------- */
3817 185 : if( bTreatAsSubdataset )
3818 : {
3819 5 : status = nc_inq_varid( cdfid, osSubdatasetName, &var);
3820 5 : if( status != NC_NOERR ) {
3821 : CPLError( CE_Warning, CPLE_AppDefined,
3822 : "%s is a netCDF file, but %s is not a variable.",
3823 : poOpenInfo->pszFilename,
3824 0 : osSubdatasetName.c_str() );
3825 :
3826 0 : nc_close( cdfid );
3827 0 : delete poDS;
3828 0 : return NULL;
3829 : }
3830 : }
3831 :
3832 185 : if( nc_inq_ndims( cdfid, &dim_count ) != NC_NOERR || dim_count < 2 )
3833 : {
3834 : CPLError( CE_Warning, CPLE_AppDefined,
3835 : "%s is a netCDF file, but not in GMT configuration.",
3836 0 : poOpenInfo->pszFilename );
3837 :
3838 0 : nc_close( cdfid );
3839 0 : delete poDS;
3840 0 : return NULL;
3841 : }
3842 :
3843 185 : CPLDebug( "GDAL_netCDF", "dim_count = %d", dim_count );
3844 :
3845 185 : if( (status = nc_get_att_text( cdfid, NC_GLOBAL, "Conventions",
3846 : attname )) != NC_NOERR ) {
3847 : CPLError( CE_Warning, CPLE_AppDefined,
3848 3 : "No UNIDATA NC_GLOBAL:Conventions attribute");
3849 : /* note that 'Conventions' is always capital 'C' in CF spec*/
3850 : }
3851 :
3852 :
3853 : /* -------------------------------------------------------------------- */
3854 : /* Create band information objects. */
3855 : /* -------------------------------------------------------------------- */
3856 185 : if ( nc_inq_nvars ( cdfid, &var_count) != NC_NOERR )
3857 : {
3858 0 : delete poDS;
3859 0 : return NULL;
3860 : }
3861 :
3862 185 : CPLDebug( "GDAL_netCDF", "var_count = %d", var_count );
3863 :
3864 : /* -------------------------------------------------------------------- */
3865 : /* Create a corresponding GDALDataset. */
3866 : /* Create Netcdf Subdataset if filename as NETCDF tag */
3867 : /* -------------------------------------------------------------------- */
3868 185 : poDS->cdfid = cdfid;
3869 :
3870 185 : poDS->ReadAttributes( cdfid, NC_GLOBAL );
3871 :
3872 : /* -------------------------------------------------------------------- */
3873 : /* Verify if only one variable has 2 dimensions */
3874 : /* -------------------------------------------------------------------- */
3875 876 : for ( j = 0; j < nvars; j++ ) {
3876 :
3877 691 : nc_inq_varndims ( cdfid, j, &ndims );
3878 691 : if( ndims >= 2 ) {
3879 228 : nVarID=j;
3880 228 : nCount++;
3881 : }
3882 : }
3883 :
3884 : /* -------------------------------------------------------------------- */
3885 : /* We have more than one variable with 2 dimensions in the */
3886 : /* file, then treat this as a subdataset container dataset. */
3887 : /* -------------------------------------------------------------------- */
3888 185 : if( (nCount > 1) && !bTreatAsSubdataset )
3889 : {
3890 18 : poDS->CreateSubDatasetList();
3891 18 : poDS->SetMetadata( poDS->papszMetadata );
3892 18 : poDS->TryLoadXML();
3893 18 : return( poDS );
3894 : }
3895 :
3896 : /* -------------------------------------------------------------------- */
3897 : /* If we are not treating things as a subdataset, then capture */
3898 : /* the name of the single available variable as the subdataset. */
3899 : /* -------------------------------------------------------------------- */
3900 167 : if( !bTreatAsSubdataset ) // nCount must be 1!
3901 : {
3902 : char szVarName[NC_MAX_NAME];
3903 :
3904 162 : nc_inq_varname( cdfid, nVarID, szVarName);
3905 162 : osSubdatasetName = szVarName;
3906 : }
3907 :
3908 : /* -------------------------------------------------------------------- */
3909 : /* Open the NETCDF subdataset NETCDF:"filename":subdataset */
3910 : /* -------------------------------------------------------------------- */
3911 167 : var=-1;
3912 167 : nc_inq_varid( cdfid, osSubdatasetName, &var);
3913 167 : nd = 0;
3914 167 : nc_inq_varndims ( cdfid, var, &nd );
3915 :
3916 167 : paDimIds = (int *)CPLCalloc(nd, sizeof( int ) );
3917 167 : panBandDimPos = ( int * ) CPLCalloc( nd, sizeof( int ) );
3918 :
3919 167 : nc_inq_vardimid( cdfid, var, paDimIds );
3920 :
3921 : /* -------------------------------------------------------------------- */
3922 : /* Check fi somebody tried to pass a variable with less than 2D */
3923 : /* -------------------------------------------------------------------- */
3924 167 : if ( nd < 2 ) {
3925 1 : CPLFree( paDimIds );
3926 1 : CPLFree( panBandDimPos );
3927 1 : delete poDS;
3928 1 : return NULL;
3929 : }
3930 :
3931 : /* -------------------------------------------------------------------- */
3932 : /* CF-1 Convention */
3933 : /* dimensions to appear in the relative order T, then Z, then Y, */
3934 : /* then X to the file. All other dimensions should, whenever */
3935 : /* possible, be placed to the left of the spatiotemporal */
3936 : /* dimensions. */
3937 : /* -------------------------------------------------------------------- */
3938 :
3939 : /* -------------------------------------------------------------------- */
3940 : /* Get X dimensions information */
3941 : /* -------------------------------------------------------------------- */
3942 166 : poDS->nXDimID = paDimIds[nd-1];
3943 166 : nc_inq_dimlen ( cdfid, poDS->nXDimID, &xdim );
3944 166 : poDS->nRasterXSize = xdim;
3945 :
3946 : /* -------------------------------------------------------------------- */
3947 : /* Get Y dimension information */
3948 : /* -------------------------------------------------------------------- */
3949 166 : poDS->nYDimID = paDimIds[nd-2];
3950 166 : nc_inq_dimlen ( cdfid, poDS->nYDimID, &ydim );
3951 166 : poDS->nRasterYSize = ydim;
3952 :
3953 :
3954 512 : for( j=0,k=0; j < nd; j++ ){
3955 346 : if( paDimIds[j] == poDS->nXDimID ){
3956 166 : panBandDimPos[0] = j; // Save Position of XDim
3957 166 : k++;
3958 : }
3959 346 : if( paDimIds[j] == poDS->nYDimID ){
3960 166 : panBandDimPos[1] = j; // Save Position of YDim
3961 166 : k++;
3962 : }
3963 : }
3964 : /* -------------------------------------------------------------------- */
3965 : /* X and Y Dimension Ids were not found! */
3966 : /* -------------------------------------------------------------------- */
3967 166 : if( k != 2 ) {
3968 0 : CPLFree( paDimIds );
3969 0 : CPLFree( panBandDimPos );
3970 0 : return NULL;
3971 : }
3972 :
3973 : /* -------------------------------------------------------------------- */
3974 : /* Read Metadata for this variable */
3975 : /* -------------------------------------------------------------------- */
3976 : /* should disable as is also done at band level, except driver needs the
3977 : variables as metadata (e.g. projection) */
3978 166 : poDS->ReadAttributes( cdfid, var );
3979 :
3980 : /* -------------------------------------------------------------------- */
3981 : /* Read Metadata for each dimension */
3982 : /* -------------------------------------------------------------------- */
3983 :
3984 512 : for( j=0; j < dim_count; j++ ){
3985 346 : nc_inq_dimname( cdfid, j, szTemp );
3986 346 : poDS->papszDimName.AddString( szTemp );
3987 346 : status = nc_inq_varid( cdfid, poDS->papszDimName[j], &nDimID );
3988 346 : if( status == NC_NOERR ) {
3989 290 : poDS->ReadAttributes( cdfid, nDimID );
3990 : }
3991 : }
3992 :
3993 166 : poDS->SetProjectionFromVar( var );
3994 166 : poDS->SetMetadata( poDS->papszMetadata );
3995 :
3996 : /* -------------------------------------------------------------------- */
3997 : /* Create bands */
3998 : /* -------------------------------------------------------------------- */
3999 166 : panBandZLev = (int *)CPLCalloc( nd-2, sizeof( int ) );
4000 :
4001 166 : nTotLevCount = 1;
4002 166 : if ( dim_count > 2 ) {
4003 10 : nDim=2;
4004 44 : for( j=0; j < nd; j++ ){
4005 58 : if( ( paDimIds[j] != poDS->nXDimID ) &&
4006 24 : ( paDimIds[j] != poDS->nYDimID ) ){
4007 14 : nc_inq_dimlen ( cdfid, paDimIds[j], &lev_count );
4008 14 : nTotLevCount *= lev_count;
4009 14 : panBandZLev[ nDim-2 ] = lev_count;
4010 14 : panBandDimPos[ nDim++ ] = j; //Save Position of ZDim
4011 : }
4012 : }
4013 : }
4014 166 : i=0;
4015 :
4016 346 : for ( unsigned int lev = 0; lev < nTotLevCount ; lev++ ) {
4017 : char ** papszToken;
4018 180 : papszToken=NULL;
4019 :
4020 : netCDFRasterBand *poBand =
4021 : new netCDFRasterBand(poDS, var, nDim, lev,
4022 : panBandZLev, panBandDimPos,
4023 180 : paDimIds, i+1 );
4024 :
4025 180 : poDS->SetBand( i+1, poBand );
4026 180 : i++;
4027 : }
4028 :
4029 166 : CPLFree( paDimIds );
4030 166 : CPLFree( panBandDimPos );
4031 166 : CPLFree( panBandZLev );
4032 :
4033 166 : poDS->nBands = i;
4034 :
4035 : // Handle angular geographic coordinates here
4036 :
4037 : /* -------------------------------------------------------------------- */
4038 : /* Initialize any PAM information. */
4039 : /* -------------------------------------------------------------------- */
4040 166 : if( bTreatAsSubdataset )
4041 : {
4042 5 : poDS->SetPhysicalFilename( poDS->osFilename );
4043 5 : poDS->SetSubdatasetName( osSubdatasetName );
4044 : }
4045 :
4046 166 : poDS->TryLoadXML();
4047 :
4048 166 : if( bTreatAsSubdataset )
4049 5 : poDS->oOvManager.Initialize( poDS, ":::VIRTUAL:::" );
4050 : else
4051 161 : poDS->oOvManager.Initialize( poDS, poDS->osFilename );
4052 :
4053 166 : return( poDS );
4054 : }
4055 :
4056 :
4057 : /************************************************************************/
4058 : /* CopyMetadata() */
4059 : /* */
4060 : /* Create a copy of metadata for NC_GLOBAL or a variable */
4061 : /************************************************************************/
4062 :
4063 124 : void CopyMetadata( void *poDS, int fpImage, int CDFVarID ) {
4064 :
4065 : char **papszMetadata;
4066 : char **papszFieldData;
4067 : const char *pszField;
4068 : char szMetaName[ NCDF_MAX_STR_LEN ];
4069 : char szMetaValue[ NCDF_MAX_STR_LEN ];
4070 : char szTemp[ NCDF_MAX_STR_LEN ];
4071 : int nItems;
4072 : int bCopyItem;
4073 :
4074 : /* Remove the following band meta but set them later from band data */
4075 : const char *papszIgnore[] = { CF_ADD_OFFSET, CF_SCALE_FACTOR,
4076 : "valid_range", "_Unsigned",
4077 124 : _FillValue, NULL };
4078 :
4079 124 : if( CDFVarID == NC_GLOBAL ) {
4080 57 : papszMetadata = GDALGetMetadata( (GDALDataset *) poDS,"");
4081 : } else {
4082 67 : papszMetadata = GDALGetMetadata( (GDALRasterBandH) poDS, NULL );
4083 : }
4084 :
4085 124 : nItems = CSLCount( papszMetadata );
4086 :
4087 714 : for(int k=0; k < nItems; k++ ) {
4088 590 : bCopyItem = TRUE;
4089 590 : pszField = CSLGetField( papszMetadata, k );
4090 : papszFieldData = CSLTokenizeString2 (pszField, "=",
4091 590 : CSLT_HONOURSTRINGS );
4092 590 : if( papszFieldData[1] != NULL ) {
4093 590 : strcpy( szMetaName, papszFieldData[ 0 ] );
4094 590 : strcpy( szMetaValue, papszFieldData[ 1 ] );
4095 :
4096 : /* Fix various issues with metadata translation */
4097 590 : if( CDFVarID == NC_GLOBAL ) {
4098 : /* Remove NC_GLOBAL prefix for netcdf global Metadata */
4099 475 : if( strncmp( szMetaName, "NC_GLOBAL#", 10 ) == 0 ) {
4100 100 : strcpy( szTemp, szMetaName+10 );
4101 100 : strcpy( szMetaName, szTemp );
4102 : }
4103 : /* GDAL Metadata renamed as GDAL-[meta] */
4104 375 : else if ( strstr( szMetaName, "#" ) == NULL ) {
4105 38 : strcpy( szTemp, "GDAL_" );
4106 38 : strcat( szTemp, szMetaName );
4107 38 : strcpy( szMetaName, szTemp );
4108 : }
4109 : /* Keep time, lev and depth information for safe-keeping */
4110 : /* Time and vertical coordinate handling need improvements */
4111 337 : else if( strncmp( szMetaName, "time#", 5 ) == 0 ) {
4112 9 : szMetaName[4] = '-';
4113 : }
4114 328 : else if( strncmp( szMetaName, "lev#", 4 ) == 0 ) {
4115 0 : szMetaName[3] = '-';
4116 : }
4117 328 : else if( strncmp( szMetaName, "depth#", 6 ) == 0 ) {
4118 0 : szMetaName[5] = '-';
4119 : }
4120 : /* Only copy data without # (previously all data was copied) */
4121 475 : if ( strstr( szMetaName, "#" ) != NULL ) {
4122 328 : bCopyItem = FALSE;
4123 : }
4124 : /* netCDF attributes do not like the '#' character. */
4125 : // for( unsigned int h=0; h < strlen( szMetaName ) -1 ; h++ ) {
4126 : // if( szMetaName[h] == '#' ) szMetaName[h] = '-';
4127 : // }
4128 : }
4129 : else {
4130 : /* for variables, don't copy varname */
4131 115 : if ( strncmp( szMetaName, "NETCDF_VARNAME", 14) == 0 )
4132 17 : bCopyItem = FALSE;
4133 : /* Don't copy band statistics */
4134 98 : else if ( strncmp( szMetaName, "STATISTICS_", 11) == 0 )
4135 0 : bCopyItem = FALSE;
4136 115 : if ( CSLFindString( (char **)papszIgnore, szMetaName ) != -1 ) {
4137 39 : bCopyItem = FALSE;
4138 : }
4139 : }
4140 :
4141 590 : if ( bCopyItem ) {
4142 206 : if ( NCDFPutAttr( fpImage, CDFVarID,szMetaName,
4143 : szMetaValue ) != CE_None )
4144 : CPLDebug( "GDAL_netCDF", "NCDFPutAttr(%d, %d, %s, %s) failed",
4145 0 : fpImage, CDFVarID,szMetaName, szMetaValue );
4146 : }
4147 :
4148 : }
4149 590 : CSLDestroy( papszFieldData );
4150 : }
4151 :
4152 : /* Set add_offset and scale_factor here if present */
4153 124 : if( CDFVarID != NC_GLOBAL ) {
4154 :
4155 : int bGotAddOffset, bGotScale;
4156 67 : GDALRasterBandH poRB = (GDALRasterBandH) poDS;
4157 67 : double dfAddOffset = GDALGetRasterOffset( poRB , &bGotAddOffset );
4158 67 : double dfScale = GDALGetRasterScale( poRB, &bGotScale );
4159 :
4160 67 : if ( bGotAddOffset && dfAddOffset != 0.0 && bGotScale && dfScale != 1.0 ) {
4161 0 : GDALSetRasterOffset( poRB, dfAddOffset );
4162 0 : GDALSetRasterScale( poRB, dfScale );
4163 : }
4164 :
4165 : }
4166 :
4167 124 : }
4168 :
4169 : /*
4170 : Driver options:
4171 :
4172 : FORMAT=NC/NC2/NC4/NC4C (COMPRESS=DEFLATE sets FORMAT=NC4C)
4173 : COMPRESS=NONE/DEFLATE (default: NONE)
4174 : ZLEVEL=[1-9] (default: 1)
4175 : WRITE_BOTTOMUP=NO/YES (default: NO)
4176 : WRITE_GDAL_TAGS=YES/NO (default: YES)
4177 : WRITE_LONLAT=YES/NO/IF_NEEDED (default: YES for geographic, NO for projected)
4178 : TYPE_LONLAT=float/double (default: double for geographic, float for projected)
4179 : PIXELTYPE=DEFAULT/SIGNEDBYTE (use SIGNEDBYTE to get a signed Byte Band)
4180 :
4181 : Config Options:
4182 :
4183 : GDAL_NETCDF_BOTTOMUP=YES/NO (default:NO) sets the default of WRITE_BOTTOMUP
4184 : and also the default behavior for import
4185 :
4186 : */
4187 :
4188 : /************************************************************************/
4189 : /* CreateLL() */
4190 : /* */
4191 : /* Shared functionality between netCDFDataset::Create() and */
4192 : /* netCDF::CreateCopy() for creating netcdf file based on a set of */
4193 : /* options and a configuration. */
4194 : /************************************************************************/
4195 :
4196 : netCDFDataset *
4197 99 : netCDFDataset::CreateLL( const char * pszFilename,
4198 : int nXSize, int nYSize, int nBands,
4199 : char ** papszOptions )
4200 : {
4201 99 : int status = NC_NOERR;
4202 : netCDFDataset *poDS;
4203 :
4204 99 : poDS = new netCDFDataset();
4205 99 : poDS->nRasterXSize = nXSize;
4206 99 : poDS->nRasterYSize = nYSize;
4207 99 : poDS->eAccess = GA_Update;
4208 198 : poDS->osFilename = pszFilename;
4209 :
4210 : /* from gtiff driver, is this ok? */
4211 : /*
4212 : poDS->nBlockXSize = nXSize;
4213 : poDS->nBlockYSize = 1;
4214 : poDS->nBlocksPerBand =
4215 : ((nYSize + poDS->nBlockYSize - 1) / poDS->nBlockYSize)
4216 : * ((nXSize + poDS->nBlockXSize - 1) / poDS->nBlockXSize);
4217 : */
4218 :
4219 : /* process options */
4220 99 : poDS->papszCreationOptions = CSLDuplicate( papszOptions );
4221 99 : poDS->ProcessCreationOptions( );
4222 :
4223 : /* -------------------------------------------------------------------- */
4224 : /* Create the dataset. */
4225 : /* -------------------------------------------------------------------- */
4226 99 : status = nc_create( pszFilename, poDS->nCreateMode, &(poDS->cdfid) );
4227 :
4228 : /* put into define mode */
4229 99 : poDS->SetDefineMode(TRUE);
4230 :
4231 99 : if( status != NC_NOERR )
4232 : {
4233 : CPLError( CE_Failure, CPLE_OpenFailed,
4234 : "Unable to create netCDF file %s (Error code %d): %s .\n",
4235 2 : pszFilename, status, nc_strerror(status) );
4236 2 : delete poDS;
4237 2 : return NULL;
4238 : }
4239 :
4240 : /* -------------------------------------------------------------------- */
4241 : /* Define dimensions */
4242 : /* -------------------------------------------------------------------- */
4243 97 : poDS->papszDimName.AddString( NCDF_DIMNAME_Y );
4244 97 : poDS->papszDimName.AddString( NCDF_DIMNAME_X );
4245 :
4246 : status = nc_def_dim( poDS->cdfid, NCDF_DIMNAME_X, nXSize,
4247 97 : &(poDS->nXDimID) );
4248 97 : NCDF_ERR(status);
4249 : CPLDebug( "GDAL_netCDF", "status nc_def_dim %s = %d, nXDimID = %d",
4250 97 : NCDF_DIMNAME_X, status, poDS->nXDimID );
4251 :
4252 : status = nc_def_dim( poDS->cdfid, NCDF_DIMNAME_Y, nYSize,
4253 97 : &(poDS->nYDimID) );
4254 97 : NCDF_ERR(status);
4255 : CPLDebug( "GDAL_netCDF", "status nc_def_dim %s = %d, nYDimID = %d",
4256 97 : NCDF_DIMNAME_Y, status, poDS->nYDimID );
4257 :
4258 97 : return poDS;
4259 :
4260 : }
4261 :
4262 : /************************************************************************/
4263 : /* Create() */
4264 : /************************************************************************/
4265 :
4266 : GDALDataset *
4267 40 : netCDFDataset::Create( const char * pszFilename,
4268 : int nXSize, int nYSize, int nBands,
4269 : GDALDataType eType,
4270 : char ** papszOptions )
4271 : {
4272 : netCDFDataset *poDS;
4273 :
4274 : CPLDebug( "GDAL_netCDF",
4275 : "\n=====\nnetCDFDataset::Create( %s, ... )\n",
4276 40 : pszFilename );
4277 :
4278 : poDS = netCDFDataset::CreateLL( pszFilename,
4279 : nXSize, nYSize, nBands,
4280 40 : papszOptions );
4281 :
4282 40 : if ( ! poDS )
4283 0 : return NULL;
4284 :
4285 : /* should we write signed or unsigned byte? */
4286 : /* TODO should this only be done in Create() */
4287 40 : poDS->bSignedData = TRUE;
4288 : const char *pszValue =
4289 40 : CSLFetchNameValue( papszOptions, "PIXELTYPE" );
4290 40 : if( pszValue == NULL )
4291 38 : pszValue = "";
4292 40 : if( eType == GDT_Byte && ( ! EQUAL(pszValue,"SIGNEDBYTE") ) )
4293 10 : poDS->bSignedData = FALSE;
4294 :
4295 : /* -------------------------------------------------------------------- */
4296 : /* Add Conventions, GDAL info and history */
4297 : /* -------------------------------------------------------------------- */
4298 40 : NCDFAddGDALHistory( poDS->cdfid, pszFilename, "", "Create" );
4299 :
4300 : /* -------------------------------------------------------------------- */
4301 : /* Define bands */
4302 : /* -------------------------------------------------------------------- */
4303 218 : for( int iBand = 1; iBand <= nBands; iBand++ )
4304 : {
4305 : poDS->SetBand( iBand, new netCDFRasterBand( poDS, eType, iBand,
4306 69 : poDS->bSignedData ) );
4307 : }
4308 :
4309 : CPLDebug( "GDAL_netCDF",
4310 : "netCDFDataset::Create( %s, ... ) done",
4311 40 : pszFilename );
4312 : /* -------------------------------------------------------------------- */
4313 : /* Return same dataset */
4314 : /* -------------------------------------------------------------------- */
4315 40 : return( poDS );
4316 :
4317 : }
4318 :
4319 :
4320 : template <class T>
4321 67 : CPLErr NCDFCopyBand( GDALRasterBand *poSrcBand, GDALRasterBand *poBand,
4322 : int nXSize, int nYSize,
4323 : GDALProgressFunc pfnProgress, void * pProgressData )
4324 : {
4325 67 : GDALDataType eDT = poSrcBand->GetRasterDataType();
4326 67 : CPLErr eErr = CE_None;
4327 67 : T *patScanline = (T *) CPLMalloc( nXSize * sizeof(T) );
4328 :
4329 3434 : for( int iLine = 0; iLine < nYSize && eErr == CE_None; iLine++ ) {
4330 3367 : eErr = poSrcBand->RasterIO( GF_Read, 0, iLine, nXSize, 1,
4331 : patScanline, nXSize, 1, eDT,
4332 : 0,0);
4333 3367 : if ( eErr == CE_None )
4334 3367 : eErr = poBand->RasterIO( GF_Write, 0, iLine, nXSize, 1,
4335 : patScanline, nXSize, 1, eDT,
4336 : 0,0);
4337 :
4338 3367 : if ( ( nYSize>10 ) && ( iLine % (nYSize/10) == 1 ) ) {
4339 474 : pfnProgress( 1.0*iLine/nYSize , NULL, pProgressData );
4340 : }
4341 : }
4342 :
4343 67 : CPLFree( patScanline );
4344 :
4345 67 : pfnProgress( 1.0, NULL, pProgressData );
4346 :
4347 67 : return eErr;
4348 : }
4349 :
4350 :
4351 : /************************************************************************/
4352 : /* CreateCopy() */
4353 : /************************************************************************/
4354 :
4355 : GDALDataset*
4356 64 : netCDFDataset::CreateCopy( const char * pszFilename, GDALDataset *poSrcDS,
4357 : int bStrict, char ** papszOptions,
4358 : GDALProgressFunc pfnProgress, void * pProgressData )
4359 : {
4360 : netCDFDataset *poDS;
4361 : void *pScaledProgress;
4362 : GDALDataType eDT;
4363 64 : CPLErr eErr = CE_None;
4364 : int nBands, nXSize, nYSize;
4365 : double adfGeoTransform[6];
4366 : const char *pszWKT;
4367 : int iBand;
4368 :
4369 : CPLDebug( "GDAL_netCDF",
4370 : "\n=====\nnetCDFDataset::CreateCopy( %s, ... )\n",
4371 64 : pszFilename );
4372 :
4373 64 : nBands = poSrcDS->GetRasterCount();
4374 64 : nXSize = poSrcDS->GetRasterXSize();
4375 64 : nYSize = poSrcDS->GetRasterYSize();
4376 :
4377 : /* -------------------------------------------------------------------- */
4378 : /* Check input bands for errors */
4379 : /* -------------------------------------------------------------------- */
4380 :
4381 64 : if (nBands == 0)
4382 : {
4383 : CPLError( CE_Failure, CPLE_NotSupported,
4384 1 : "NetCDF driver does not support source dataset with zero band.\n");
4385 1 : return NULL;
4386 : }
4387 :
4388 132 : for( iBand=1; iBand <= nBands; iBand++ )
4389 : {
4390 73 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand( iBand );
4391 73 : eDT = poSrcBand->GetRasterDataType();
4392 73 : if (eDT == GDT_Unknown || GDALDataTypeIsComplex(eDT))
4393 : {
4394 : CPLError( CE_Failure, CPLE_NotSupported,
4395 4 : "NetCDF driver does not support source dataset with band of complex type.");
4396 4 : return NULL;
4397 : }
4398 : }
4399 :
4400 59 : if( !pfnProgress( 0.0, NULL, pProgressData ) )
4401 0 : return NULL;
4402 :
4403 : /* same as in Create() */
4404 : poDS = netCDFDataset::CreateLL( pszFilename,
4405 : nXSize, nYSize, nBands,
4406 59 : papszOptions );
4407 59 : if ( ! poDS )
4408 2 : return NULL;
4409 :
4410 : /* -------------------------------------------------------------------- */
4411 : /* Copy global metadata */
4412 : /* Add Conventions, GDAL info and history */
4413 : /* -------------------------------------------------------------------- */
4414 57 : CopyMetadata((void *) poSrcDS, poDS->cdfid, NC_GLOBAL );
4415 : NCDFAddGDALHistory( poDS->cdfid, pszFilename,
4416 57 : poSrcDS->GetMetadataItem("NC_GLOBAL#history",""),
4417 114 : "CreateCopy" );
4418 :
4419 57 : pfnProgress( 0.1, NULL, pProgressData );
4420 :
4421 : /* -------------------------------------------------------------------- */
4422 : /* Define Bands */
4423 : /* -------------------------------------------------------------------- */
4424 248 : for( iBand=1; iBand <= nBands; iBand++ ) {
4425 : CPLDebug( "GDAL_netCDF", "creating band # %d/%d ",
4426 67 : iBand,nBands );
4427 :
4428 : char szBandName[ NC_MAX_NAME ];
4429 : char szLongName[ NC_MAX_NAME ];
4430 : const char *tmpMetadata;
4431 : netCDFRasterBand *poBand;
4432 67 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand( iBand );
4433 67 : eDT = poSrcBand->GetRasterDataType();
4434 :
4435 : /* Get var name from NETCDF_VARNAME */
4436 67 : tmpMetadata = poSrcBand->GetMetadataItem("NETCDF_VARNAME");
4437 67 : if( tmpMetadata != NULL) {
4438 17 : if( nBands > 1 ) sprintf(szBandName,"%s%d",tmpMetadata,iBand);
4439 17 : else strcpy( szBandName, tmpMetadata );
4440 : // poSrcBand->SetMetadataItem("NETCDF_VARNAME","");
4441 : }
4442 : else
4443 50 : szBandName[0]='\0';
4444 :
4445 : /* Get long_name from <var>#long_name */
4446 : sprintf(szLongName,"%s#%s",
4447 67 : poSrcBand->GetMetadataItem("NETCDF_VARNAME"),
4448 134 : CF_LNG_NAME);
4449 67 : tmpMetadata = poSrcDS->GetMetadataItem(szLongName);
4450 67 : if( tmpMetadata != NULL)
4451 16 : strcpy( szLongName, tmpMetadata);
4452 : else
4453 51 : szLongName[0]='\0';
4454 :
4455 67 : int bSignedData = TRUE;
4456 67 : if ( eDT == GDT_Byte ) {
4457 : /* GDAL defaults to unsigned bytes, but check if metadata says its
4458 : signed, as NetCDF can support this for certain formats. */
4459 52 : bSignedData = FALSE;
4460 : tmpMetadata = poSrcBand->GetMetadataItem("PIXELTYPE",
4461 52 : "IMAGE_STRUCTURE");
4462 52 : if ( tmpMetadata && EQUAL(tmpMetadata,"SIGNEDBYTE") )
4463 1 : bSignedData = TRUE;
4464 : }
4465 :
4466 : poBand = new netCDFRasterBand( poDS, eDT, iBand,
4467 : bSignedData,
4468 67 : szBandName, szLongName );
4469 67 : poDS->SetBand( iBand, poBand );
4470 :
4471 : /* Copy Metadata for band */
4472 67 : poBand->SetNoDataValue( poSrcBand->GetNoDataValue(0) );
4473 : CopyMetadata( (void *) GDALGetRasterBand( poSrcDS, iBand ),
4474 67 : poDS->cdfid, poBand->nZId );
4475 :
4476 : }
4477 :
4478 57 : pfnProgress( 0.2, NULL, pProgressData );
4479 :
4480 : /* -------------------------------------------------------------------- */
4481 : /* Copy GeoTransform and Projection */
4482 : /* -------------------------------------------------------------------- */
4483 57 : eErr = poSrcDS->GetGeoTransform( adfGeoTransform );
4484 57 : if ( eErr == CE_None ) {
4485 57 : poDS->SetGeoTransform( adfGeoTransform );
4486 : // set this is to disable AddProjectionVars() from being called
4487 57 : poDS->bSetGeoTransform = FALSE;
4488 : }
4489 :
4490 57 : pszWKT = poSrcDS->GetProjectionRef( );
4491 57 : if ( pszWKT ) {
4492 57 : poDS->SetProjection( pszWKT );
4493 : // now we can call AddProjectionVars()
4494 57 : poDS->bSetGeoTransform = TRUE;
4495 : pScaledProgress = GDALCreateScaledProgress( 0.20, 0.50, pfnProgress,
4496 57 : pProgressData );
4497 57 : poDS->AddProjectionVars( GDALScaledProgress, pScaledProgress );
4498 57 : GDALDestroyScaledProgress( pScaledProgress );
4499 :
4500 : }
4501 :
4502 57 : pfnProgress( 0.5, NULL, pProgressData );
4503 :
4504 : /* -------------------------------------------------------------------- */
4505 : /* Write Bands */
4506 : /* -------------------------------------------------------------------- */
4507 : /* make sure we are in data mode */
4508 57 : poDS->SetDefineMode( FALSE );
4509 :
4510 : double dfTemp,dfTemp2;
4511 57 : dfTemp = dfTemp2 = 0.5;
4512 :
4513 124 : for( iBand=1; iBand <= nBands; iBand++ ) {
4514 :
4515 67 : dfTemp2 = dfTemp + 0.4/nBands;
4516 : pScaledProgress =
4517 : GDALCreateScaledProgress( dfTemp, dfTemp2,
4518 67 : pfnProgress, pProgressData );
4519 67 : dfTemp = dfTemp2;
4520 :
4521 67 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand( iBand );
4522 :
4523 : /* new */
4524 : GDALDataType eDT;
4525 67 : CPLErr eErr = CE_None;
4526 :
4527 : CPLDebug( "GDAL_netCDF", "copying band data # %d/%d ",
4528 67 : iBand,nBands );
4529 :
4530 : // hBand = GDALGetRasterBand( poSrcDS, i );
4531 67 : GDALRasterBand *poBand = poDS->GetRasterBand( iBand );
4532 67 : eDT = poSrcBand->GetRasterDataType();
4533 67 : eErr = CE_None;
4534 :
4535 : /* -------------------------------------------------------------------- */
4536 : /* Copy Band data */
4537 : /* -------------------------------------------------------------------- */
4538 67 : if( eDT == GDT_Byte ) {
4539 52 : CPLDebug( "GDAL_netCDF", "GByte Band#%d", iBand );
4540 : NCDFCopyBand<GByte>( poSrcBand, poBand, nXSize, nYSize,
4541 52 : GDALScaledProgress, pScaledProgress );
4542 : }
4543 18 : else if( ( eDT == GDT_UInt16 ) || ( eDT == GDT_Int16 ) ) {
4544 3 : CPLDebug( "GDAL_netCDF", "GInt16 Band#%d", iBand );
4545 : NCDFCopyBand<GInt16>( poSrcBand, poBand, nXSize, nYSize,
4546 3 : GDALScaledProgress, pScaledProgress );
4547 : }
4548 15 : else if( (eDT == GDT_UInt32) || (eDT == GDT_Int32) ) {
4549 3 : CPLDebug( "GDAL_netCDF", "GInt16 Band#%d", iBand );
4550 : NCDFCopyBand<GInt32>( poSrcBand, poBand, nXSize, nYSize,
4551 3 : GDALScaledProgress, pScaledProgress );
4552 : }
4553 9 : else if( (eDT == GDT_Float32) ) {
4554 7 : CPLDebug( "GDAL_netCDF", "float Band#%d", iBand);
4555 : NCDFCopyBand<float>( poSrcBand, poBand, nXSize, nYSize,
4556 7 : GDALScaledProgress, pScaledProgress );
4557 : }
4558 2 : else if( (eDT == GDT_Float64) ) {
4559 2 : CPLDebug( "GDAL_netCDF", "double Band#%d", iBand);
4560 : NCDFCopyBand<double>( poSrcBand, poBand, nXSize, nYSize,
4561 2 : GDALScaledProgress, pScaledProgress );
4562 : }
4563 : else {
4564 : CPLError( CE_Failure, CPLE_NotSupported,
4565 : "The NetCDF driver does not support GDAL data type %d",
4566 0 : eDT );
4567 : }
4568 :
4569 67 : GDALDestroyScaledProgress( pScaledProgress );
4570 :
4571 : }
4572 :
4573 : /* -------------------------------------------------------------------- */
4574 : /* Cleanup and close. */
4575 : /* -------------------------------------------------------------------- */
4576 57 : delete( poDS );
4577 : // CPLFree(pszProj4Defn );
4578 :
4579 57 : pfnProgress( 0.95, NULL, pProgressData );
4580 :
4581 : /* -------------------------------------------------------------------- */
4582 : /* Re-open dataset, and copy any auxilary pam information. */
4583 : /* Disable PAM, at least temporarily. See bug #4244 */
4584 : /* -------------------------------------------------------------------- */
4585 57 : poDS = (netCDFDataset *) GDALOpen( pszFilename, GA_ReadOnly );
4586 :
4587 : // if( poDS )
4588 : // poDS->CloneInfo( poSrcDS, GCIF_PAM_DEFAULT );
4589 :
4590 57 : pfnProgress( 1.0, NULL, pProgressData );
4591 :
4592 57 : return poDS;
4593 : }
4594 :
4595 :
4596 : /* note: some logic depends on bIsProjected and bIsGeoGraphic */
4597 : /* which may not be known when Create() is called */
4598 : void
4599 99 : netCDFDataset::ProcessCreationOptions( )
4600 : {
4601 : const char *pszValue;
4602 :
4603 : /* File format */
4604 99 : nFormat = NCDF_FORMAT_NC;
4605 99 : pszValue = CSLFetchNameValue( papszCreationOptions, "FORMAT" );
4606 99 : if ( pszValue != NULL ) {
4607 0 : if ( EQUAL( pszValue, "NC" ) ) {
4608 0 : nFormat = NCDF_FORMAT_NC;
4609 : }
4610 : #ifdef NETCDF_HAS_NC2
4611 0 : else if ( EQUAL( pszValue, "NC2" ) ) {
4612 0 : nFormat = NCDF_FORMAT_NC2;
4613 : }
4614 : #endif
4615 : #ifdef NETCDF_HAS_NC4
4616 : else if ( EQUAL( pszValue, "NC4" ) ) {
4617 : nFormat = NCDF_FORMAT_NC4;
4618 : }
4619 : else if ( EQUAL( pszValue, "NC4C" ) ) {
4620 : nFormat = NCDF_FORMAT_NC4C;
4621 : }
4622 : #endif
4623 : else {
4624 : CPLError( CE_Failure, CPLE_NotSupported,
4625 0 : "FORMAT=%s in not supported, using the default NC format.", pszValue );
4626 : }
4627 : }
4628 :
4629 : /* compression only available for NC4 */
4630 : #ifdef NETCDF_HAS_NC4
4631 :
4632 : /* COMPRESS option */
4633 : pszValue = CSLFetchNameValue( papszCreationOptions, "COMPRESS" );
4634 : if ( pszValue != NULL ) {
4635 : if ( EQUAL( pszValue, "NONE" ) ) {
4636 : nCompress = NCDF_COMPRESS_NONE;
4637 : }
4638 : else if ( EQUAL( pszValue, "DEFLATE" ) ) {
4639 : nCompress = NCDF_COMPRESS_DEFLATE;
4640 : if ( !((nFormat == NCDF_FORMAT_NC4) || (nFormat == NCDF_FORMAT_NC4C)) ) {
4641 : CPLError( CE_Warning, CPLE_IllegalArg,
4642 : "NOTICE: Format set to NC4C because compression is set to DEFLATE." );
4643 : nFormat = NCDF_FORMAT_NC4C;
4644 : }
4645 : }
4646 : else {
4647 : CPLError( CE_Failure, CPLE_NotSupported,
4648 : "COMPRESS=%s is not supported.", pszValue );
4649 : }
4650 : }
4651 :
4652 : /* ZLEVEL option */
4653 : pszValue = CSLFetchNameValue( papszCreationOptions, "ZLEVEL" );
4654 : if( pszValue != NULL )
4655 : {
4656 : nZLevel = atoi( pszValue );
4657 : if (!(nZLevel >= 1 && nZLevel <= 9))
4658 : {
4659 : CPLError( CE_Warning, CPLE_IllegalArg,
4660 : "ZLEVEL=%s value not recognised, ignoring.",
4661 : pszValue );
4662 : nZLevel = NCDF_DEFLATE_LEVEL;
4663 : }
4664 : }
4665 :
4666 : #endif
4667 :
4668 : /* set nCreateMode based on nFormat */
4669 99 : switch ( nFormat ) {
4670 : #ifdef NETCDF_HAS_NC2
4671 : case NCDF_FORMAT_NC2:
4672 0 : nCreateMode = NC_CLOBBER|NC_64BIT_OFFSET;
4673 0 : break;
4674 : #endif
4675 : #ifdef NETCDF_HAS_NC4
4676 : case NCDF_FORMAT_NC4:
4677 : nCreateMode = NC_CLOBBER|NC_NETCDF4;
4678 : break;
4679 : case NCDF_FORMAT_NC4C:
4680 : nCreateMode = NC_CLOBBER|NC_NETCDF4|NC_CLASSIC_MODEL;
4681 : break;
4682 : #endif
4683 : case NCDF_FORMAT_NC:
4684 : default:
4685 99 : nCreateMode = NC_CLOBBER;
4686 : break;
4687 : }
4688 :
4689 : CPLDebug( "GDAL_netCDF",
4690 : "file options: format=%d compress=%d zlevel=%d",
4691 99 : nFormat, nCompress, nZLevel );
4692 :
4693 : /* netcdf standard is bottom-up */
4694 : /* overriden by config option GDAL_NETCDF_BOTTOMUP and -co option WRITE_BOTTOMUP */
4695 : // bBottomUp = CSLTestBoolean( CPLGetConfigOption( "GDAL_NETCDF_BOTTOMUP", "YES" ) );
4696 99 : bBottomUp = CSLTestBoolean( CPLGetConfigOption( "GDAL_NETCDF_BOTTOMUP", "NO" ) );
4697 99 : bBottomUp = CSLFetchBoolean( papszCreationOptions, "WRITE_BOTTOMUP", bBottomUp );
4698 :
4699 : /* TODO could add a config option GDAL_NETCDF_PREF=GDAL/CF */
4700 :
4701 : /* moved over projection stuff to setprojection */
4702 :
4703 99 : }
4704 :
4705 136 : int netCDFDataset::DefVarDeflate( int nVarId, int bChunking )
4706 : {
4707 : #ifdef NETCDF_HAS_NC4
4708 : // must set chunk size to avoid huge performace hit
4709 : // perhaps another solution it to change the chunk cache?
4710 : // http://www.unidata.ucar.edu/software/netcdf/docs/netcdf.html#Chunk-Cache
4711 : if ( nCompress == NCDF_COMPRESS_DEFLATE ) {
4712 : status = nc_def_var_deflate(cdfid,nVarId,1,1,nZLevel);
4713 : NCDF_ERR(status);
4714 : if ( status != NC_NOERR && bChunking ) {
4715 : //TODO make sure this is ok
4716 : size_t chunksize[] = { 1, nRasterXSize };
4717 : status = nc_def_var_chunking( cdfid, nVarId,
4718 : NC_CHUNKED, chunksize );
4719 : NCDF_ERR(status);
4720 : }
4721 : return status;
4722 : }
4723 : #endif
4724 136 : return NC_NOERR;
4725 : }
4726 :
4727 :
4728 : /************************************************************************/
4729 : /* GDALRegister_netCDF() */
4730 : /************************************************************************/
4731 :
4732 558 : void GDALRegister_netCDF()
4733 :
4734 : {
4735 558 : if (! GDAL_CHECK_VERSION("netCDF driver"))
4736 0 : return;
4737 :
4738 558 : if( GDALGetDriverByName( "netCDF" ) == NULL )
4739 : {
4740 : GDALDriver *poDriver;
4741 : char szCreateOptions[3072];
4742 :
4743 537 : poDriver = new GDALDriver( );
4744 :
4745 : /* -------------------------------------------------------------------- */
4746 : /* Build full creation option list. */
4747 : /* -------------------------------------------------------------------- */
4748 : sprintf( szCreateOptions, "%s",
4749 : "<CreationOptionList>"
4750 : " <Option name='FORMAT' type='string-select' default='NC'>"
4751 : " <Value>NC</Value>"
4752 : #ifdef NETCDF_HAS_NC2
4753 : " <Value>NC2</Value>"
4754 : #endif
4755 : #ifdef NETCDF_HAS_NC4
4756 : " <Value>NC4</Value>"
4757 : " <Value>NC4C</Value>"
4758 : #endif
4759 : " </Option>"
4760 : #ifdef NETCDF_HAS_NC4
4761 : " <Option name='COMPRESS' type='string-select' default='NONE'>"
4762 : " <Value>NONE</Value>"
4763 : " <Value>DEFLATE</Value>"
4764 : " </Option>"
4765 : " <Option name='ZLEVEL' type='int' description='DEFLATE compression level 1-9' default='1'/>"
4766 : #endif
4767 : " <Option name='WRITE_BOTTOMUP' type='boolean' default='NO'>"
4768 : " </Option>"
4769 : " <Option name='WRITE_GDAL_TAGS' type='boolean' default='YES'>"
4770 : " </Option>"
4771 : " <Option name='WRITE_LONLAT' type='string-select'>"
4772 : " <Value>YES</Value>"
4773 : " <Value>NO</Value>"
4774 : " <Value>IF_NEEDED</Value>"
4775 : " </Option>"
4776 : " <Option name='TYPE_LONLAT' type='string-select'>"
4777 : " <Value>float</Value>"
4778 : " <Value>double</Value>"
4779 : " </Option>"
4780 : " <Option name='PIXELTYPE' type='string-select' description='only used in Create()'>"
4781 : " <Value>DEFAULT</Value>"
4782 : " <Value>SIGNEDBYTE</Value>"
4783 : " </Option>"
4784 537 : "</CreationOptionList>" );
4785 :
4786 :
4787 : /* -------------------------------------------------------------------- */
4788 : /* Set the driver details. */
4789 : /* -------------------------------------------------------------------- */
4790 537 : poDriver->SetDescription( "netCDF" );
4791 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
4792 537 : "Network Common Data Format" );
4793 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
4794 537 : "frmt_netcdf.html" );
4795 537 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "nc" );
4796 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONOPTIONLIST,
4797 537 : szCreateOptions );
4798 :
4799 : /* make driver config and capabilities available */
4800 537 : poDriver->SetMetadataItem( "NETCDF_VERSION", nc_inq_libvers() );
4801 537 : poDriver->SetMetadataItem( "NETCDF_CONVENTIONS", NCDF_CONVENTIONS_CF );
4802 : #ifdef NETCDF_HAS_NC2
4803 537 : poDriver->SetMetadataItem( "NETCDF_HAS_NC2", "YES" );
4804 : #endif
4805 : #ifdef NETCDF_HAS_NC4
4806 : poDriver->SetMetadataItem( "NETCDF_HAS_NC4", "YES" );
4807 : #endif
4808 : #ifdef NETCDF_HAS_HDF4
4809 : poDriver->SetMetadataItem( "NETCDF_HAS_HDF4", "YES" );
4810 : #endif
4811 : #ifdef HAVE_HDF4
4812 537 : poDriver->SetMetadataItem( "GDAL_HAS_HDF4", "YES" );
4813 : #endif
4814 : #ifdef HAVE_HDF5
4815 537 : poDriver->SetMetadataItem( "GDAL_HAS_HDF5", "YES" );
4816 : #endif
4817 :
4818 : /* set pfns and register driver */
4819 537 : poDriver->pfnOpen = netCDFDataset::Open;
4820 537 : poDriver->pfnCreateCopy = netCDFDataset::CreateCopy;
4821 537 : poDriver->pfnCreate = netCDFDataset::Create;
4822 537 : poDriver->pfnIdentify = netCDFDataset::Identify;
4823 :
4824 537 : GetGDALDriverManager( )->RegisterDriver( poDriver );
4825 : }
4826 : }
4827 :
4828 : /************************************************************************/
4829 : /* New functions */
4830 : /************************************************************************/
4831 :
4832 : /* Test for GDAL version string >= target */
4833 147 : int NCDFIsGDALVersionGTE(const char* pszVersion, int nTarget)
4834 : {
4835 147 : int nVersion = 0;
4836 147 : int nVersions [] = {0,0,0,0};
4837 : char **papszTokens;
4838 :
4839 : /* Valid strings are "GDAL 1.9dev, released 2011/01/18" and "GDAL 1.8.1 " */
4840 147 : if ( pszVersion == NULL || EQUAL( pszVersion, "" ) )
4841 0 : return FALSE;
4842 147 : else if ( ! EQUALN("GDAL ", pszVersion, 5) )
4843 0 : return FALSE;
4844 147 : else if ( EQUALN("GDAL 1.9dev", pszVersion,11 ) )
4845 147 : return nTarget <= 1900;
4846 0 : else if ( EQUALN("GDAL 1.8dev", pszVersion,11 ) )
4847 0 : return nTarget <= 1800;
4848 :
4849 0 : papszTokens = CSLTokenizeString2( pszVersion+5, ".", 0 );
4850 :
4851 0 : for ( int iToken = 0; papszTokens && papszTokens[iToken]; iToken++ ) {
4852 0 : nVersions[iToken] = atoi( papszTokens[iToken] );
4853 : }
4854 : /* (GDAL_VERSION_MAJOR*1000+GDAL_VERSION_MINOR*100+GDAL_VERSION_REV*10+GDAL_VERSION_BUILD) */
4855 0 : nVersion = nVersions[0]*1000 + nVersions[1]*100 +
4856 0 : nVersions[2]*10 + nVersions[3];
4857 :
4858 0 : CSLDestroy( papszTokens );
4859 0 : return nTarget <= nVersion;
4860 : }
4861 :
4862 : /* Add Conventions, GDAL version and history */
4863 97 : void NCDFAddGDALHistory( int fpImage,
4864 : const char * pszFilename, const char *pszOldHist,
4865 : const char * pszFunctionName)
4866 : {
4867 : char szTemp[NC_MAX_NAME];
4868 :
4869 : nc_put_att_text( fpImage, NC_GLOBAL, "Conventions",
4870 : strlen(NCDF_CONVENTIONS_CF),
4871 97 : NCDF_CONVENTIONS_CF );
4872 :
4873 : nc_put_att_text( fpImage, NC_GLOBAL, "GDAL",
4874 97 : strlen(NCDF_GDAL), NCDF_GDAL );
4875 :
4876 : /* Add history */
4877 : #ifdef GDAL_SET_CMD_LINE_DEFINED_TMP
4878 : if ( ! EQUAL(GDALGetCmdLine(), "" ) )
4879 : strcpy( szTemp, GDALGetCmdLine() );
4880 : else
4881 : sprintf( szTemp, "GDAL %s( %s, ... )",pszFunctionName,pszFilename );
4882 : #else
4883 97 : sprintf( szTemp, "GDAL %s( %s, ... )",pszFunctionName,pszFilename );
4884 : #endif
4885 :
4886 97 : NCDFAddHistory( fpImage, szTemp, pszOldHist );
4887 :
4888 97 : }
4889 :
4890 : /* code taken from cdo and libcdi, used for writing the history attribute */
4891 : //void cdoDefHistory(int fileID, char *histstring)
4892 97 : void NCDFAddHistory(int fpImage, const char *pszAddHist, const char *pszOldHist)
4893 : {
4894 : char strtime[32];
4895 : time_t tp;
4896 : struct tm *ltime;
4897 :
4898 97 : char *pszNewHist = NULL;
4899 97 : size_t nNewHistSize = 0;
4900 97 : int disableHistory = FALSE;
4901 : int status;
4902 :
4903 : /* Check pszOldHist - as if there was no previous history, it will be
4904 : a null pointer - if so set as empty. */
4905 97 : if (NULL == pszOldHist) {
4906 41 : pszOldHist = "";
4907 : }
4908 :
4909 97 : tp = time(NULL);
4910 97 : if ( tp != -1 )
4911 : {
4912 97 : ltime = localtime(&tp);
4913 97 : (void) strftime(strtime, sizeof(strtime), "%a %b %d %H:%M:%S %Y: ", ltime);
4914 : }
4915 :
4916 : // status = nc_get_att_text( fpImage, NC_GLOBAL,
4917 : // "history", pszOldHist );
4918 : // printf("status: %d pszOldHist: [%s]\n",status,pszOldHist);
4919 :
4920 97 : nNewHistSize = strlen(pszOldHist)+strlen(strtime)+strlen(pszAddHist)+1+1;
4921 97 : pszNewHist = (char *) CPLMalloc(nNewHistSize * sizeof(char));
4922 :
4923 97 : strcpy(pszNewHist, strtime);
4924 97 : strcat(pszNewHist, pszAddHist);
4925 :
4926 97 : if ( disableHistory == FALSE && pszNewHist )
4927 : {
4928 97 : if ( ! EQUAL(pszOldHist,"") )
4929 16 : strcat(pszNewHist, "\n");
4930 97 : strcat(pszNewHist, pszOldHist);
4931 : }
4932 :
4933 : status = nc_put_att_text( fpImage, NC_GLOBAL,
4934 : "history", strlen(pszNewHist),
4935 97 : pszNewHist );
4936 97 : NCDF_ERR(status);
4937 :
4938 97 : CPLFree(pszNewHist);
4939 97 : }
4940 :
4941 :
4942 41 : int NCDFIsCfProjection( const char* pszProjection )
4943 : {
4944 : /* Find the appropriate mapping */
4945 823 : for (int iMap = 0; poNetcdfSRS_PT[iMap].WKT_SRS != NULL; iMap++ ) {
4946 : // printf("now at %d, proj=%s\n",i, poNetcdfSRS_PT[i].GDAL_SRS);
4947 823 : if ( EQUAL( pszProjection, poNetcdfSRS_PT[iMap].WKT_SRS ) ) {
4948 41 : if ( poNetcdfSRS_PT[iMap].mappings != NULL )
4949 41 : return TRUE;
4950 : else
4951 0 : return FALSE;
4952 : }
4953 : }
4954 0 : return FALSE;
4955 : }
4956 :
4957 :
4958 : /* Write any needed projection attributes *
4959 : * poPROJCS: ptr to proj crd system
4960 : * pszProjection: name of projection system in GDAL WKT
4961 : * fpImage: open NetCDF file in writing mode
4962 : * NCDFVarID: NetCDF Var Id of proj system we're writing in to
4963 : *
4964 : * The function first looks for the oNetcdfSRS_PP mapping object
4965 : * that corresponds to the input projection name. If none is found
4966 : * the generic mapping is used. In the case of specific mappings,
4967 : * the driver looks for each attribute listed in the mapping object
4968 : * and then looks up the value within the OGR_SRSNode. In the case
4969 : * of the generic mapping, the lookup is reversed (projection params,
4970 : * then mapping). For more generic code, GDAL->NETCDF
4971 : * mappings and the associated value are saved in std::map objects.
4972 : */
4973 :
4974 : /* NOTE modifications by ET to combine the specific and generic mappings */
4975 :
4976 41 : void NCDFWriteProjAttribs( const OGR_SRSNode *poPROJCS,
4977 : const char* pszProjection,
4978 : const int fpImage, const int NCDFVarID )
4979 : {
4980 : double dfStdP[2];
4981 41 : int bFoundStdP1=FALSE,bFoundStdP2=FALSE;
4982 41 : double dfValue=0.0;
4983 : const char *pszParamStr, *pszParamVal;
4984 : const std::string *pszNCDFAtt, *pszGDALAtt;
4985 : static const oNetcdfSRS_PP *poMap = NULL;
4986 41 : int nMapIndex = -1;
4987 41 : int bWriteVal = FALSE;
4988 :
4989 : //Attribute <GDAL,NCDF> and Value <NCDF,value> mappings
4990 41 : std::map< std::string, std::string > oAttMap;
4991 41 : std::map< std::string, std::string >::iterator oAttIter;
4992 41 : std::map< std::string, double > oValMap;
4993 41 : std::map< std::string, double >::iterator oValIter, oValIter2;
4994 : //results to write
4995 41 : std::vector< std::pair<std::string,double> > oOutList;
4996 :
4997 : /* Find the appropriate mapping */
4998 823 : for (int iMap = 0; poNetcdfSRS_PT[iMap].WKT_SRS != NULL; iMap++ ) {
4999 823 : if ( EQUAL( pszProjection, poNetcdfSRS_PT[iMap].WKT_SRS ) ) {
5000 41 : nMapIndex = iMap;
5001 41 : poMap = poNetcdfSRS_PT[iMap].mappings;
5002 41 : break;
5003 : }
5004 : }
5005 :
5006 : //ET TODO if projection name is not found, should we do something special?
5007 41 : if ( nMapIndex == -1 ) {
5008 : CPLError( CE_Warning, CPLE_AppDefined,
5009 : "projection name %s not found in the lookup tables!!!",
5010 0 : pszProjection);
5011 : }
5012 : /* if no mapping was found or assigned, set the generic one */
5013 41 : if ( !poMap ) {
5014 : CPLError( CE_Warning, CPLE_AppDefined,
5015 : "projection name %s in not part of the CF standard, will not be supported by CF!",
5016 0 : pszProjection);
5017 0 : poMap = poGenericMappings;
5018 : }
5019 :
5020 : /* initialize local map objects */
5021 470 : for ( int iMap = 0; poMap[iMap].WKT_ATT != NULL; iMap++ ) {
5022 194 : oAttMap[poMap[iMap].WKT_ATT] = poMap[iMap].CF_ATT;
5023 : }
5024 :
5025 415 : for( int iChild = 0; iChild < poPROJCS->GetChildCount(); iChild++ ) {
5026 :
5027 : const OGR_SRSNode *poNode;
5028 :
5029 374 : poNode = poPROJCS->GetChild( iChild );
5030 374 : if( !EQUAL(poNode->GetValue(),"PARAMETER")
5031 : || poNode->GetChildCount() != 2 )
5032 180 : continue;
5033 194 : pszParamStr = poNode->GetChild(0)->GetValue();
5034 194 : pszParamVal = poNode->GetChild(1)->GetValue();
5035 :
5036 194 : oValMap[pszParamStr] = atof(pszParamVal);
5037 : }
5038 :
5039 : /* Lookup mappings and fill output vector */
5040 41 : if ( poMap != poGenericMappings ) { /* specific mapping, loop over mapping values */
5041 :
5042 235 : for ( oAttIter = oAttMap.begin(); oAttIter != oAttMap.end(); oAttIter++ ) {
5043 :
5044 194 : pszGDALAtt = &(oAttIter->first);
5045 194 : pszNCDFAtt = &(oAttIter->second);
5046 194 : oValIter = oValMap.find( *pszGDALAtt );
5047 :
5048 194 : if ( oValIter != oValMap.end() ) {
5049 :
5050 194 : dfValue = oValIter->second;
5051 194 : bWriteVal = TRUE;
5052 :
5053 : /* special case for PS (Polar Stereographic) grid
5054 : See comments in netcdfdataset.h for this projection. */
5055 194 : if ( EQUAL( SRS_PP_LATITUDE_OF_ORIGIN, pszGDALAtt->c_str() ) &&
5056 : EQUAL(pszProjection, SRS_PT_POLAR_STEREOGRAPHIC) ) {
5057 3 : double dfLatPole = 0.0;
5058 3 : if ( dfValue > 0.0) dfLatPole = 90.0;
5059 3 : else dfLatPole = -90.0;
5060 : oOutList.push_back( std::make_pair( CF_PP_LAT_PROJ_ORIGIN,
5061 3 : dfLatPole ) );
5062 : }
5063 : /* special case for LCC-1SP
5064 : See comments in netcdfdataset.h for this projection. */
5065 191 : else if ( EQUAL( SRS_PP_SCALE_FACTOR, pszGDALAtt->c_str() ) &&
5066 : EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_1SP) ) {
5067 : /* default is to not write as it is not CF-1 */
5068 0 : bWriteVal = FALSE;
5069 : /* test if there is no standard_parallel1 */
5070 0 : if ( oValMap.find( std::string(CF_PP_STD_PARALLEL_1) ) == oValMap.end() ) {
5071 : /* if scale factor != 1.0 write value for GDAL, but this is not supported by CF-1 */
5072 0 : if ( !CPLIsEqual(dfValue,1.0) ) {
5073 : CPLError( CE_Failure, CPLE_NotSupported,
5074 : "NetCDF driver export of LCC-1SP with scale factor != 1.0 "
5075 : "and no standard_parallel1 is not CF-1 (bug #3324).\n"
5076 0 : "Use the 2SP variant which is supported by CF." );
5077 0 : bWriteVal = TRUE;
5078 : }
5079 : /* else copy standard_parallel1 from latitude_of_origin, because scale_factor=1.0 */
5080 : else {
5081 0 : oValIter2 = oValMap.find( std::string(SRS_PP_LATITUDE_OF_ORIGIN) );
5082 0 : if (oValIter2 != oValMap.end() ) {
5083 : oOutList.push_back( std::make_pair( CF_PP_STD_PARALLEL_1,
5084 0 : oValIter2->second) );
5085 : }
5086 : else {
5087 : CPLError( CE_Failure, CPLE_NotSupported,
5088 : "NetCDF driver export of LCC-1SP with no standard_parallel1 "
5089 0 : "and no latitude_of_origin is not suuported (bug #3324).");
5090 : }
5091 : }
5092 : }
5093 : }
5094 194 : if ( bWriteVal )
5095 194 : oOutList.push_back( std::make_pair( *pszNCDFAtt, dfValue ) );
5096 :
5097 : }
5098 : // else printf("NOT FOUND!!!\n");
5099 : }
5100 :
5101 : }
5102 : else { /* generic mapping, loop over projected values */
5103 :
5104 0 : for ( oValIter = oValMap.begin(); oValIter != oValMap.end(); oValIter++ ) {
5105 :
5106 0 : pszGDALAtt = &(oValIter->first);
5107 0 : dfValue = oValIter->second;
5108 :
5109 0 : oAttIter = oAttMap.find( *pszGDALAtt );
5110 :
5111 0 : if ( oAttIter != oAttMap.end() ) {
5112 0 : oOutList.push_back( std::make_pair( oAttIter->second, dfValue ) );
5113 : }
5114 : /* for SRS_PP_SCALE_FACTOR write 2 mappings */
5115 0 : else if ( EQUAL(pszGDALAtt->c_str(), SRS_PP_SCALE_FACTOR) ) {
5116 : oOutList.push_back( std::make_pair( CF_PP_SCALE_FACTOR_MERIDIAN,
5117 0 : dfValue ) );
5118 : oOutList.push_back( std::make_pair( CF_PP_SCALE_FACTOR_ORIGIN,
5119 0 : dfValue ) );
5120 : }
5121 : /* if not found insert the GDAL name */
5122 : else {
5123 0 : oOutList.push_back( std::make_pair( *pszGDALAtt, dfValue ) );
5124 : }
5125 : }
5126 : }
5127 :
5128 : /* Write all the values that were found */
5129 : // std::vector< std::pair<std::string,double> >::reverse_iterator it;
5130 : // for (it = oOutList.rbegin(); it != oOutList.rend(); it++ ) {
5131 41 : std::vector< std::pair<std::string,double> >::iterator it;
5132 238 : for (it = oOutList.begin(); it != oOutList.end(); it++ ) {
5133 197 : pszParamVal = (it->first).c_str();
5134 197 : dfValue = it->second;
5135 : /* Handle the STD_PARALLEL attrib */
5136 197 : if( EQUAL( pszParamVal, CF_PP_STD_PARALLEL_1 ) ) {
5137 16 : bFoundStdP1 = TRUE;
5138 16 : dfStdP[0] = dfValue;
5139 : }
5140 181 : else if( EQUAL( pszParamVal, CF_PP_STD_PARALLEL_2 ) ) {
5141 7 : bFoundStdP2 = TRUE;
5142 7 : dfStdP[1] = dfValue;
5143 : }
5144 : else {
5145 : nc_put_att_double( fpImage, NCDFVarID, pszParamVal,
5146 174 : NC_DOUBLE, 1,&dfValue );
5147 : }
5148 : }
5149 : /* Now write the STD_PARALLEL attrib */
5150 41 : if ( bFoundStdP1 ) {
5151 : /* one value or equal values */
5152 25 : if ( !bFoundStdP2 || dfStdP[0] == dfStdP[1] ) {
5153 : nc_put_att_double( fpImage, NCDFVarID, CF_PP_STD_PARALLEL,
5154 9 : NC_DOUBLE, 1, &dfStdP[0] );
5155 : }
5156 : else { /* two values */
5157 : nc_put_att_double( fpImage, NCDFVarID, CF_PP_STD_PARALLEL,
5158 7 : NC_DOUBLE, 2, dfStdP );
5159 : }
5160 41 : }
5161 41 : }
5162 :
5163 2483 : CPLErr NCDFSafeStrcat(char** ppszDest, char* pszSrc, size_t* nDestSize)
5164 : {
5165 : /* Reallocate the data string until the content fits */
5166 7092 : while(*nDestSize < (strlen(*ppszDest) + strlen(pszSrc) + 1)) {
5167 2126 : (*nDestSize) *= 2;
5168 2126 : *ppszDest = (char*) CPLRealloc((void*) *ppszDest, *nDestSize);
5169 : }
5170 2483 : strcat(*ppszDest, pszSrc);
5171 :
5172 2483 : return CE_None;
5173 : }
5174 :
5175 : /* helper function for NCDFGetAttr() */
5176 : /* sets pdfValue to first value returned */
5177 : /* and if bSetPszValue=True sets szValue with all values attribute values */
5178 : /* pszValue is the responsibility of the caller and must be freed */
5179 4651 : CPLErr NCDFGetAttr1( int nCdfId, int nVarId, const char *pszAttrName,
5180 : double *pdfValue, char **pszValue, int bSetPszValue )
5181 : {
5182 4651 : nc_type nAttrType = NC_NAT;
5183 4651 : size_t nAttrLen = 0;
5184 : size_t nAttrValueSize;
5185 4651 : int status = 0; /*rename this */
5186 : size_t m;
5187 : char szTemp[ NCDF_MAX_STR_LEN ];
5188 4651 : char *pszAttrValue = NULL;
5189 4651 : double dfValue = 0.0;
5190 :
5191 4651 : status = nc_inq_att( nCdfId, nVarId, pszAttrName, &nAttrType, &nAttrLen);
5192 4651 : if ( status != NC_NOERR )
5193 0 : return CE_Failure;
5194 :
5195 : /* Allocate guaranteed minimum size */
5196 4651 : nAttrValueSize = nAttrLen + 1;
5197 4651 : pszAttrValue = (char *) CPLCalloc( nAttrValueSize, sizeof( char ));
5198 4651 : *pszAttrValue = '\0';
5199 :
5200 4651 : if ( nAttrLen > 1 && nAttrType != NC_CHAR )
5201 279 : NCDFSafeStrcat(&pszAttrValue, (char *)"{ ", &nAttrValueSize);
5202 :
5203 4651 : switch (nAttrType) {
5204 : case NC_CHAR:
5205 3029 : nc_get_att_text( nCdfId, nVarId, pszAttrName, pszAttrValue );
5206 3029 : pszAttrValue[nAttrLen]='\0';
5207 3029 : dfValue = 0.0;
5208 3029 : break;
5209 : /* TODO support NC_UBYTE */
5210 : case NC_BYTE:
5211 : signed char *pscTemp;
5212 285 : pscTemp = (signed char *) CPLCalloc( nAttrLen, sizeof( signed char ) );
5213 285 : nc_get_att_schar( nCdfId, nVarId, pszAttrName, pscTemp );
5214 285 : dfValue = (double)pscTemp[0];
5215 289 : for(m=0; m < nAttrLen-1; m++) {
5216 4 : sprintf( szTemp, "%d, ", pscTemp[m] );
5217 4 : NCDFSafeStrcat(&pszAttrValue, szTemp, &nAttrValueSize);
5218 : }
5219 285 : sprintf( szTemp, "%d", pscTemp[m] );
5220 285 : NCDFSafeStrcat(&pszAttrValue, szTemp, &nAttrValueSize);
5221 285 : CPLFree(pscTemp);
5222 285 : break;
5223 : case NC_SHORT:
5224 : short *psTemp;
5225 213 : psTemp = (short *) CPLCalloc( nAttrLen, sizeof( short ) );
5226 213 : nc_get_att_short( nCdfId, nVarId, pszAttrName, psTemp );
5227 213 : dfValue = (double)psTemp[0];
5228 399 : for(m=0; m < nAttrLen-1; m++) {
5229 186 : sprintf( szTemp, "%hd, ", psTemp[m] );
5230 186 : NCDFSafeStrcat(&pszAttrValue, szTemp, &nAttrValueSize);
5231 : }
5232 213 : sprintf( szTemp, "%hd", psTemp[m] );
5233 213 : NCDFSafeStrcat(&pszAttrValue, szTemp, &nAttrValueSize);
5234 213 : CPLFree(psTemp);
5235 213 : break;
5236 : case NC_INT:
5237 : int *pnTemp;
5238 110 : pnTemp = (int *) CPLCalloc( nAttrLen, sizeof( int ) );
5239 110 : nc_get_att_int( nCdfId, nVarId, pszAttrName, pnTemp );
5240 110 : dfValue = (double)pnTemp[0];
5241 173 : for(m=0; m < nAttrLen-1; m++) {
5242 63 : sprintf( szTemp, "%d, ", pnTemp[m] );
5243 63 : NCDFSafeStrcat(&pszAttrValue, szTemp, &nAttrValueSize);
5244 : }
5245 110 : sprintf( szTemp, "%d", pnTemp[m] );
5246 110 : NCDFSafeStrcat(&pszAttrValue, szTemp, &nAttrValueSize);
5247 110 : CPLFree(pnTemp);
5248 110 : break;
5249 : case NC_FLOAT:
5250 : float *pfTemp;
5251 144 : pfTemp = (float *) CPLCalloc( nAttrLen, sizeof( float ) );
5252 144 : nc_get_att_float( nCdfId, nVarId, pszAttrName, pfTemp );
5253 144 : dfValue = (double)pfTemp[0];
5254 156 : for(m=0; m < nAttrLen-1; m++) {
5255 12 : sprintf( szTemp, "%.8g, ", pfTemp[m] );
5256 12 : NCDFSafeStrcat(&pszAttrValue, szTemp, &nAttrValueSize);
5257 : }
5258 144 : sprintf( szTemp, "%.8g", pfTemp[m] );
5259 144 : NCDFSafeStrcat(&pszAttrValue,szTemp, &nAttrValueSize);
5260 144 : CPLFree(pfTemp);
5261 144 : break;
5262 : case NC_DOUBLE:
5263 : double *pdfTemp;
5264 870 : pdfTemp = (double *) CPLCalloc(nAttrLen, sizeof(double));
5265 870 : nc_get_att_double( nCdfId, nVarId, pszAttrName, pdfTemp );
5266 870 : dfValue = pdfTemp[0];
5267 908 : for(m=0; m < nAttrLen-1; m++) {
5268 38 : sprintf( szTemp, "%.16g, ", pdfTemp[m] );
5269 38 : NCDFSafeStrcat(&pszAttrValue, szTemp, &nAttrValueSize);
5270 : }
5271 870 : sprintf( szTemp, "%.16g", pdfTemp[m] );
5272 870 : NCDFSafeStrcat(&pszAttrValue, szTemp, &nAttrValueSize);
5273 870 : CPLFree(pdfTemp);
5274 870 : break;
5275 : default:
5276 : CPLDebug( "GDAL_netCDF", "NCDFGetAttr unsupported type %d for attribute %s",
5277 0 : nAttrType,pszAttrName);
5278 0 : CPLFree( pszAttrValue );
5279 0 : pszAttrValue = NULL;
5280 : break;
5281 : }
5282 :
5283 4651 : if ( nAttrLen > 1 && nAttrType!= NC_CHAR )
5284 279 : NCDFSafeStrcat(&pszAttrValue, (char *)" }", &nAttrValueSize);
5285 :
5286 : // CPLDebug( "GDAL_netCDF", "NCDFGetAttr got %s=%s / %f",
5287 : // pszAttrName,pszAttrValue,dfValue);
5288 :
5289 : /* set return values */
5290 4651 : if ( bSetPszValue == TRUE ) *pszValue = pszAttrValue;
5291 160 : else CPLFree ( pszAttrValue );
5292 4651 : if ( pdfValue ) *pdfValue = dfValue;
5293 :
5294 4651 : return CE_None;
5295 : }
5296 :
5297 :
5298 : /* sets pdfValue to first value found */
5299 160 : CPLErr NCDFGetAttr( int nCdfId, int nVarId, const char *pszAttrName,
5300 : double *pdfValue )
5301 : {
5302 160 : return NCDFGetAttr1( nCdfId, nVarId, pszAttrName, pdfValue, NULL, FALSE );
5303 : }
5304 :
5305 : /* pszValue is the responsibility of the caller and must be freed */
5306 4491 : CPLErr NCDFGetAttr( int nCdfId, int nVarId, const char *pszAttrName,
5307 : char **pszValue )
5308 : {
5309 4491 : return NCDFGetAttr1( nCdfId, nVarId, pszAttrName, FALSE, pszValue, TRUE );
5310 : }
5311 :
5312 :
5313 : /* By default write NC_CHAR, but detect for int/float/double */
5314 206 : CPLErr NCDFPutAttr( int nCdfId, int nVarId,
5315 : const char *pszAttrName, const char *pszValue )
5316 : {
5317 206 : nc_type nAttrType = NC_CHAR;
5318 206 : nc_type nTmpAttrType = NC_CHAR;
5319 206 : size_t nAttrLen = 0;
5320 206 : int status = 0;
5321 : size_t i;
5322 : char szTemp[ NCDF_MAX_STR_LEN ];
5323 206 : char *pszTemp = NULL;
5324 206 : char **papszValues = NULL;
5325 :
5326 206 : int bIsArray = FALSE;
5327 206 : int nValue = 0;
5328 206 : float fValue = 0.0f;
5329 206 : double dfValue = 0.0;
5330 :
5331 206 : if ( EQUAL( pszValue, "" ) )
5332 0 : return CE_Failure;
5333 :
5334 206 : strcpy( szTemp,pszValue );
5335 :
5336 : /* tokenize to find multiple values */
5337 206 : int last_char = strlen(pszValue) - 1;
5338 206 : if ( ( szTemp[0] == '{' ) && ( szTemp[last_char] == '}' ) ) {
5339 6 : bIsArray = TRUE;
5340 6 : szTemp[0] = ' ';
5341 6 : szTemp[last_char] = ' ';
5342 : }
5343 : papszValues = CSLTokenizeString2( szTemp, ",", CSLT_STRIPLEADSPACES |
5344 206 : CSLT_STRIPENDSPACES );
5345 206 : if ( papszValues == NULL )
5346 0 : return CE_Failure;
5347 :
5348 206 : nAttrLen = CSLCount(papszValues);
5349 :
5350 : /* first detect type */
5351 206 : nAttrType = NC_CHAR;
5352 539 : for ( i=0; i<nAttrLen; i++ ) {
5353 333 : nTmpAttrType = NC_CHAR;
5354 333 : errno = 0;
5355 333 : nValue = strtol( papszValues[i], &pszTemp, 10 );
5356 : /* test for int */
5357 : /* TODO test for Byte and short - can this be done safely? */
5358 364 : if ( (errno == 0) && (papszValues[i] != pszTemp) && (*pszTemp == 0) ) {
5359 31 : nTmpAttrType = NC_INT;
5360 : }
5361 : else {
5362 : /* test for double */
5363 302 : errno = 0;
5364 302 : dfValue = strtod( papszValues[i], &pszTemp );
5365 302 : if ( (errno == 0) && (papszValues[i] != pszTemp) && (*pszTemp == 0) ) {
5366 : /* test for float instead of double */
5367 : /* strtof() is C89, which is not available in MSVC */
5368 : /* see if we loose precision if we cast to float and write to char* */
5369 22 : fValue = (float)dfValue;
5370 22 : sprintf( szTemp,"%.8g",fValue);
5371 22 : if ( EQUAL(szTemp, papszValues[i] ) )
5372 19 : nTmpAttrType = NC_FLOAT;
5373 : else
5374 3 : nTmpAttrType = NC_DOUBLE;
5375 : }
5376 : }
5377 333 : if ( nTmpAttrType > nAttrType )
5378 27 : nAttrType = nTmpAttrType;
5379 : }
5380 :
5381 : /* now write the data */
5382 390 : if ( !bIsArray && nAttrType == NC_CHAR ) {
5383 : status = nc_put_att_text( nCdfId, nVarId, pszAttrName,
5384 184 : strlen( pszValue ), pszValue );
5385 184 : NCDF_ERR(status);
5386 : }
5387 : else {
5388 :
5389 22 : switch( nAttrType ) {
5390 : case NC_INT:
5391 : int *pnTemp;
5392 12 : pnTemp = (int *) CPLCalloc( nAttrLen, sizeof( int ) );
5393 28 : for(i=0; i < nAttrLen; i++) {
5394 16 : pnTemp[i] = strtol( papszValues[i], &pszTemp, 10 );
5395 : }
5396 : status = nc_put_att_int( nCdfId, nVarId, pszAttrName,
5397 12 : NC_INT, nAttrLen, pnTemp );
5398 12 : NCDF_ERR(status);
5399 12 : CPLFree(pnTemp);
5400 12 : break;
5401 : case NC_FLOAT:
5402 : float *pfTemp;
5403 8 : pfTemp = (float *) CPLCalloc( nAttrLen, sizeof( float ) );
5404 114 : for(i=0; i < nAttrLen; i++) {
5405 106 : pfTemp[i] = (float)strtod( papszValues[i], &pszTemp );
5406 : }
5407 : status = nc_put_att_float( nCdfId, nVarId, pszAttrName,
5408 8 : NC_FLOAT, nAttrLen, pfTemp );
5409 8 : NCDF_ERR(status);
5410 8 : CPLFree(pfTemp);
5411 8 : break;
5412 : case NC_DOUBLE:
5413 : double *pdfTemp;
5414 2 : pdfTemp = (double *) CPLCalloc( nAttrLen, sizeof( double ) );
5415 5 : for(i=0; i < nAttrLen; i++) {
5416 3 : pdfTemp[i] = strtod( papszValues[i], &pszTemp );
5417 : }
5418 : status = nc_put_att_double( nCdfId, nVarId, pszAttrName,
5419 2 : NC_DOUBLE, nAttrLen, pdfTemp );
5420 2 : NCDF_ERR(status);
5421 2 : CPLFree(pdfTemp);
5422 2 : break;
5423 : default:
5424 0 : if ( papszValues ) CSLDestroy( papszValues );
5425 0 : return CE_Failure;
5426 : break;
5427 : }
5428 : }
5429 :
5430 206 : if ( papszValues ) CSLDestroy( papszValues );
5431 :
5432 206 : return CE_None;
5433 : }
|
