Index: frmts/netcdf/netcdfdataset.cpp =================================================================== --- frmts/netcdf/netcdfdataset.cpp (revision 13399) +++ frmts/netcdf/netcdfdataset.cpp (working copy) @@ -30,8 +30,12 @@ #include "netcdfdataset.h" CPL_CVSID("$Id$"); +static bool approx_equal(double a, double b) +{ + const double epsilon = 1e-5; + return (fabs(a-b) <= epsilon); +} - /************************************************************************/ /* ==================================================================== */ /* netCDFRasterBand */ @@ -43,33 +47,88 @@ nc_type nc_datatype; int nZId; int nZDim; - int nLevel; + int nLevel; int nBandXPos; int nBandYPos; int *panBandZPos; int *panBandZLev; int bNoDataSet; double dfNoDataValue; - CPLErr CreateBandMetadata( ); - - public: + int bApplyScaleOffset; + double dfOffset; + double dfScale; + CPLErr CreateBandMetadata( ); - netCDFRasterBand( netCDFDataset *poDS, - int nZId, - int nZDim, - int nLevel, - int *panBandZLen, - int *panBandPos, - int nBand ); - ~netCDFRasterBand( ); - virtual double GetNoDataValue( int * ); - virtual CPLErr SetNoDataValue( double ); - virtual CPLErr IReadBlock( int, int, void * ); + public: + netCDFRasterBand( netCDFDataset *poDS, + int nZId, + int nZDim, + int nLevel, + int *panBandZLen, + int *panBandPos, + int nBand ); + ~netCDFRasterBand( ); + virtual double GetNoDataValue( int * ); + virtual CPLErr SetNoDataValue( double ); + virtual double GetOffset( int *pbSuccess = NULL ); + virtual CPLErr SetOffset( double dfNewValue ); + virtual double GetScale( int *pbSuccess = NULL ); + virtual CPLErr SetScale( double dfNewValue ); + virtual CPLErr IReadBlock( int, int, void * ); + }; /************************************************************************/ +/* GetOffset() */ +/************************************************************************/ + +double netCDFRasterBand::GetOffset( int *pbSuccess ) + +{ + if( pbSuccess != NULL ) + *pbSuccess = TRUE; + + return dfOffset; +} + +/************************************************************************/ +/* SetOffset() */ +/************************************************************************/ + +CPLErr netCDFRasterBand::SetOffset( double dfNewOffset ) + +{ + dfOffset = dfNewOffset; + return CE_None; +} + +/************************************************************************/ +/* GetScale() */ +/************************************************************************/ + +double netCDFRasterBand::GetScale( int *pbSuccess ) + +{ + if( pbSuccess != NULL ) + *pbSuccess = TRUE; + + return dfScale; +} + +/************************************************************************/ +/* SetScale() */ +/************************************************************************/ + +CPLErr netCDFRasterBand::SetScale( double dfNewScale ) + +{ + dfScale = dfNewScale; + return CE_None; +} + +/************************************************************************/ /* GetMetadata() */ /************************************************************************/ char **netCDFDataset::GetMetadata( const char *pszDomain ) @@ -88,9 +147,9 @@ const char * netCDFDataset::GetProjectionRef() { if( bGotGeoTransform ) - return pszProjection; + return pszProjection; else - return ""; + return ""; } /************************************************************************/ @@ -125,17 +184,17 @@ netCDFRasterBand::~netCDFRasterBand() { - if( panBandZPos ) - CPLFree( panBandZPos ); + if( panBandZPos ) + CPLFree( panBandZPos ); if( panBandZLev ) - CPLFree( panBandZLev ); + CPLFree( panBandZLev ); } /************************************************************************/ /* CreateBandMetadata() */ /************************************************************************/ -CPLErr netCDFRasterBand::CreateBandMetadata( ) +CPLErr netCDFRasterBand::CreateBandMetadata( ) { char szVarName[NC_MAX_NAME]; char szMetaName[NC_MAX_NAME]; @@ -153,135 +212,135 @@ char szTemp[NC_MAX_NAME]; const char *pszValue; - nc_type nVarType; + nc_type nVarType; netCDFDataset *poDS; poDS = (netCDFDataset *) this->poDS; -/* -------------------------------------------------------------------- */ -/* Compute all dimensions from Band number and save in Metadata */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Compute all dimensions from Band number and save in Metadata */ + /* -------------------------------------------------------------------- */ nc_inq_varname( poDS->cdfid, nZId, szVarName ); nc_inq_varndims( poDS->cdfid, nZId, &nd ); -/* -------------------------------------------------------------------- */ -/* Compute multidimention band position */ -/* */ -/* BandPosition = (Total - sum(PastBandLevels) - 1)/sum(remainingLevels)*/ -/* if Data[2,3,4,x,y] */ -/* */ -/* BandPos0 = (nBand ) / (3*4) */ -/* BandPos1 = (nBand - BandPos0*(3*4) ) / (4) */ -/* BandPos2 = (nBand - BandPos0*(3*4) ) % (4) */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Compute multidimention band position */ + /* */ + /* BandPosition = (Total - sum(PastBandLevels) - 1)/sum(remainingLevels)*/ + /* if Data[2,3,4,x,y] */ + /* */ + /* BandPos0 = (nBand ) / (3*4) */ + /* BandPos1 = (nBand - BandPos0*(3*4) ) / (4) */ + /* BandPos2 = (nBand - BandPos0*(3*4) ) % (4) */ + /* -------------------------------------------------------------------- */ sprintf( szMetaName,"NETCDF_VARNAME"); sprintf( szMetaTemp,"%s",szVarName); SetMetadataItem( szMetaName, szMetaTemp ); if( nd == 3 ) { - Sum *= panBandZLev[0]; + Sum *= panBandZLev[0]; } for( i=0; i < nd-2 ; i++ ) { - if( i != nd - 2 -1 ) { - for( j=i+1; j < nd-2; j++ ) { - Sum *= panBandZLev[j]; - } - result = (int) ( ( nLevel-Taken ) / Sum ); - } - else { - result = (int) ( ( nLevel-Taken ) % Sum ); - } + if( i != nd - 2 -1 ) { + for( j=i+1; j < nd-2; j++ ) { + Sum *= panBandZLev[j]; + } + result = (int) ( ( nLevel-Taken ) / Sum ); + } + else { + result = (int) ( ( nLevel-Taken ) % Sum ); + } - strcpy(szVarName, poDS->papszDimName[poDS->paDimIds[ - panBandZPos[i]]] ); + strcpy(szVarName, poDS->papszDimName[poDS->paDimIds[ + panBandZPos[i]]] ); - sprintf( szMetaName,"NETCDF_DIMENSION_%s", szVarName ); + sprintf( szMetaName,"NETCDF_DIMENSION_%s", szVarName ); - status=nc_inq_varid(poDS->cdfid, - szVarName, - &nVarID ); + status=nc_inq_varid(poDS->cdfid, + szVarName, + &nVarID ); -/* -------------------------------------------------------------------- */ -/* Try to uppercase the first letter of the variable */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Try to uppercase the first letter of the variable */ + /* -------------------------------------------------------------------- */ - if( status != NC_NOERR ) { - szVarName[0]=toupper(szVarName[0]); - status=nc_inq_varid(poDS->cdfid, - szVarName, - &nVarID ); - } + if( status != NC_NOERR ) { + szVarName[0]=toupper(szVarName[0]); + status=nc_inq_varid(poDS->cdfid, + szVarName, + &nVarID ); + } - status = nc_inq_vartype( poDS->cdfid, nVarID, &nVarType ); + status = nc_inq_vartype( poDS->cdfid, nVarID, &nVarType ); - nDims = 0; - status = nc_inq_varndims( poDS->cdfid, nVarID, &nDims ); + nDims = 0; + status = nc_inq_varndims( poDS->cdfid, nVarID, &nDims ); - if( nDims == 1 ) { - count[0]=1; - start[0]=result; - switch( nVarType ) { - case NC_SHORT: - short sData; - status = nc_get_vara_short( poDS->cdfid, nVarID, - start, - count, &sData ); - sprintf( szMetaTemp,"%d", sData ); - break; - case NC_INT: - int nData; - status = nc_get_vara_int( poDS->cdfid, nVarID, - start, - count, &nData ); - sprintf( szMetaTemp,"%d", nData ); - break; - case NC_FLOAT: - float fData; - status = nc_get_vara_float( poDS->cdfid, nVarID, - start, - count, &fData ); - sprintf( szMetaTemp,"%f", fData ); - break; - case NC_DOUBLE: - double dfData; - status = nc_get_vara_double( poDS->cdfid, nVarID, - start, - count, &dfData); - sprintf( szMetaTemp,"%g", dfData ); - break; - default: - break; - } - } - else - sprintf( szMetaTemp,"%d", result+1); - + if( nDims == 1 ) { + count[0]=1; + start[0]=result; + switch( nVarType ) { + case NC_SHORT: + short sData; + status = nc_get_vara_short( poDS->cdfid, nVarID, + start, + count, &sData ); + sprintf( szMetaTemp,"%d", sData ); + break; + case NC_INT: + int nData; + status = nc_get_vara_int( poDS->cdfid, nVarID, + start, + count, &nData ); + sprintf( szMetaTemp,"%d", nData ); + break; + case NC_FLOAT: + float fData; + status = nc_get_vara_float( poDS->cdfid, nVarID, + start, + count, &fData ); + sprintf( szMetaTemp,"%f", fData ); + break; + case NC_DOUBLE: + double dfData; + status = nc_get_vara_double( poDS->cdfid, nVarID, + start, + count, &dfData); + sprintf( szMetaTemp,"%g", dfData ); + break; + default: + break; + } + } + else + sprintf( szMetaTemp,"%d", result+1); + SetMetadataItem( szMetaName, szMetaTemp ); -/* -------------------------------------------------------------------- */ -/* Fetch dimension units */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Fetch dimension units */ + /* -------------------------------------------------------------------- */ - strcpy( szTemp, szVarName ); - strcat( szTemp, "#units" ); - pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); - if( pszValue != NULL ) { - if( EQUAL( pszValue, "T") ) { - strcpy( szTemp, szVarName ); - strcat( szTemp, "#original_units" ); - pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); - strcpy( szTemp, "NETCDF_"); - strcat( szTemp, szVarName ); - strcat( szTemp, "_original_units" ); - SetMetadataItem( szTemp, pszValue ); - } - else { - strcpy( szTemp, "NETCDF_"); - strcat( szTemp, szVarName ); - strcat( szTemp, "_units" ); - SetMetadataItem( szTemp, pszValue ); - } - } - Taken += result * Sum; + strcpy( szTemp, szVarName ); + strcat( szTemp, "#units" ); + pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); + if( pszValue != NULL ) { + if( EQUAL( pszValue, "T") ) { + strcpy( szTemp, szVarName ); + strcat( szTemp, "#original_units" ); + pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); + strcpy( szTemp, "NETCDF_"); + strcat( szTemp, szVarName ); + strcat( szTemp, "_original_units" ); + SetMetadataItem( szTemp, pszValue ); + } + else { + strcpy( szTemp, "NETCDF_"); + strcat( szTemp, szVarName ); + strcat( szTemp, "_units" ); + SetMetadataItem( szTemp, pszValue ); + } + } + Taken += result * Sum; } @@ -293,17 +352,18 @@ /* netCDFRasterBand() */ /************************************************************************/ -netCDFRasterBand::netCDFRasterBand( netCDFDataset *poDS, - int nZId, - int nZDim, - int nLevel, - int *panBandZLev, - int *panBandDimPos, - int nBand) +netCDFRasterBand::netCDFRasterBand( netCDFDataset *poDS, + int nZId, + int nZDim, + int nLevel, + int *panBandZLev, + int *panBandDimPos, + int nBand) { double dfNoData; int bNoDataSet = FALSE; + int bApplyScaleOffset = FALSE; nc_type vartype=NC_NAT; nc_type atttype=NC_NAT; size_t attlen; @@ -321,19 +381,19 @@ this->nBandXPos = panBandDimPos[0]; this->nBandYPos = panBandDimPos[1]; -/* -------------------------------------------------------------------- */ -/* Take care of all other dimmensions */ -/* ------------------------------------------------------------------ */ + /* -------------------------------------------------------------------- */ + /* Take care of all other dimmensions */ + /* ------------------------------------------------------------------ */ if( nZDim > 2 ) { - this->panBandZPos = - (int *) CPLCalloc( nZDim-2, sizeof( int ) ); - this->panBandZLev = - (int *) CPLCalloc( nZDim-2, sizeof( int ) ); + this->panBandZPos = + (int *) CPLCalloc( nZDim-2, sizeof( int ) ); + this->panBandZLev = + (int *) CPLCalloc( nZDim-2, sizeof( int ) ); - for ( int i=0; i < nZDim - 2; i++ ){ - this->panBandZPos[i] = panBandDimPos[i+2]; - this->panBandZLev[i] = panBandZLev[i]; - } + for ( int i=0; i < nZDim - 2; i++ ){ + this->panBandZPos[i] = panBandDimPos[i+2]; + this->panBandZLev[i] = panBandZLev[i]; + } } CreateBandMetadata(); bNoDataSet = FALSE; @@ -342,13 +402,14 @@ nBlockXSize = poDS->GetRasterXSize( ); nBlockYSize = 1; -/* -------------------------------------------------------------------- */ -/* Get the type of the "z" variable, our target raster array. */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Get the type of the "z" variable, our target raster array. */ + /* -------------------------------------------------------------------- */ if( nc_inq_var( poDS->cdfid, nZId, NULL, &nc_datatype, NULL, NULL, - NULL ) != NC_NOERR ){ - CPLError( CE_Failure, CPLE_AppDefined, - "Error in nc_var_inq() on 'z'." ); + NULL ) != NC_NOERR ) + { + CPLError( CE_Failure, CPLE_AppDefined, + "Error in nc_var_inq() on 'z'." ); return; } @@ -365,103 +426,125 @@ else { if( nBand == 1 ) - CPLError( CE_Warning, CPLE_AppDefined, - "Unsupported netCDF datatype (%d), treat as Float32.", - (int) nc_datatype ); + CPLError( CE_Warning, CPLE_AppDefined, + "Unsupported netCDF datatype (%d), treat as Float32.", + (int) nc_datatype ); eDataType = GDT_Float32; } -/* -------------------------------------------------------------------- */ -/* Find out what is No Data for this variable */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Find out what is No Data for this variable */ + /* -------------------------------------------------------------------- */ - status = nc_inq_att( poDS->cdfid, nZId, - _FillValue, &atttype, &attlen); + status = nc_inq_att( poDS->cdfid, nZId, + _FillValue, &atttype, &attlen); -/* -------------------------------------------------------------------- */ -/* Look for either Missing_Value or _FillValue attributes */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Look for either Missing_Value or _FillValue attributes */ + /* -------------------------------------------------------------------- */ if( status == NC_NOERR ) { - strcpy(szNoValueName, _FillValue ); + strcpy(szNoValueName, _FillValue ); } else { - status = nc_inq_att( poDS->cdfid, nZId, - "missing_value", &atttype, &attlen ); - if( status == NC_NOERR ) { + status = nc_inq_att( poDS->cdfid, nZId, + "missing_value", &atttype, &attlen ); + if( status == NC_NOERR ) { - strcpy( szNoValueName, "missing_value" ); - } + strcpy( szNoValueName, "missing_value" ); + } } nc_inq_vartype( poDS->cdfid, nZId, &vartype ); if( status == NC_NOERR ) { - switch( atttype ) { - case NC_CHAR: - char *fillc; - fillc = (char *) CPLCalloc( attlen+1, sizeof(char) ); - status=nc_get_att_text( poDS->cdfid, nZId, - szNoValueName, fillc ); - dfNoData = atof( fillc ); - CPLFree(fillc); - break; - case NC_SHORT: - short sNoData; - status = nc_get_att_short( poDS->cdfid, nZId, - szNoValueName, &sNoData ); - dfNoData = (double) sNoData; - break; - case NC_INT: - int nNoData; - status = nc_get_att_int( poDS->cdfid, nZId, - szNoValueName, &nNoData ); - dfNoData = (double) nNoData; - break; - case NC_FLOAT: - float fNoData; - status = nc_get_att_float( poDS->cdfid, nZId, - szNoValueName, &fNoData ); - dfNoData = (double) fNoData; - break; - case NC_DOUBLE: - status = nc_get_att_double( poDS->cdfid, nZId, - szNoValueName, &dfNoData ); - break; - default: - break; - } - status = nc_get_att_double( poDS->cdfid, nZId, - szNoValueName, &dfNoData ); - + switch( atttype ) { + case NC_CHAR: + char *fillc; + fillc = (char *) CPLCalloc( attlen+1, sizeof(char) ); + status=nc_get_att_text( poDS->cdfid, nZId, + szNoValueName, fillc ); + dfNoData = atof( fillc ); + CPLFree(fillc); + break; + case NC_SHORT: + short sNoData; + status = nc_get_att_short( poDS->cdfid, nZId, + szNoValueName, &sNoData ); + dfNoData = (double) sNoData; + break; + case NC_INT: + int nNoData; + status = nc_get_att_int( poDS->cdfid, nZId, + szNoValueName, &nNoData ); + dfNoData = (double) nNoData; + break; + case NC_FLOAT: + float fNoData; + status = nc_get_att_float( poDS->cdfid, nZId, + szNoValueName, &fNoData ); + dfNoData = (double) fNoData; + break; + case NC_DOUBLE: + status = nc_get_att_double( poDS->cdfid, nZId, + szNoValueName, &dfNoData ); + break; + default: + break; + } + status = nc_get_att_double( poDS->cdfid, nZId, + szNoValueName, &dfNoData ); + } else { - switch( vartype ) { - case NC_BYTE: - /* don't do default fill-values for bytes, too risky */ - dfNoData = 0.0; - break; - case NC_CHAR: - dfNoData = NC_FILL_CHAR; - break; - case NC_SHORT: - dfNoData = NC_FILL_SHORT; - break; - case NC_INT: - dfNoData = NC_FILL_INT; - break; - case NC_FLOAT: - dfNoData = NC_FILL_FLOAT; - break; - case NC_DOUBLE: - dfNoData = NC_FILL_DOUBLE; - break; - default: - dfNoData = 0.0; - break; - } - bNoDataSet = TRUE; + switch( vartype ) { + case NC_BYTE: + /* don't do default fill-values for bytes, too risky */ + dfNoData = 0.0; + break; + case NC_CHAR: + dfNoData = NC_FILL_CHAR; + break; + case NC_SHORT: + dfNoData = NC_FILL_SHORT; + break; + case NC_INT: + dfNoData = NC_FILL_INT; + break; + case NC_FLOAT: + dfNoData = NC_FILL_FLOAT; + break; + case NC_DOUBLE: + dfNoData = NC_FILL_DOUBLE; + break; + default: + dfNoData = 0.0; + break; + } } + /* -------------------------------------------------------------------- */ + /* SH: Try to find and offset and scale_factor. */ + /* eDataType should be set to float when bApplyScaleOffset */ + /* -------------------------------------------------------------------- */ + if ( nc_inq_attid ( poDS->cdfid, nZId, "add_offset", NULL) == NC_NOERR){ + double add_offset=0.0; + double scale_factor=1.0; + nc_get_att_double( poDS->cdfid, nZId, "add_offset", &add_offset ); + + if ( nc_inq_attid ( poDS->cdfid, nZId, "scale_factor", NULL) == NC_NOERR){ + nc_get_att_double( poDS->cdfid, nZId, "scale_factor", &scale_factor ); + } + SetOffset( add_offset ); + SetScale( scale_factor ); + // Override the data type to float + if ( eDataType != GDT_Float64 ) + eDataType = GDT_Float32; + // A note on missing value. The original fill value is kept. In the unfortunate + // case that a value is the missing value after the application of scale and + // offset, then you have a problem. + bApplyScaleOffset = TRUE; + } SetNoDataValue( dfNoData ); - + bNoDataSet = TRUE; + } /************************************************************************/ @@ -469,7 +552,7 @@ /************************************************************************/ CPLErr netCDFRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff, - void * pImage ) + void * pImage ) { int nErr=-1; @@ -481,20 +564,23 @@ int Sum=-1; int Taken=-1; int nd; + ouble dfScale, dfOffset; + nc_type nVarType; *pszName='\0'; memset( start, 0, sizeof( start ) ); memset( edge, 0, sizeof( edge ) ); nc_inq_varndims ( cdfid, nZId, &nd ); + nc_inq_vartype( cdfid, nZId, &nVarType ); -/* -------------------------------------------------------------------- */ -/* Locate X, Y and Z position in the array */ -/* -------------------------------------------------------------------- */ - + /* -------------------------------------------------------------------- */ + /* Locate X, Y and Z position in the array */ + /* -------------------------------------------------------------------- */ + start[nBandXPos] = 0; // x dim can move arround in array start[nBandYPos] = nBlockYOff; // y - - edge[nBandXPos] = nBlockXSize; + + edge[nBandXPos] = nBlockXSize; edge[nBandYPos] = 1; if( nd == 3 ) { @@ -502,70 +588,159 @@ edge [panBandZPos[0]] = 1; } -/* -------------------------------------------------------------------- */ -/* Compute multidimention band position */ -/* */ -/* BandPosition = (Total - sum(PastBandLevels) - 1)/sum(remainingLevels)*/ -/* if Data[2,3,4,x,y] */ -/* */ -/* BandPos0 = (nBand ) / (3*4) */ -/* BandPos1 = (nBand - (3*4) ) / (4) */ -/* BandPos2 = (nBand - (3*4) ) % (4) */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Compute multidimention band position */ + /* */ + /* BandPosition = (Total - sum(PastBandLevels) - 1)/sum(remainingLevels)*/ + /* if Data[2,3,4,x,y] */ + /* */ + /* BandPos0 = (nBand ) / (3*4) */ + /* BandPos1 = (nBand - (3*4) ) / (4) */ + /* BandPos2 = (nBand - (3*4) ) % (4) */ + /* -------------------------------------------------------------------- */ if( nd > 3 ) { - for( i=0; i < nd-2-1 ; i++ ) { - Sum = 1; - Taken = 0; - for( j=i+1; j < nd-2; j++ ) { - Sum *= panBandZLev[j]; - } - start[panBandZPos[i]] = ( nLevel-Taken ) / Sum; - edge[panBandZPos[i]] = 1; - Taken += Sum; - } - start[panBandZPos[i]] = ( nLevel-( Taken-Sum ) ) % Sum; - edge[panBandZPos[i]] = 1; + for( i=0; i < nd-2-1 ; i++ ) { + Sum = 1; + Taken = 0; + for( j=i+1; j < nd-2; j++ ) { + Sum *= panBandZLev[j]; + } + start[panBandZPos[i]] = ( nLevel-Taken ) / Sum; + edge[panBandZPos[i]] = 1; + Taken += Sum; + } + start[panBandZPos[i]] = ( nLevel-( Taken-Sum ) ) % Sum; + edge[panBandZPos[i]] = 1; } + dfScale = GetScale(); + dfOffset = GetOffset(); + if( nVarType == NC_BYTE ) + if ( !bApplyScaleOffset ) + nErr = nc_get_vara_uchar( cdfid, nZId, start, edge, + (unsigned char *) pImage ); + else + { + unsigned char *data_values = new unsigned char[nBlockXSize * nBlockYSize]; - if( eDataType == GDT_Byte ) - nErr = nc_get_vara_uchar( cdfid, nZId, start, edge, - (unsigned char *) pImage ); - else if( eDataType == GDT_Int16 ) - nErr = nc_get_vara_short( cdfid, nZId, start, edge, + nErr = nc_get_vara_uchar( cdfid, nZId, start, edge, + data_values ); + // loop, scale, offset, write to pImage + for (int i = 0; i < nBlockXSize * nBlockYSize; i++) + { + if ( approx_equal(data_values[i], dfNoDataValue) ) + ((float *)pImage)[i] = dfNoDataValue; + else + ((float *)pImage)[i] = data_values[i] * dfScale + dfOffset; + } + delete data_values; + } + else if( nVarType == NC_SHORT ) + if ( !bApplyScaleOffset ) + nErr = nc_get_vara_short( cdfid, nZId, start, edge, (short int *) pImage ); - else if( eDataType == GDT_Int32 ) - { + else + { + short int *data_values = new short int[nBlockXSize * nBlockYSize]; + + nErr = nc_get_vara_short( cdfid, nZId, start, edge, + data_values ); + // loop, scale, offset, write to pImage + for (int i = 0; i < nBlockXSize * nBlockYSize; i++) + { + if ( approx_equal(data_values[i], dfNoDataValue) ) + ((float *)pImage)[i] =dfNoDataValue; + else + ((float *)pImage)[i] = data_values[i] * dfScale + dfOffset; + } + delete data_values; + } + + else if( nVarType == NC_INT ) if( sizeof(long) == 4 ) - nErr = nc_get_vara_long( cdfid, nZId, start, edge, - (long *) pImage ); + if ( !bApplyScaleOffset ) + nErr = nc_get_vara_long( cdfid, nZId, start, edge, + (long *) pImage ); + else + { + long *data_values = new long[nBlockXSize * nBlockYSize]; + nErr = nc_get_vara_long( cdfid, nZId, start, edge, + data_values ); + // loop, scale, offset, write to pImage + for (int i = 0; i < nBlockXSize * nBlockYSize; i++) + { + if ( approx_equal(data_values[i], dfNoDataValue) ) + ((float *)pImage)[i] = dfNoDataValue; + else + ((float *)pImage)[i] = data_values[i] * dfScale + dfOffset; + } + delete data_values; + } else - nErr = nc_get_vara_int( cdfid, nZId, start, edge, - (int *) pImage ); + if ( !bApplyScaleOffset ) + nErr = nc_get_vara_int( cdfid, nZId, start, edge, + (int *) pImage ); + else + { + int *data_values = new int[nBlockXSize * nBlockYSize]; + + nErr = nc_get_vara_int( cdfid, nZId, start, edge, + data_values ); + // loop, scale, offset, write to pImage + for (int i = 0; i < nBlockXSize * nBlockYSize; i++) + { + if ( approx_equal(data_values[i], dfNoDataValue) ) + ((float *)pImage)[i] = dfNoDataValue; + else + ((float *)pImage)[i] = data_values[i] * dfScale + dfOffset; + } + delete data_values; + } + else if( nVarType == NC_FLOAT ) + { + // Always loop over data_values because of the NaN check + float *data_values = new float[nBlockXSize * nBlockYSize]; + + nErr = nc_get_vara_float( cdfid, nZId, start, edge, + data_values ); + // loop, scale, offset, write to pImage + for (int i = 0; i < nBlockXSize * nBlockYSize; i++) + { + if ( approx_equal(data_values[i], dfNoDataValue) || CPLIsNan( data_values[i] ) ) + ((float *)pImage)[i] = dfNoDataValue; + else if ( bApplyScaleOffset ) + ((float *) pImage)[i] = data_values[i] * dfScale + dfOffset; + else + ((float *) pImage)[i] = data_values[i]; + } + delete data_values; } - else if( eDataType == GDT_Float32 ){ - nErr = nc_get_vara_float( cdfid, nZId, start, edge, - (float *) pImage ); - for( i=0; iadfGeoTransform[0] = 0.0; poDS->adfGeoTransform[1] = 1.0; @@ -698,12 +863,12 @@ poDS->adfGeoTransform[4] = 0.0; poDS->adfGeoTransform[5] = 1.0; poDS->pszProjection = NULL; - -/* -------------------------------------------------------------------- */ -/* Look for grid_mapping metadata */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Look for grid_mapping metadata */ + /* -------------------------------------------------------------------- */ + strcpy( szGridMappingValue, "" ); strcpy( szGridMappingName, "" ); @@ -712,321 +877,497 @@ strcat(szTemp,"#grid_mapping"); pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); if( pszValue ) { - strcpy(szGridMappingName,szTemp); - strcpy(szGridMappingValue,pszValue); + strcpy(szGridMappingName,szTemp); + strcpy(szGridMappingValue,pszValue); } -/* -------------------------------------------------------------------- */ -/* Look for dimension: lon */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Look for dimension: lon */ + /* -------------------------------------------------------------------- */ - for( i = 0; (i < strlen( poDS->papszDimName[ poDS->nDimXid ] ) && - i < 3 ); i++ ) { - szDimNameX[i] = tolower( ( poDS->papszDimName[poDS->nDimXid] )[i] ); - } - szDimNameX[3] = '\0'; + for( i = 0; (i < strlen( poDS->papszDimName[ poDS->nDimXid ] ) && + i < 3 ); i++ ) { + szDimNameX[i] = tolower( ( poDS->papszDimName[poDS->nDimXid] )[i] ); + } + szDimNameX[3] = '\0'; -/* -------------------------------------------------------------------- */ -/* Read grid_mappinginformation and set projections */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Read grid_mappinginformation and set projections */ + /* -------------------------------------------------------------------- */ - if( !( EQUAL(szGridMappingName,"" ) ) ) { - nc_inq_varid( cdfid, szGridMappingValue, &nVarProjectionID ); - poDS->ReadAttributes( cdfid, nVarProjectionID ); - - strcpy( szTemp, szGridMappingValue ); - strcat( szTemp, "#" ); - strcat( szTemp, GRD_MAPPING_NAME ); - pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); + if( !( EQUAL(szGridMappingName,"" ) ) ) { + nc_inq_varid( cdfid, szGridMappingValue, &nVarProjectionID ); + poDS->ReadAttributes( cdfid, nVarProjectionID ); - if( pszValue != NULL ) { -/* -------------------------------------------------------------------- */ -/* Transverse Mercator */ -/* -------------------------------------------------------------------- */ + strcpy( szTemp, szGridMappingValue ); + strcat( szTemp, "#" ); + strcat( szTemp, GRD_MAPPING_NAME ); + pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); + papszGeolocation = CSLSetNameValue( papszGeolocation, + "grid_mapping", + pszValue); + if( pszValue != NULL ) { + /* -------------------------------------------------------------------- */ + /* Transverse Mercator */ + /* -------------------------------------------------------------------- */ - if( EQUAL( pszValue, TM ) ) { + if( EQUAL( pszValue, TM ) ) { - dfScale = - poDS->FetchCopyParm( szGridMappingValue, - SCALE_FACTOR, 1.0 ); + strcpy(szTemp,szGridMappingValue); + strcat( szTemp, "#" ); + strcat( szTemp, SCALE_FACTOR ); + pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); + papszGeolocation = CSLSetNameValue( papszGeolocation, + SCALE_FACTOR, + pszValue); + if( pszValue ) + dfScale = atof( pszValue ); - dfCenterLon = - poDS->FetchCopyParm( szGridMappingValue, - LONG_CENTRAL_MERIDIAN, 0.0 ); + strcpy(szTemp,szGridMappingValue); + strcat( szTemp, "#" ); + strcat( szTemp, LONG_CENTRAL_MERIDIAN ); + papszGeolocation = CSLSetNameValue( papszGeolocation, + LONG_CENTRAL_MERIDIAN, + pszValue); + pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); + if( pszValue ) + dfCenterLon = atof( pszValue ); - dfCenterLat = - poDS->FetchCopyParm( szGridMappingValue, - LAT_PROJ_ORIGIN, 0.0 ); + strcpy(szTemp,szGridMappingValue); + strcat( szTemp, "#" ); + strcat( szTemp, LAT_PROJ_ORIGIN ); + papszGeolocation = CSLSetNameValue( papszGeolocation, + LAT_PROJ_ORIGIN, + pszValue); + pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); + if( pszValue ) + dfCenterLat = atof( pszValue ); - dfFalseEasting = - poDS->FetchCopyParm( szGridMappingValue, - FALSE_EASTING, 0.0 ); + strcpy(szTemp,szGridMappingValue); + strcat( szTemp, "#" ); + strcat( szTemp, FALSE_EASTING ); + papszGeolocation = CSLSetNameValue( papszGeolocation, + FALSE_EASTING, + pszValue); + pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); + if( pszValue ) + dfFalseEasting = atof( pszValue ); - dfFalseNorthing = - poDS->FetchCopyParm( szGridMappingValue, - FALSE_NORTHING, 0.0 ); + strcpy(szTemp,szGridMappingValue); + strcat( szTemp, "#" ); + strcat( szTemp, FALSE_NORTHING ); + papszGeolocation = CSLSetNameValue( papszGeolocation, + FALSE_NORTHING, + pszValue); + pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); + if( pszValue ) + dfFalseNorthing = atof (pszValue); - oSRS.SetTM( dfCenterLat, - dfCenterLon, - dfScale, - dfFalseEasting, - dfFalseNorthing ); - oSRS.SetWellKnownGeogCS( "WGS84" ); - } -/* -------------------------------------------------------------------- */ -/* Lambert conformal conic */ -/* -------------------------------------------------------------------- */ - else if( EQUAL( pszValue, L_C_CONIC ) ) { - - dfStdP1 = - poDS->FetchCopyParm( szGridMappingValue, - STD_PARALLEL_1, 0.0 ); + oSRS.SetTM( dfCenterLat, + dfCenterLon, + dfScale, + dfFalseEasting, + dfFalseNorthing ); + oSRS.SetWellKnownGeogCS( "WGS84" ); + } + /* -------------------------------------------------------------------- */ + /* Lambert conformal conic */ + /* -------------------------------------------------------------------- */ + else if( EQUAL( pszValue, L_C_CONIC ) ) { - dfStdP2 = - poDS->FetchCopyParm( szGridMappingValue, - STD_PARALLEL_2, 0.0 ); + strcpy( szTemp,szGridMappingValue ); + strcat( szTemp, "#" ); + strcat( szTemp, STD_PARALLEL_1 ); + pszValue = CSLFetchNameValue( poDS->papszMetadata, szTemp ); + papszGeolocation = CSLSetNameValue( papszGeolocation, + STD_PARALLEL_1, + pszValue); - dfCenterLon = - poDS->FetchCopyParm( szGridMappingValue, - LONG_CENTRAL_MERIDIAN, 0.0 ); - dfCenterLat = - poDS->FetchCopyParm( szGridMappingValue, - LAT_PROJ_ORIGIN, 0.0 ); + if( pszValue ) { + dfStdP1 = atof( pszValue ); + } - dfFalseEasting = - poDS->FetchCopyParm( szGridMappingValue, - FALSE_EASTING, 0.0 ); + strcpy( szTemp,szGridMappingValue ); + strcat( szTemp, "#" ); + strcat( szTemp, STD_PARALLEL_1 ); + pszValue = CSLFetchNameValue( poDS->papszMetadata, szTemp ); + papszGeolocation = CSLSetNameValue( papszGeolocation, + STD_PARALLEL_2, + pszValue); - dfFalseNorthing = - poDS->FetchCopyParm( szGridMappingValue, - FALSE_NORTHING, 0.0 ); - - oSRS.SetLCC( dfStdP1, dfStdP2, dfCenterLat, dfCenterLon, - dfFalseEasting, dfFalseNorthing ); - oSRS.SetWellKnownGeogCS( "WGS84" ); + if( pszValue ) { + dfStdP2 = atof( pszValue ); + } - } -/* -------------------------------------------------------------------- */ -/* Is this Latitude/Longitude Grid */ -/* -------------------------------------------------------------------- */ - - } else if( EQUAL( szDimNameX,"lon" ) ) { - oSRS.SetWellKnownGeogCS( "WGS84" ); + strcpy( szTemp,szGridMappingValue ); + strcat( szTemp, "#" ); + strcat( szTemp, LONG_CENTRAL_MERIDIAN ); + pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); + papszGeolocation = CSLSetNameValue( papszGeolocation, + LONG_CENTRAL_MERIDIAN, + pszValue); - } else { + if( pszValue ) + dfCenterLon = atof(pszValue); + + + strcpy(szTemp,szGridMappingValue); + strcat( szTemp, "#" ); + strcat( szTemp, LAT_PROJ_ORIGIN ); + pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); + papszGeolocation = CSLSetNameValue( papszGeolocation, + LAT_PROJ_ORIGIN, + pszValue); + + if( pszValue ) + dfCenterLat = atof(pszValue); + + strcpy(szTemp,szGridMappingValue); + strcat( szTemp, "#" ); + strcat( szTemp, FALSE_EASTING ); + pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); + papszGeolocation = CSLSetNameValue( papszGeolocation, + FALSE_EASTING, + pszValue); + + if( pszValue ) + dfFalseEasting = atof(pszValue); + + strcpy(szTemp,szGridMappingValue); + strcat( szTemp, "#" ); + strcat( szTemp, FALSE_NORTHING ); + pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); + papszGeolocation = CSLSetNameValue( papszGeolocation, + FALSE_NORTHING, + pszValue); + + if( pszValue ) + dfFalseNorthing = atof(pszValue); + + oSRS.SetLCC( dfStdP1, dfStdP2, dfCenterLat, dfCenterLon, + dfFalseEasting, dfFalseNorthing ); + oSRS.SetWellKnownGeogCS( "WGS84" ); + + } + /* -------------------------------------------------------------------- */ + /* Is this Latitude/Longitude Grid */ + /* -------------------------------------------------------------------- */ + + } else if( EQUAL( szDimNameX,"lon" ) ) { + papszGeolocation = CSLSetNameValue( papszGeolocation, + "grid_mapping", + "Latitude Longitude" ); + + strcpy( szTemp, poDS->papszDimName[nDimXid] ); + strcat( szTemp, "#" ); + strcat( szTemp, UNITS ); + + pszValue = CSLFetchNameValue( poDS->papszMetadata, szTemp ); + + strcpy( szTemp, "Lon_"); + strcat( szTemp, UNITS ); + + papszGeolocation = CSLSetNameValue( papszGeolocation, + szTemp, + pszValue); + + strcpy(szTemp, poDS->papszDimName[nDimXid]); + strcat(szTemp, "#" ); + strcat(szTemp, AXIS ); + + pszValue = CSLFetchNameValue( poDS->papszMetadata, szTemp ); + + strcpy( szTemp, "Lon_"); + strcat( szTemp, AXIS ); + + papszGeolocation = CSLSetNameValue( papszGeolocation, + szTemp, + pszValue); + + strcpy(szTemp, poDS->papszDimName[nDimXid]); + strcat(szTemp, "#" ); + strcat(szTemp, BOUNDS ); + + pszValue = CSLFetchNameValue( poDS->papszMetadata, szTemp ); + + strcpy( szTemp, "Lon_"); + strcat( szTemp, BOUNDS ); + + papszGeolocation = CSLSetNameValue( papszGeolocation, + szTemp, + pszValue); + strcpy(szTemp, poDS->papszDimName[nDimXid]); + strcat(szTemp, "#" ); + strcat(szTemp, ORIG_AXIS ); + + pszValue = CSLFetchNameValue( poDS->papszMetadata, szTemp ); + + strcpy( szTemp, "Lon_"); + strcat( szTemp, ORIG_AXIS ); + + papszGeolocation = CSLSetNameValue( papszGeolocation, + szTemp, + pszValue); + + + strcpy( szTemp, poDS->papszDimName[nDimYid] ); + strcat( szTemp, "#" ); + strcat( szTemp, UNITS ); + + pszValue = CSLFetchNameValue( poDS->papszMetadata, szTemp ); + + strcpy( szTemp, "Lat_"); + strcat( szTemp, UNITS ); + + papszGeolocation = CSLSetNameValue( papszGeolocation, + szTemp, + pszValue); + + strcpy(szTemp, poDS->papszDimName[nDimYid]); + strcat(szTemp, "#" ); + strcat(szTemp, AXIS ); + + pszValue = CSLFetchNameValue( poDS->papszMetadata, szTemp ); + + strcpy( szTemp, "Lat_"); + strcat( szTemp, AXIS ); + + papszGeolocation = CSLSetNameValue( papszGeolocation, + szTemp, + pszValue); + + strcpy(szTemp, poDS->papszDimName[nDimYid]); + strcat(szTemp, "#" ); + strcat(szTemp, BOUNDS ); + + pszValue = CSLFetchNameValue( poDS->papszMetadata, szTemp ); + + strcpy( szTemp, "Lat_"); + strcat( szTemp, BOUNDS ); + + papszGeolocation = CSLSetNameValue( papszGeolocation, + szTemp, + pszValue) + ; + strcpy(szTemp, poDS->papszDimName[nDimYid]); + strcat(szTemp, "#" ); + strcat(szTemp, ORIG_AXIS ); + + pszValue = CSLFetchNameValue( poDS->papszMetadata, szTemp ); + + strcpy( szTemp, "Lat_"); + strcat( szTemp, ORIG_AXIS ); + + papszGeolocation = CSLSetNameValue( papszGeolocation, + szTemp, + pszValue); + + + + oSRS.SetWellKnownGeogCS( "WGS84" ); + } else { // This would be too indiscrimant. But we should set // it if we know the data is geographic. //oSRS.SetWellKnownGeogCS( "WGS84" ); - } - } -/* -------------------------------------------------------------------- */ -/* Read projection coordinates */ -/* -------------------------------------------------------------------- */ + } + } + /* -------------------------------------------------------------------- */ + /* Read projection coordinates */ + /* -------------------------------------------------------------------- */ - nc_inq_varid( cdfid, poDS->papszDimName[nDimXid], &nVarDimXID ); - nc_inq_varid( cdfid, poDS->papszDimName[nDimYid], &nVarDimYID ); - - if( ( nVarDimXID != -1 ) && ( nVarDimYID != -1 ) ) { - pdfXCoord = (double *) CPLCalloc( xdim, sizeof(double) ); - pdfYCoord = (double *) CPLCalloc( ydim, sizeof(double) ); - -/* -------------------------------------------------------------------- */ -/* Is pixel spacing is uniform accross the map? */ -/* -------------------------------------------------------------------- */ - start[0] = 0; - edge[0] = xdim; - - status = nc_get_vara_double( cdfid, nVarDimXID, - start, edge, pdfXCoord); - edge[0] = ydim; - status = nc_get_vara_double( cdfid, nVarDimYID, - start, edge, pdfYCoord); + nc_inq_varid( cdfid, poDS->papszDimName[nDimXid], &nVarDimXID ); + nc_inq_varid( cdfid, poDS->papszDimName[nDimYid], &nVarDimYID ); -/* -------------------------------------------------------------------- */ -/* Check Longitude */ -/* -------------------------------------------------------------------- */ + if( ( nVarDimXID != -1 ) && ( nVarDimYID != -1 ) ) { + pdfXCoord = (double *) CPLCalloc( xdim, sizeof(double) ); + pdfYCoord = (double *) CPLCalloc( ydim, sizeof(double) ); - nSpacingBegin = (int) poDS->rint((pdfXCoord[1]-pdfXCoord[0]) * 1000); - - nSpacingMiddle = (int) poDS->rint((pdfXCoord[xdim / 2] - - pdfXCoord[(xdim / 2) + 1]) * 1000); - - nSpacingLast = (int) poDS->rint((pdfXCoord[xdim - 2] - - pdfXCoord[xdim-1]) * 1000); - - if( ( abs( nSpacingBegin ) == abs( nSpacingLast ) ) && - ( abs( nSpacingBegin ) == abs( nSpacingMiddle ) ) && - ( abs( nSpacingMiddle ) == abs( nSpacingLast ) ) ) { + /* -------------------------------------------------------------------- */ + /* Is pixel spacing is uniform accross the map? */ + /* -------------------------------------------------------------------- */ + start[0] = 0; + edge[0] = xdim; -/* -------------------------------------------------------------------- */ -/* Longitude is equaly spaced, check lattitde */ -/* -------------------------------------------------------------------- */ - nSpacingBegin = (int) poDS->rint((pdfYCoord[1]-pdfYCoord[0]) * - 1000); - - nSpacingMiddle = (int) poDS->rint((pdfYCoord[ydim / 2] - - pdfYCoord[(ydim / 2) + 1]) * - 1000); - - nSpacingLast = (int) poDS->rint((pdfYCoord[ydim - 2] - - pdfYCoord[ydim-1]) * - 1000); + status = nc_get_vara_double( cdfid, nVarDimXID, + start, edge, pdfXCoord); + edge[0] = ydim; + status = nc_get_vara_double( cdfid, nVarDimYID, + start, edge, pdfYCoord); - -/* -------------------------------------------------------------------- */ -/* For Latitude we allow an error of 0.1 degrees for gaussion */ -/* gridding */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Check Longitude */ + /* -------------------------------------------------------------------- */ - if((( abs( abs(nSpacingBegin) - abs(nSpacingLast) ) ) < 100 ) && - (( abs( abs(nSpacingBegin) - abs(nSpacingMiddle) ) ) < 100 ) && - (( abs( abs(nSpacingMiddle) - abs(nSpacingLast) ) ) < 100) ) { + nSpacingBegin = (int) poDS->rint((pdfXCoord[1]-pdfXCoord[0]) * 1000); - if( ( abs( nSpacingBegin ) != abs( nSpacingLast ) ) || - ( abs( nSpacingBegin ) != abs( nSpacingMiddle ) ) || - ( abs( nSpacingMiddle ) != abs( nSpacingLast ) ) ) { - - CPLError(CE_Warning, 1,"Latitude grid not spaced evenly.\nSeting projection for grid spacing is within 0.1 degrees threshold.\n"); + nSpacingMiddle = (int) poDS->rint((pdfXCoord[xdim / 2] - + pdfXCoord[(xdim / 2) + 1]) * 1000); - } -/* -------------------------------------------------------------------- */ -/* We have gridded data s we can set the Gereferencing info. */ -/* -------------------------------------------------------------------- */ + nSpacingLast = (int) poDS->rint((pdfXCoord[xdim - 2] - + pdfXCoord[xdim-1]) * 1000); -/* -------------------------------------------------------------------- */ -/* Enable GeoTransform */ -/* -------------------------------------------------------------------- */ - poDS->bGotGeoTransform = TRUE; - - poDS->adfGeoTransform[0] = pdfXCoord[0]; - poDS->adfGeoTransform[3] = pdfYCoord[0]; - poDS->adfGeoTransform[2] = 0; - poDS->adfGeoTransform[4] = 0; - poDS->adfGeoTransform[1] = (( pdfXCoord[xdim-1] - - pdfXCoord[0] ) / - ( poDS->nRasterXSize - 1 )); + if( ( abs( nSpacingBegin ) == abs( nSpacingLast )) && + ( abs( nSpacingBegin ) == abs( nSpacingMiddle )) && + ( abs( nSpacingMiddle ) == abs( nSpacingLast )) ) { - poDS->adfGeoTransform[5] = (( pdfYCoord[ydim-1] - - pdfYCoord[0] ) / - ( poDS->nRasterYSize - 1 )); -/* -------------------------------------------------------------------- */ -/* Compute the center of the pixel */ -/* -------------------------------------------------------------------- */ - poDS->adfGeoTransform[0] = pdfXCoord[0] - - (poDS->adfGeoTransform[1] / 2); + /* -------------------------------------------------------------------- */ + /* Longitude is equaly spaced, check lattitde */ + /* -------------------------------------------------------------------- */ + nSpacingBegin = (int) poDS->rint((pdfYCoord[1]-pdfYCoord[0]) * + 1000); - poDS->adfGeoTransform[3] = pdfYCoord[0] - - (poDS->adfGeoTransform[5] / 2); + nSpacingMiddle = (int) poDS->rint((pdfYCoord[ydim / 2] - + pdfYCoord[(ydim / 2) + 1]) * + 1000); - oSRS.exportToWkt( &(poDS->pszProjection) ); - } - } + nSpacingLast = (int) poDS->rint((pdfYCoord[ydim - 2] - + pdfYCoord[ydim-1]) * + 1000); - CPLFree( pdfXCoord ); - CPLFree( pdfYCoord ); - } + if( ( abs( nSpacingBegin ) == abs( nSpacingLast )) && + ( abs( nSpacingBegin ) == abs( nSpacingMiddle )) && + ( abs( nSpacingMiddle ) == abs( nSpacingLast )) ) { + /* -------------------------------------------------------------------- */ + /* We have gridded data s we can set the Gereferencing info. */ + /* -------------------------------------------------------------------- */ -/* -------------------------------------------------------------------- */ -/* Is this a netCDF file created by GDAL? */ -/* -------------------------------------------------------------------- */ - if( !EQUAL( szGridMappingValue, "" ) ) { - strcpy( szTemp,szGridMappingValue); - strcat( szTemp, "#" ); - strcat( szTemp, "spatial_ref"); - pszWKT = CSLFetchNameValue(poDS->papszMetadata, szTemp); - - if( pszWKT != NULL ) { - - pszProjection = strdup( pszWKT ); - - strcpy(szTemp,szGridMappingValue); - strcat( szTemp, "#" ); - strcat( szTemp, "GeoTransform"); - - pszGeoTransform = CSLFetchNameValue(poDS->papszMetadata, szTemp); - -/* -------------------------------------------------------------------- */ -/* Look for GeoTransform Array */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Enable GeoTransform */ + /* -------------------------------------------------------------------- */ + poDS->bGotGeoTransform = TRUE; - if( pszGeoTransform != NULL ) { - papszGeoTransform = CSLTokenizeString2( pszGeoTransform, - " ", - CSLT_HONOURSTRINGS ); - poDS->bGotGeoTransform = TRUE; - - poDS->adfGeoTransform[0] = atof( papszGeoTransform[0] ); - poDS->adfGeoTransform[1] = atof( papszGeoTransform[1] ); - poDS->adfGeoTransform[2] = atof( papszGeoTransform[2] ); - poDS->adfGeoTransform[3] = atof( papszGeoTransform[3] ); - poDS->adfGeoTransform[4] = atof( papszGeoTransform[4] ); - poDS->adfGeoTransform[5] = atof( papszGeoTransform[5] ); -/* -------------------------------------------------------------------- */ -/* Look for corner array values */ -/* -------------------------------------------------------------------- */ - } else { - double dfNN, dfSN, dfEE, dfWE; - strcpy(szTemp,szGridMappingValue); - strcat( szTemp, "#" ); - strcat( szTemp, "Northernmost_Northing"); - pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); - - if( pszValue != NULL ) { - dfNN = atof( pszValue ); - } - strcpy(szTemp,szGridMappingValue); - strcat( szTemp, "#" ); - strcat( szTemp, "Southernmost_Northing"); - pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); - - if( pszValue != NULL ) { - dfSN = atof( pszValue ); - } - - strcpy(szTemp,szGridMappingValue); - strcat( szTemp, "#" ); - strcat( szTemp, "Easternmost_Easting"); - pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); - - if( pszValue != NULL ) { - dfEE = atof( pszValue ); - } - - strcpy(szTemp,szGridMappingValue); - strcat( szTemp, "#" ); - strcat( szTemp, "Westernmost_Easting"); - pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); - - if( pszValue != NULL ) { - dfWE = atof( pszValue ); - } - - adfGeoTransform[0] = dfWE; - adfGeoTransform[1] = (dfEE - dfWE) / - ( poDS->GetRasterXSize() - 1 ); - adfGeoTransform[2] = 0.0; - adfGeoTransform[3] = dfNN; - adfGeoTransform[4] = 0.0; - adfGeoTransform[5] = (dfSN - dfNN) / - ( poDS->GetRasterYSize() - 1 ); -/* -------------------------------------------------------------------- */ -/* Compute the center of the pixel */ -/* -------------------------------------------------------------------- */ - adfGeoTransform[0] = dfWE - - (adfGeoTransform[1] / 2); + poDS->adfGeoTransform[0] = pdfXCoord[0]; + poDS->adfGeoTransform[3] = pdfYCoord[0]; + poDS->adfGeoTransform[2] = 0; + poDS->adfGeoTransform[4] = 0; + poDS->adfGeoTransform[1] = (( pdfXCoord[xdim-1] - + pdfXCoord[0] ) / + ( poDS->nRasterXSize - 1 )); - adfGeoTransform[3] = dfNN - - (adfGeoTransform[5] / 2); + poDS->adfGeoTransform[5] = (( pdfYCoord[ydim-1] - + pdfYCoord[0] ) / + ( poDS->nRasterYSize - 1 )); + /* -------------------------------------------------------------------- */ + /* Compute the center of the pixel */ + /* -------------------------------------------------------------------- */ + poDS->adfGeoTransform[0] = pdfXCoord[0] + - (poDS->adfGeoTransform[1] / 2); + poDS->adfGeoTransform[3] = pdfYCoord[0] + - (poDS->adfGeoTransform[5] / 2); - bGotGeoTransform = TRUE; - } - CSLDestroy( papszGeoTransform ); - } - } + oSRS.exportToWkt( &(poDS->pszProjection) ); + } + } + CPLFree( pdfXCoord ); + CPLFree( pdfYCoord ); + + } + + + /* -------------------------------------------------------------------- */ + /* Is this a netCDF file created by GDAL? */ + /* -------------------------------------------------------------------- */ + if( !EQUAL( szGridMappingValue, "" ) ) { + strcpy( szTemp,szGridMappingValue); + strcat( szTemp, "#" ); + strcat( szTemp, "spatial_ref"); + pszWKT = CSLFetchNameValue(poDS->papszMetadata, szTemp); + + if( pszWKT != NULL ) { + + pszProjection = strdup( pszWKT ); + + strcpy(szTemp,szGridMappingValue); + strcat( szTemp, "#" ); + strcat( szTemp, "GeoTransform"); + + pszGeoTransform = CSLFetchNameValue(poDS->papszMetadata, szTemp); + + /* -------------------------------------------------------------------- */ + /* Look for GeoTransform Array */ + /* -------------------------------------------------------------------- */ + + if( pszGeoTransform != NULL ) { + papszGeoTransform = CSLTokenizeString2( pszGeoTransform, + " ", + CSLT_HONOURSTRINGS ); + poDS->bGotGeoTransform = TRUE; + + poDS->adfGeoTransform[0] = atof( papszGeoTransform[0] ); + poDS->adfGeoTransform[1] = atof( papszGeoTransform[1] ); + poDS->adfGeoTransform[2] = atof( papszGeoTransform[2] ); + poDS->adfGeoTransform[3] = atof( papszGeoTransform[3] ); + poDS->adfGeoTransform[4] = atof( papszGeoTransform[4] ); + poDS->adfGeoTransform[5] = atof( papszGeoTransform[5] ); + /* -------------------------------------------------------------------- */ + /* Look for corner array values */ + /* -------------------------------------------------------------------- */ + } else { + double dfNN, dfSN, dfEE, dfWE; + strcpy(szTemp,szGridMappingValue); + strcat( szTemp, "#" ); + strcat( szTemp, "Northernmost_Northing"); + pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); + + if( pszValue != NULL ) { + dfNN = atof( pszValue ); + } + strcpy(szTemp,szGridMappingValue); + strcat( szTemp, "#" ); + strcat( szTemp, "Southernmost_Northing"); + pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); + + if( pszValue != NULL ) { + dfSN = atof( pszValue ); + } + + strcpy(szTemp,szGridMappingValue); + strcat( szTemp, "#" ); + strcat( szTemp, "Easternmost_Easting"); + pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); + + if( pszValue != NULL ) { + dfEE = atof( pszValue ); + } + + strcpy(szTemp,szGridMappingValue); + strcat( szTemp, "#" ); + strcat( szTemp, "Westernmost_Easting"); + pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); + + if( pszValue != NULL ) { + dfWE = atof( pszValue ); + } + + adfGeoTransform[0] = dfWE; + adfGeoTransform[1] = (dfEE - dfWE) / + ( poDS->GetRasterXSize() - 1 ); + adfGeoTransform[2] = 0.0; + adfGeoTransform[3] = dfNN; + adfGeoTransform[4] = 0.0; + adfGeoTransform[5] = (dfSN - dfNN) / + ( poDS->GetRasterYSize() - 1 ); + /* -------------------------------------------------------------------- */ + /* Compute the center of the pixel */ + /* -------------------------------------------------------------------- */ + adfGeoTransform[0] = dfWE + - (adfGeoTransform[1] / 2); + + adfGeoTransform[3] = dfNN + - (adfGeoTransform[5] / 2); + + + bGotGeoTransform = TRUE; + } + CSLDestroy( papszGeoTransform ); + + } + } + } /************************************************************************/ /* GetGeoTransform() */ @@ -1085,93 +1426,93 @@ nc_inq_varnatts( cdfid, var, &nbAttr ); if( var == NC_GLOBAL ) { - strcpy( szVarName,"NC_GLOBAL" ); + strcpy( szVarName,"NC_GLOBAL" ); } else { - nc_inq_varname( cdfid, var, szVarName ); + nc_inq_varname( cdfid, var, szVarName ); } for( int l=0; l < nbAttr; l++){ - - nc_inq_attname( cdfid, var, l, szAttrName); - sprintf( szMetaName, "%s#%s", szVarName, szAttrName ); - *szMetaTemp='\0'; - nc_inq_att( cdfid, var, szAttrName, &nAttrType, &nAttrLen ); - - - switch (nAttrType) { - case NC_CHAR: - char *pszTemp; - pszTemp = (char *) CPLCalloc( nAttrLen+1, sizeof( char ) ); - nc_get_att_text( cdfid, var, szAttrName,pszTemp ); - pszTemp[nAttrLen]='\0'; - strcpy(szMetaTemp,pszTemp); - CPLFree(pszTemp); - break; - case NC_SHORT: - short *psTemp; - - psTemp = (short *) CPLCalloc( nAttrLen, sizeof( short ) ); - nc_get_att_short( cdfid, var, szAttrName, psTemp ); - for(m=0; m < nAttrLen-1; m++) { - sprintf( szTemp, "%d, ",psTemp[m] ); - strcat(szMetaTemp,szTemp); - } - sprintf( szTemp, "%d",psTemp[m] ); - CPLFree(psTemp); - strcat(szMetaTemp,szTemp); - - break; - case NC_INT: - int *pnTemp; - - pnTemp = (int *) CPLCalloc( nAttrLen, sizeof( int ) ); - nc_get_att_int( cdfid, var, szAttrName, pnTemp ); - for(m=0; m < nAttrLen-1; m++) { - sprintf( szTemp, "%d",pnTemp[m] ); - strcat(szMetaTemp,szTemp); - } - sprintf( szTemp, "%d",pnTemp[m] ); - CPLFree(pnTemp); - strcat(szMetaTemp,szTemp); - break; - case NC_FLOAT: - float *pfTemp; - pfTemp = (float *) CPLCalloc( nAttrLen, sizeof( float ) ); - nc_get_att_float( cdfid, var, szAttrName, pfTemp ); - for(m=0; m < nAttrLen-1; m++) { - sprintf( szTemp, "%e",pfTemp[m] ); - strcat(szMetaTemp,szTemp); - } - sprintf( szTemp, "%e",pfTemp[m] ); - CPLFree(pfTemp); - strcat(szMetaTemp,szTemp); - - break; - case NC_DOUBLE: - double *pdfTemp; - pdfTemp = (double *) CPLCalloc(nAttrLen, sizeof(double)); - nc_get_att_double( cdfid, var, szAttrName, pdfTemp ); - for(m=0; m < nAttrLen-1; m++) { - sprintf( szTemp, "%g",pdfTemp[m] ); - strcat(szMetaTemp,szTemp); - } - sprintf( szTemp, "%g",pdfTemp[m] ); - CPLFree(pdfTemp); - strcat(szMetaTemp,szTemp); - - break; - default: - break; - } - papszMetadata = CSLSetNameValue(papszMetadata, - szMetaName, - szMetaTemp); - + nc_inq_attname( cdfid, var, l, szAttrName); + sprintf( szMetaName, "%s#%s", szVarName, szAttrName ); + *szMetaTemp='\0'; + nc_inq_att( cdfid, var, szAttrName, &nAttrType, &nAttrLen ); + + + switch (nAttrType) { + case NC_CHAR: + char *pszTemp; + pszTemp = (char *) CPLCalloc( nAttrLen+1, sizeof( char ) ); + nc_get_att_text( cdfid, var, szAttrName,pszTemp ); + pszTemp[nAttrLen]='\0'; + strcpy(szMetaTemp,pszTemp); + CPLFree(pszTemp); + break; + case NC_SHORT: + short *psTemp; + + psTemp = (short *) CPLCalloc( nAttrLen, sizeof( short ) ); + nc_get_att_short( cdfid, var, szAttrName, psTemp ); + for(m=0; m < nAttrLen-1; m++) { + sprintf( szTemp, "%d, ",psTemp[m] ); + strcat(szMetaTemp,szTemp); + } + sprintf( szTemp, "%d",psTemp[m] ); + CPLFree(psTemp); + strcat(szMetaTemp,szTemp); + + break; + case NC_INT: + int *pnTemp; + + pnTemp = (int *) CPLCalloc( nAttrLen, sizeof( int ) ); + nc_get_att_int( cdfid, var, szAttrName, pnTemp ); + for(m=0; m < nAttrLen-1; m++) { + sprintf( szTemp, "%d",pnTemp[m] ); + strcat(szMetaTemp,szTemp); + } + sprintf( szTemp, "%d",pnTemp[m] ); + CPLFree(pnTemp); + strcat(szMetaTemp,szTemp); + break; + case NC_FLOAT: + float *pfTemp; + pfTemp = (float *) CPLCalloc( nAttrLen, sizeof( float ) ); + nc_get_att_float( cdfid, var, szAttrName, pfTemp ); + for(m=0; m < nAttrLen-1; m++) { + sprintf( szTemp, "%e",pfTemp[m] ); + strcat(szMetaTemp,szTemp); + } + sprintf( szTemp, "%e",pfTemp[m] ); + CPLFree(pfTemp); + strcat(szMetaTemp,szTemp); + + break; + case NC_DOUBLE: + double *pdfTemp; + pdfTemp = (double *) CPLCalloc(nAttrLen, sizeof(double)); + nc_get_att_double( cdfid, var, szAttrName, pdfTemp ); + for(m=0; m < nAttrLen-1; m++) { + sprintf( szTemp, "%g",pdfTemp[m] ); + strcat(szMetaTemp,szTemp); + } + sprintf( szTemp, "%g",pdfTemp[m] ); + CPLFree(pdfTemp); + strcat(szMetaTemp,szTemp); + + break; + default: + break; + } + + papszMetadata = CSLSetNameValue(papszMetadata, + szMetaName, + szMetaTemp); + } - + return CE_None; } @@ -1205,79 +1546,79 @@ nc_inq_nvars ( cdfid, &nVarCount ); for ( nVar = 0; nVar < nVarCount; nVar++ ) { - nc_inq_varndims ( cdfid, nVar, &nDims ); - if( nDims >= 2 ) { - ponDimIds = (int *) CPLCalloc( nDims, sizeof( int ) ); - nc_inq_vardimid ( cdfid, nVar, ponDimIds ); - -/* -------------------------------------------------------------------- */ -/* Create Sub dataset list */ -/* -------------------------------------------------------------------- */ - szDim[0]='\0'; - for( i = 0; i < nDims; i++ ) { - nc_inq_dimlen ( cdfid, ponDimIds[i], &nDimLen ); - sprintf(szTemp, "%d", (int) nDimLen); - strcat(szTemp, "x" ); - strcat(szDim, szTemp); - } + nc_inq_varndims ( cdfid, nVar, &nDims ); + if( nDims >= 2 ) { + ponDimIds = (int *) CPLCalloc( nDims, sizeof( int ) ); + nc_inq_vardimid ( cdfid, nVar, ponDimIds ); - nc_inq_vartype( cdfid, nVar, &nVarType ); -/* -------------------------------------------------------------------- */ -/* Get rid of the last "x" character */ -/* -------------------------------------------------------------------- */ - szDim[strlen(szDim) - 1] = '\0'; - - switch( nVarType ) { - case NC_BYTE: - strcpy(szType, "8-bit character"); - break; - case NC_SHORT: - strcpy(szType, "8-bit integer"); - break; - case NC_INT: - strcpy(szType, "16-bit integer"); - break; - case NC_FLOAT: - strcpy(szType, "32-bit floating-point"); - break; - case NC_DOUBLE: - strcpy(szType, "64-bit floating-point"); - break; - default: - break; - } - nc_inq_varname ( cdfid, nVar, szName); - nc_inq_att( cdfid, nVar, "standard_name", &nAttype, &nAttlen); - if( nc_get_att_text ( cdfid, nVar, "standard_name", - szVarStdName ) == NC_NOERR ) { - szVarStdName[nAttlen] = '\0'; - } - else { - strcpy( szVarStdName, szName ); - } - - sprintf( szTemp, "SUBDATASET_%d_NAME", nSub) ; - - poDS->papszSubDatasets = - CSLSetNameValue( poDS->papszSubDatasets, szTemp, - CPLSPrintf( "NETCDF:\"%s\":%s", - poDS->pszFilename, - szName) ) ; - - sprintf( szTemp, "SUBDATASET_%d_DESC", nSub++ ); + /* -------------------------------------------------------------------- */ + /* Create Sub dataset list */ + /* -------------------------------------------------------------------- */ + szDim[0]='\0'; + for( i = 0; i < nDims; i++ ) { + nc_inq_dimlen ( cdfid, ponDimIds[i], &nDimLen ); + sprintf(szTemp, "%d", (int) nDimLen); + strcat(szTemp, "x" ); + strcat(szDim, szTemp); + } - poDS->papszSubDatasets = - CSLSetNameValue( poDS->papszSubDatasets, szTemp, - CPLSPrintf( "[%s] %s (%s)", - szDim, - szVarStdName, - szType ) ); - CPLFree(ponDimIds); - } + nc_inq_vartype( cdfid, nVar, &nVarType ); + /* -------------------------------------------------------------------- */ + /* Get rid of the last "x" character */ + /* -------------------------------------------------------------------- */ + szDim[strlen(szDim) - 1] = '\0'; + + switch( nVarType ) { + case NC_BYTE: + strcpy(szType, "8-bit character"); + break; + case NC_SHORT: + strcpy(szType, "8-bit integer"); + break; + case NC_INT: + strcpy(szType, "16-bit integer"); + break; + case NC_FLOAT: + strcpy(szType, "32-bit floating-point"); + break; + case NC_DOUBLE: + strcpy(szType, "64-bit floating-point"); + break; + default: + break; + } + nc_inq_varname ( cdfid, nVar, szName); + nc_inq_att( cdfid, nVar, "standard_name", &nAttype, &nAttlen); + if( nc_get_att_text ( cdfid, nVar, "standard_name", + szVarStdName ) == NC_NOERR ) { + szVarStdName[nAttlen] = '\0'; + } + else { + strcpy( szVarStdName, szName ); + } + + sprintf( szTemp, "SUBDATASET_%d_NAME", nSub) ; + + poDS->papszSubDatasets = + CSLSetNameValue( poDS->papszSubDatasets, szTemp, + CPLSPrintf( "NETCDF:\"%s\":%s", + poDS->pszFilename, + szName) ) ; + + sprintf( szTemp, "SUBDATASET_%d_DESC", nSub++ ); + + poDS->papszSubDatasets = + CSLSetNameValue( poDS->papszSubDatasets, szTemp, + CPLSPrintf( "[%s] %s (%s)", + szDim, + szVarStdName, + szType ) ); + CPLFree(ponDimIds); + } } } - + /************************************************************************/ /* Open() */ /************************************************************************/ @@ -1297,73 +1638,70 @@ int nDimID; char attname[NC_MAX_NAME]; int ndims, nvars, ngatts, unlimdimid; - int nCount=0; - int nVarID; -/* -------------------------------------------------------------------- */ -/* Does this appear to be a netcdf file? */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Does this appear to be a netcdf file? */ + /* -------------------------------------------------------------------- */ if( !EQUALN(poOpenInfo->pszFilename,"NETCDF:",7) - && ( poOpenInfo->nHeaderBytes < 5 - || !EQUALN((const char *) (poOpenInfo->pabyHeader),"CDF\001",5))) + && ( poOpenInfo->nHeaderBytes < 5 + || !EQUALN((const char *) (poOpenInfo->pabyHeader),"CDF\001",5))) return NULL; -/* -------------------------------------------------------------------- */ -/* Check if filename start with NETCDF: tag */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Check if filename start with NETCDF: tag */ + /* -------------------------------------------------------------------- */ netCDFDataset *poDS; poDS = new netCDFDataset(); - poDS->papszName = - CSLTokenizeString2( poOpenInfo->pszFilename, - ":", CSLT_HONOURSTRINGS|CSLT_PRESERVEESCAPES ); + poDS->papszName = CSLTokenizeString2( poOpenInfo->pszFilename, + ":", CSLT_HONOURSTRINGS ); if( EQUAL( poDS->papszName[0], "NETCDF" ) ) { - if( ( CSLCount(poDS->papszName) == 3 ) ){ - poDS->pszFilename = strdup( poDS->papszName[1] ); - } + if( ( CSLCount(poDS->papszName) == 3 ) ){ + poDS->pszFilename = strdup( poDS->papszName[1] ); + } } else { - poDS->pszFilename = strdup( poOpenInfo->pszFilename ); + poDS->pszFilename = strdup( poOpenInfo->pszFilename ); } -/* -------------------------------------------------------------------- */ -/* Try opening the dataset. */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Try opening the dataset. */ + /* -------------------------------------------------------------------- */ if( nc_open( poDS->pszFilename, NC_NOWRITE, &cdfid ) != NC_NOERR ) { - delete poDS; + delete poDS; return NULL; } -/* -------------------------------------------------------------------- */ -/* Is this a real netCDF file? */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Is this a real netCDF file? */ + /* -------------------------------------------------------------------- */ status = nc_inq(cdfid, &ndims, &nvars, &ngatts, &unlimdimid); if( status != NC_NOERR ) { - CPLFree( poDS->pszFilename ); - CSLDestroy( poDS->papszName ); + CPLFree( poDS->pszFilename ); + CSLDestroy( poDS->papszName ); return NULL; } -/* -------------------------------------------------------------------- */ -/* Does the request variable exist? */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Does the request variable exist? */ + /* -------------------------------------------------------------------- */ if( ( CSLCount(poDS->papszName) == 3 ) ){ - status = nc_inq_varid( cdfid, poDS->papszName[2], &var); - if( status != NC_NOERR ) { - CPLError( CE_Warning, CPLE_AppDefined, - "%s is a netCDF file, but %s is not a variable.", - poOpenInfo->pszFilename, - poDS->papszName[2] ); - - nc_close( cdfid ); - return NULL; - } + status = nc_inq_varid( cdfid, poDS->papszName[2], &var); + if( status != NC_NOERR ) { + CPLError( CE_Warning, CPLE_AppDefined, + "%s is a netCDF file, but %s is not a variable.", + poOpenInfo->pszFilename, + poDS->papszName[2] ); + + nc_close( cdfid ); + return NULL; + } } if( nc_inq_ndims( cdfid, &dim_count ) != NC_NOERR || dim_count < 2 ) { - CPLError( CE_Warning, CPLE_AppDefined, + CPLError( CE_Warning, CPLE_AppDefined, "%s is a netCDF file, but not in GMT configuration.", poOpenInfo->pszFilename ); @@ -1374,212 +1712,186 @@ CPLDebug( "GDAL_netCDF", "dim_count = %d\n", dim_count ); if( (status = nc_get_att_text( cdfid, NC_GLOBAL, "Conventions", - attname )) != NC_NOERR ) { - CPLError( CE_Warning, CPLE_AppDefined, - "No UNIDATA NC_GLOBAL:Conventions attribute"); - /* note that 'Conventions' is always capital 'C' in CF spec*/ - } + attname )) != NC_NOERR ) { + CPLError( CE_Warning, CPLE_AppDefined, + "No UNIDATA NC_GLOBAL:Conventions attribute"); + /* note that 'Conventions' is always capital 'C' in CF spec*/ + } -/* -------------------------------------------------------------------- */ -/* Create band information objects. */ -/* -------------------------------------------------------------------- */ - if ( nc_inq_nvars ( cdfid, &var_count) != NC_NOERR ) - return NULL; - - CPLDebug( "GDAL_netCDF", "var_count = %d\n", var_count ); - + /* -------------------------------------------------------------------- */ + /* Create band information objects. */ + /* -------------------------------------------------------------------- */ + if ( nc_inq_nvars ( cdfid, &var_count) != NC_NOERR ) + return NULL; -/* -------------------------------------------------------------------- */ -/* Create a corresponding GDALDataset. */ -/* Create Netcdf Subdataset if filename as NETCDF tag */ -/* -------------------------------------------------------------------- */ - poDS->cdfid = cdfid; + CPLDebug( "GDAL_netCDF", "var_count = %d\n", var_count ); - poDS->ReadAttributes( cdfid, NC_GLOBAL ); -/* -------------------------------------------------------------------- */ -/* Verify if only one variable has 2 dimensions */ -/* -------------------------------------------------------------------- */ - for ( j = 0; j < nvars; j++ ) { + /* -------------------------------------------------------------------- */ + /* Create a corresponding GDALDataset. */ + /* Create Netcdf Subdataset if filename as NETCDF tag */ + /* -------------------------------------------------------------------- */ + poDS->cdfid = cdfid; - nc_inq_varndims ( cdfid, j, &ndims ); - if( ndims >= 2 ) { - nVarID=j; - nCount++; - } - } + poDS->ReadAttributes( cdfid, NC_GLOBAL ); -/* -------------------------------------------------------------------- */ -/* We have more than one variable with 2 dimensions in the file */ -/* -------------------------------------------------------------------- */ - if( (nCount > 1) && ( !EQUAL( poDS->papszName[0], "NETCDF" ) ) ) { - poDS->CreateSubDatasetList( ); - poDS->SetMetadata( poDS->papszMetadata ); - return( poDS ); - } -/* -------------------------------------------------------------------- */ -/* If we have only one varialbe */ -/* Generate a new filename with format NETCDF:"filename":subdataset */ -/* -------------------------------------------------------------------- */ - if( ( nCount == 1 ) && ( !EQUAL( poDS->papszName[0], "NETCDF" ) ) ){ - CPLFree( poDS->pszFilename ); - CSLDestroy( poDS->papszName ); + if( ( !EQUAL( poDS->papszName[0], "NETCDF" ) ) ) { + poDS->CreateSubDatasetList( ); + poDS->SetMetadata( poDS->papszMetadata ); + return( poDS ); + } - char pszNETCDFFilename[NC_MAX_NAME]; - char szVarName[NC_MAX_NAME]; - nc_inq_varname( cdfid, nVarID, szVarName); + /* -------------------------------------------------------------------- */ + /* Open the NETCDF subdataset NETCDF:"filename":subdataset */ + /* -------------------------------------------------------------------- */ - strcpy( pszNETCDFFilename,"NETCDF:\"" ); - strcat( pszNETCDFFilename,poOpenInfo->pszFilename ); - strcat( pszNETCDFFilename,"\":" ); - strcat( pszNETCDFFilename, szVarName ); + var=-1; + nc_inq_varid( cdfid, poDS->papszName[2], &var); + nc_inq_varndims ( cdfid, var, &nd ); - CPLDebug( "GDAL_netCDF", "NETCDFFilename = %s\n", - pszNETCDFFilename); + poDS->paDimIds = (int *)CPLCalloc(nd, sizeof( int ) ); + poDS->panBandDimPos = ( int * ) CPLCalloc( nd, sizeof( int ) ); - poDS->papszName = CSLTokenizeString2( pszNETCDFFilename, - ":", CSLT_HONOURSTRINGS ); - poDS->pszFilename = strdup( pszNETCDFFilename ); + nc_inq_vardimid( cdfid, var, poDS->paDimIds ); - } + /* -------------------------------------------------------------------- */ + /* Check fi somebody tried to pass a variable with less than 2D */ + /* -------------------------------------------------------------------- */ -/* -------------------------------------------------------------------- */ -/* Open the NETCDF subdataset NETCDF:"filename":subdataset */ -/* -------------------------------------------------------------------- */ + if ( nd < 2 ) { + CPLFree( poDS->paDimIds ); + CPLFree( poDS->panBandDimPos ); + return NULL; + } - var=-1; - nc_inq_varid( cdfid, poDS->papszName[2], &var); - nc_inq_varndims ( cdfid, var, &nd ); + /* -------------------------------------------------------------------- */ + /* CF-1 Convention */ + /* dimensions to appear in the relative order T, then Z, then Y, */ + /* then X to the file. All other dimensions should, whenever */ + /* possible, be placed to the left of the spatiotemporal */ + /* dimensions. */ + /* -------------------------------------------------------------------- */ - poDS->paDimIds = (int *)CPLCalloc(nd, sizeof( int ) ); - poDS->panBandDimPos = ( int * ) CPLCalloc( nd, sizeof( int ) ); + /* -------------------------------------------------------------------- */ + /* Get X dimensions information */ + /* -------------------------------------------------------------------- */ + poDS->nDimXid = poDS->paDimIds[nd-1]; + nc_inq_dimlen ( cdfid, poDS->nDimXid, &poDS->xdim ); + poDS->nRasterXSize = poDS->xdim; - nc_inq_vardimid( cdfid, var, poDS->paDimIds ); + /* -------------------------------------------------------------------- */ + /* Get Y dimension information */ + /* -------------------------------------------------------------------- */ + poDS->nDimYid = poDS->paDimIds[nd-2]; + nc_inq_dimlen ( cdfid, poDS->nDimYid, &poDS->ydim ); + poDS->nRasterYSize = poDS->ydim; - -/* -------------------------------------------------------------------- */ -/* Check fi somebody tried to pass a variable with less than 2D */ -/* -------------------------------------------------------------------- */ - if ( nd < 2 ) { - CPLFree( poDS->paDimIds ); - CPLFree( poDS->panBandDimPos ); - return NULL; - } + for( j=0,k=0; j < nd; j++ ){ + if( poDS->paDimIds[j] == poDS->nDimXid ){ + poDS->panBandDimPos[0] = j; // Save Position of XDim + k++; + } + if( poDS->paDimIds[j] == poDS->nDimYid ){ + poDS->panBandDimPos[1] = j; // Save Position of YDim + k++; + } + } + /* -------------------------------------------------------------------- */ + /* X and Y Dimension Ids were not found! */ + /* -------------------------------------------------------------------- */ + if( k != 2 ) { + CPLFree( poDS->paDimIds ); + CPLFree( poDS->panBandDimPos ); + return NULL; + } -/* -------------------------------------------------------------------- */ -/* CF-1 Convention */ -/* dimensions to appear in the relative order T, then Z, then Y, */ -/* then X to the file. All other dimensions should, whenever */ -/* possible, be placed to the left of the spatiotemporal */ -/* dimensions. */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Read Metadata for this variable */ + /* -------------------------------------------------------------------- */ + poDS->ReadAttributes( cdfid, var ); -/* -------------------------------------------------------------------- */ -/* Get X dimensions information */ -/* -------------------------------------------------------------------- */ - poDS->nDimXid = poDS->paDimIds[nd-1]; - nc_inq_dimlen ( cdfid, poDS->nDimXid, &poDS->xdim ); - poDS->nRasterXSize = poDS->xdim; + /* -------------------------------------------------------------------- */ + /* Read Metadata for each dimension */ + /* -------------------------------------------------------------------- */ -/* -------------------------------------------------------------------- */ -/* Get Y dimension information */ -/* -------------------------------------------------------------------- */ - poDS->nDimYid = poDS->paDimIds[nd-2]; - nc_inq_dimlen ( cdfid, poDS->nDimYid, &poDS->ydim ); - poDS->nRasterYSize = poDS->ydim; + for( j=0; j < dim_count; j++ ){ + nc_inq_dimname( cdfid, j, poDS->papszDimName[j] ); + status = nc_inq_varid( cdfid, poDS->papszDimName[j], &nDimID ); + if( status == NC_NOERR ) { + poDS->ReadAttributes( cdfid, nDimID ); + } + } + poDS->SetProjection( var ); + /* -------------------------------------------------------------------- */ + /* Create bands */ + /* -------------------------------------------------------------------- */ + poDS->panBandZLev = (int *)CPLCalloc( sizeof( nd ) - 2, + sizeof( int ) ); - for( j=0,k=0; j < nd; j++ ){ - if( poDS->paDimIds[j] == poDS->nDimXid ){ - poDS->panBandDimPos[0] = j; // Save Position of XDim - k++; - } - if( poDS->paDimIds[j] == poDS->nDimYid ){ - poDS->panBandDimPos[1] = j; // Save Position of YDim - k++; - } - } -/* -------------------------------------------------------------------- */ -/* X and Y Dimension Ids were not found! */ -/* -------------------------------------------------------------------- */ - if( k != 2 ) { - CPLFree( poDS->paDimIds ); - CPLFree( poDS->panBandDimPos ); - return NULL; - } - -/* -------------------------------------------------------------------- */ -/* Read Metadata for this variable */ -/* -------------------------------------------------------------------- */ - poDS->ReadAttributes( cdfid, var ); - -/* -------------------------------------------------------------------- */ -/* Read Metadata for each dimension */ -/* -------------------------------------------------------------------- */ - - for( j=0; j < dim_count; j++ ){ - nc_inq_dimname( cdfid, j, poDS->papszDimName[j] ); - status = nc_inq_varid( cdfid, poDS->papszDimName[j], &nDimID ); - if( status == NC_NOERR ) { - poDS->ReadAttributes( cdfid, nDimID ); - } - } + nTotLevCount = 1; + if ( dim_count > 2 ) { + nDim=2; + for( j=0; j < nd; j++ ){ + if( ( poDS->paDimIds[j] != poDS->nDimXid ) && + ( poDS->paDimIds[j] != poDS->nDimYid ) ){ + nc_inq_dimlen ( cdfid, poDS->paDimIds[j], &lev_count ); + nTotLevCount *= lev_count; + poDS->panBandZLev[ nDim-2 ] = lev_count; + poDS->panBandDimPos[ nDim++ ] = j; //Save Position of ZDim + } + } + } + i=0; - poDS->SetProjection( var ); -/* -------------------------------------------------------------------- */ -/* Create bands */ -/* -------------------------------------------------------------------- */ - poDS->panBandZLev = (int *)CPLCalloc( sizeof( nd ) - 2, - sizeof( int ) ); - - nTotLevCount = 1; - if ( dim_count > 2 ) { - nDim=2; - for( j=0; j < nd; j++ ){ - if( ( poDS->paDimIds[j] != poDS->nDimXid ) && - ( poDS->paDimIds[j] != poDS->nDimYid ) ){ - nc_inq_dimlen ( cdfid, poDS->paDimIds[j], &lev_count ); - nTotLevCount *= lev_count; - poDS->panBandZLev[ nDim-2 ] = lev_count; - poDS->panBandDimPos[ nDim++ ] = j; //Save Position of ZDim - } - } - } - i=0; + for ( unsigned int lev = 0; lev < nTotLevCount ; lev++ ) { + char szBandItemName[ MAX_STR_LEN ]; + char ** papszToken; + const char * pszGeo; + papszToken=NULL; - for ( unsigned int lev = 0; lev < nTotLevCount ; lev++ ) { - char ** papszToken; - papszToken=NULL; + netCDFRasterBand *poBand=new netCDFRasterBand(poDS, + var, + nDim, + lev, + poDS->panBandZLev, + poDS->panBandDimPos, + i+1 ); - netCDFRasterBand *poBand = - new netCDFRasterBand(poDS, var, nDim, lev, - poDS->panBandZLev, poDS->panBandDimPos, i+1 ); + for( j = 0; j < CSLCount( poDS->papszGeolocation ); j++ ) { + pszGeo = CSLGetField( poDS->papszGeolocation, j ); + papszToken = CSLTokenizeString2( pszGeo, + "=", CSLT_HONOURSTRINGS ); - poDS->SetBand( i+1, poBand ); - i++; - } + strcpy( szBandItemName, "NETCDF_GEOLOCATION_" ); + strcat( szBandItemName, papszToken[0] ); + poBand->SetMetadataItem( szBandItemName, papszToken[1] ); + CSLDestroy( papszToken ); + } - CPLFree( poDS->paDimIds ); - CPLFree( poDS->panBandDimPos ); - CPLFree( poDS->panBandZLev ); - - poDS->nBands = i; + poDS->SetBand( i+1, poBand ); + i++; + } - // Handle angular geographic coordinates here + CPLFree( poDS->paDimIds ); + CPLFree( poDS->panBandDimPos ); + CPLFree( poDS->panBandZLev ); -/* -------------------------------------------------------------------- */ -/* Initialize any PAM information. */ -/* -------------------------------------------------------------------- */ - poDS->SetMetadata( poDS->papszMetadata ); + poDS->nBands = i; - poDS->SetDescription( poOpenInfo->pszFilename ); - poDS->TryLoadXML( ); + // Handle angular geographic coordinates here - return( poDS ); + /* -------------------------------------------------------------------- */ + /* Initialize any PAM information. */ + /* -------------------------------------------------------------------- */ + poDS->SetDescription( poOpenInfo->pszFilename ); + poDS->TryLoadXML( ); + poDS->SetMetadata( poDS->papszMetadata ); + return( poDS ); } @@ -1591,9 +1903,9 @@ void CopyMetadata( void *poDS, int fpImage, int CDFVarID ) { -/* -------------------------------------------------------------------- */ -/* Add CF-1.0 Conventions Global attribute */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Add CF-1.0 Conventions Global attribute */ + /* -------------------------------------------------------------------- */ char **papszMetadata; char **papszFieldData; const char *pszField; @@ -1602,54 +1914,54 @@ int nDataLength; int nItems; -/* -------------------------------------------------------------------- */ -/* Global metadata are set with NC_GLOBAL as the varid */ -/* -------------------------------------------------------------------- */ - + /* -------------------------------------------------------------------- */ + /* Global metadata are set with NC_GLOBAL as the varid */ + /* -------------------------------------------------------------------- */ + if( CDFVarID == NC_GLOBAL ) { - papszMetadata = GDALGetMetadata( (GDALDataset *) poDS,""); + papszMetadata = GDALGetMetadata( (GDALDataset *) poDS,""); - nc_put_att_text( fpImage, - NC_GLOBAL, - "Conventions", - 6, - "CF-1.0" ); + nc_put_att_text( fpImage, + NC_GLOBAL, + "Conventions", + 6, + "CF-1.0" ); } else { - papszMetadata = GDALGetMetadata( (GDALRasterBandH) poDS, NULL ); + papszMetadata = GDALGetMetadata( (GDALRasterBandH) poDS, NULL ); } - nItems = CSLCount( papszMetadata ); - + nItems = CSLCount( papszMetadata ); + for(int k=0; k < nItems; k++ ) { - pszField = CSLGetField( papszMetadata, k ); - papszFieldData = CSLTokenizeString2 (pszField, "=", - CSLT_HONOURSTRINGS ); - if( papszFieldData[1] != NULL ) { - strcpy( szMetaName, papszFieldData[ 0 ] ); - strcpy( szMetaValue, papszFieldData[ 1 ] ); + pszField = CSLGetField( papszMetadata, k ); + papszFieldData = CSLTokenizeString2 (pszField, "=", + CSLT_HONOURSTRINGS ); + if( papszFieldData[1] != NULL ) { + strcpy( szMetaName, papszFieldData[ 0 ] ); + strcpy( szMetaValue, papszFieldData[ 1 ] ); -/* -------------------------------------------------------------------- */ -/* netCDF attributes do not like the '#' character. */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* netCDF attributes do not like the '#' character. */ + /* -------------------------------------------------------------------- */ - for( unsigned int h=0; h < strlen( szMetaName ) -1 ; h++ ) { - if( szMetaName[h] == '#' ) szMetaName[h] = '-'; - } - - nDataLength = strlen( szMetaValue ); - nc_put_att_text( fpImage, - CDFVarID, - szMetaName, - nDataLength, - szMetaValue ); + for( unsigned int h=0; h < strlen( szMetaName ) -1 ; h++ ) { + if( szMetaName[h] == '#' ) szMetaName[h] = '-'; + } - - } - CSLDestroy( papszFieldData ); + nDataLength = strlen( szMetaValue ); + nc_put_att_text( fpImage, + CDFVarID, + szMetaName, + nDataLength, + szMetaValue ); + + } + CSLDestroy( papszFieldData ); + } } /************************************************************************/ @@ -1658,9 +1970,9 @@ static GDALDataset* -NCDFCreateCopy( const char * pszFilename, GDALDataset *poSrcDS, - int bStrict, char ** papszOptions, - GDALProgressFunc pfnProgress, void * pProgressData ) + NCDFCreateCopy( const char * pszFilename, GDALDataset *poSrcDS, + int bStrict, char ** papszOptions, + GDALProgressFunc pfnProgress, void * pProgressData ) { int nBands = poSrcDS->GetRasterCount(); @@ -1681,9 +1993,9 @@ bProgressive = CSLFetchBoolean( papszOptions, "PROGRESSIVE", FALSE ); -/* -------------------------------------------------------------------- */ -/* Create the dataset. */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Create the dataset. */ + /* -------------------------------------------------------------------- */ int fpImage; int status; int nXDimID = 0; @@ -1693,8 +2005,8 @@ if( status != NC_NOERR ) { - CPLError( CE_Failure, CPLE_OpenFailed, - "Unable to create netCDF file %s.\n", + CPLError( CE_Failure, CPLE_OpenFailed, + "Unable to create netCDF file %s.\n", pszFilename ); return NULL; } @@ -1715,535 +2027,535 @@ CopyMetadata((void *) poSrcDS, fpImage, NC_GLOBAL ); -/* -------------------------------------------------------------------- */ -/* Set Projection for netCDF data CF-1 Convention */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Set Projection for netCDF data CF-1 Convention */ + /* -------------------------------------------------------------------- */ OGRSpatialReference oSRS; char *pszWKT = (char *) poSrcDS->GetProjectionRef(); if( pszWKT != NULL ) - oSRS.importFromWkt( &pszWKT ); + oSRS.importFromWkt( &pszWKT ); if( oSRS.IsGeographic() ) { - int status; - int i; - int NCDFVarID; + int status; + int i; + int NCDFVarID; - double dfNN=0.0; - double dfSN=0.0; - double dfEE=0.0; - double dfWE=0.0; - double adfGeoTransform[6]; - char szGeoTransform[ MAX_STR_LEN ]; - char szTemp[ MAX_STR_LEN ]; + double dfNN=0.0; + double dfSN=0.0; + double dfEE=0.0; + double dfWE=0.0; + double adfGeoTransform[6]; + char szGeoTransform[ MAX_STR_LEN ]; + char szTemp[ MAX_STR_LEN ]; -/* -------------------------------------------------------------------- */ -/* Copy GeoTransform array from source */ -/* -------------------------------------------------------------------- */ - poSrcDS->GetGeoTransform( adfGeoTransform ); + /* -------------------------------------------------------------------- */ + /* Copy GeoTransform array from source */ + /* -------------------------------------------------------------------- */ + poSrcDS->GetGeoTransform( adfGeoTransform ); bWriteGeoTransform = TRUE; - *szGeoTransform = '\0'; - for( i=0; i<6; i++ ) { - sprintf( szTemp, "%g ", - adfGeoTransform[i] ); - strcat( szGeoTransform, szTemp ); - } - CPLDebug( "GDAL_netCDF", "szGeoTranform = %s", szGeoTransform ); - strcpy( pszNetcdfProjection, "GDAL_Geographics" ); + *szGeoTransform = '\0'; + for( i=0; i<6; i++ ) { + sprintf( szTemp, "%g ", + adfGeoTransform[i] ); + strcat( szGeoTransform, szTemp ); + } + CPLDebug( "GDAL_netCDF", "szGeoTranform = %s", szGeoTransform ); + strcpy( pszNetcdfProjection, "GDAL_Geographics" ); - status = nc_def_var( fpImage, - pszNetcdfProjection, - NC_CHAR, - 0, NULL, &NCDFVarID ); - - dfNN = adfGeoTransform[3]; - dfSN = ( adfGeoTransform[5] * nYSize ) + dfNN; - dfWE = adfGeoTransform[0]; - dfEE = ( adfGeoTransform[1] * nXSize ) + dfWE; + status = nc_def_var( fpImage, + pszNetcdfProjection, + NC_CHAR, + 0, NULL, &NCDFVarID ); + dfNN = adfGeoTransform[3]; + dfSN = ( adfGeoTransform[5] * nYSize ) + dfNN; + dfWE = adfGeoTransform[0]; + dfEE = ( adfGeoTransform[1] * nXSize ) + dfWE; + status = nc_put_att_double( fpImage, - NCDFVarID, - "Northernmost_Northing", - NC_DOUBLE, - 1, - &dfNN ); + NCDFVarID, + "Northernmost_Northing", + NC_DOUBLE, + 1, + &dfNN ); status = nc_put_att_double( fpImage, - NCDFVarID, - "Southernmost_Northing", - NC_DOUBLE, - 1, - &dfSN ); - status = nc_put_att_double( fpImage, - NCDFVarID, - "Easternmost_Easting", - NC_DOUBLE, - 1, - &dfEE ); - status = nc_put_att_double( fpImage, - NCDFVarID, - "Westernmost_Easting", - NC_DOUBLE, - 1, - &dfWE ); - pszWKT = (char *) poSrcDS->GetProjectionRef() ; + NCDFVarID, + "Southernmost_Northing", + NC_DOUBLE, + 1, + &dfSN ); + status = nc_put_att_double( fpImage, + NCDFVarID, + "Easternmost_Easting", + NC_DOUBLE, + 1, + &dfEE ); + status = nc_put_att_double( fpImage, + NCDFVarID, + "Westernmost_Easting", + NC_DOUBLE, + 1, + &dfWE ); + pszWKT = (char *) poSrcDS->GetProjectionRef() ; - nc_put_att_text( fpImage, - NCDFVarID, - "spatial_ref", - strlen( pszWKT ), - pszWKT ); + nc_put_att_text( fpImage, + NCDFVarID, + "spatial_ref", + strlen( pszWKT ), + pszWKT ); - nc_put_att_text( fpImage, - NCDFVarID, - "GeoTransform", - strlen( szGeoTransform ), - szGeoTransform ); + nc_put_att_text( fpImage, + NCDFVarID, + "GeoTransform", + strlen( szGeoTransform ), + szGeoTransform ); - nc_put_att_text( fpImage, - NCDFVarID, - GRD_MAPPING_NAME, - 30, - "Geographics Coordinate System" ); - nc_put_att_text( fpImage, - NCDFVarID, - LNG_NAME, - 14, - "Grid_latitude" ); + nc_put_att_text( fpImage, + NCDFVarID, + GRD_MAPPING_NAME, + 30, + "Geographics Coordinate System" ); + nc_put_att_text( fpImage, + NCDFVarID, + LNG_NAME, + 14, + "Grid_latitude" ); } if( oSRS.IsProjected() ) { - const char *pszParamStr, *pszParamVal; + const char *pszParamStr, *pszParamVal; const OGR_SRSNode *poPROJCS = oSRS.GetAttrNode( "PROJCS" ); - int status; - int i; - int NCDFVarID; - const char *pszProjection; - double dfNN=0.0; - double dfSN=0.0; - double dfEE=0.0; - double dfWE=0.0; - double adfGeoTransform[6]; - char szGeoTransform[ MAX_STR_LEN ]; - char szTemp[ MAX_STR_LEN ]; + int status; + int i; + int NCDFVarID; + const char *pszProjection; + double dfNN=0.0; + double dfSN=0.0; + double dfEE=0.0; + double dfWE=0.0; + double adfGeoTransform[6]; + char szGeoTransform[ MAX_STR_LEN ]; + char szTemp[ MAX_STR_LEN ]; - poSrcDS->GetGeoTransform( adfGeoTransform ); + poSrcDS->GetGeoTransform( adfGeoTransform ); - *szGeoTransform = '\0'; - for( i=0; i<6; i++ ) { - sprintf( szTemp, "%g ", - adfGeoTransform[i] ); - strcat( szGeoTransform, szTemp ); - } + *szGeoTransform = '\0'; + for( i=0; i<6; i++ ) { + sprintf( szTemp, "%g ", + adfGeoTransform[i] ); + strcat( szGeoTransform, szTemp ); + } - CPLDebug( "GDAL_netCDF", "szGeoTranform = %s", szGeoTransform ); + CPLDebug( "GDAL_netCDF", "szGeoTranform = %s", szGeoTransform ); pszProjection = oSRS.GetAttrValue( "PROJECTION" ); bWriteGeoTransform = TRUE; - for(i=0; poNetcdfSRS[i].netCDFSRS != NULL; i++ ) { - if( EQUAL( poNetcdfSRS[i].SRS, pszProjection ) ) { - CPLDebug( "GDAL_netCDF", "PROJECTION = %s", - poNetcdfSRS[i].netCDFSRS); - strcpy( pszNetcdfProjection, poNetcdfSRS[i].netCDFSRS ); + for(i=0; poNetcdfSRS[i].netCDFSRS != NULL; i++ ) { + if( EQUAL( poNetcdfSRS[i].SRS, pszProjection ) ) { + CPLDebug( "GDAL_netCDF", "PROJECTION = %s", + poNetcdfSRS[i].netCDFSRS); + strcpy( pszNetcdfProjection, poNetcdfSRS[i].netCDFSRS ); - break; - } - } + break; + } + } - status = nc_def_var( fpImage, - poNetcdfSRS[i].netCDFSRS, - NC_CHAR, - 0, NULL, &NCDFVarID ); + status = nc_def_var( fpImage, + poNetcdfSRS[i].netCDFSRS, + NC_CHAR, + 0, NULL, &NCDFVarID ); - dfNN = adfGeoTransform[3]; - dfSN = ( adfGeoTransform[5] * nYSize ) + dfNN; - dfWE = adfGeoTransform[0]; - dfEE = ( adfGeoTransform[1] * nXSize ) + dfWE; + dfNN = adfGeoTransform[3]; + dfSN = ( adfGeoTransform[5] * nYSize ) + dfNN; + dfWE = adfGeoTransform[0]; + dfEE = ( adfGeoTransform[1] * nXSize ) + dfWE; status = nc_put_att_double( fpImage, - NCDFVarID, - "Northernmost_Northing", - NC_DOUBLE, - 1, - &dfNN ); + NCDFVarID, + "Northernmost_Northing", + NC_DOUBLE, + 1, + &dfNN ); status = nc_put_att_double( fpImage, - NCDFVarID, - "Southernmost_Northing", - NC_DOUBLE, - 1, - &dfSN ); - status = nc_put_att_double( fpImage, - NCDFVarID, - "Easternmost_Easting", - NC_DOUBLE, - 1, - &dfEE ); - status = nc_put_att_double( fpImage, - NCDFVarID, - "Westernmost_Easting", - NC_DOUBLE, - 1, - &dfWE ); - pszWKT = (char *) poSrcDS->GetProjectionRef() ; + NCDFVarID, + "Southernmost_Northing", + NC_DOUBLE, + 1, + &dfSN ); + status = nc_put_att_double( fpImage, + NCDFVarID, + "Easternmost_Easting", + NC_DOUBLE, + 1, + &dfEE ); + status = nc_put_att_double( fpImage, + NCDFVarID, + "Westernmost_Easting", + NC_DOUBLE, + 1, + &dfWE ); + pszWKT = (char *) poSrcDS->GetProjectionRef() ; - nc_put_att_text( fpImage, - NCDFVarID, - "spatial_ref", - strlen( pszWKT ), - pszWKT ); + nc_put_att_text( fpImage, + NCDFVarID, + "spatial_ref", + strlen( pszWKT ), + pszWKT ); - nc_put_att_text( fpImage, - NCDFVarID, - "GeoTransform", - strlen( szGeoTransform ), - szGeoTransform ); + nc_put_att_text( fpImage, + NCDFVarID, + "GeoTransform", + strlen( szGeoTransform ), + szGeoTransform ); - nc_put_att_text( fpImage, - NCDFVarID, - GRD_MAPPING_NAME, - strlen( pszNetcdfProjection ), - pszNetcdfProjection ); + nc_put_att_text( fpImage, + NCDFVarID, + GRD_MAPPING_NAME, + strlen( pszNetcdfProjection ), + pszNetcdfProjection ); for( int iChild = 0; iChild < poPROJCS->GetChildCount(); iChild++ ) { const OGR_SRSNode *poNode; - float fValue; + float fValue; poNode = poPROJCS->GetChild( iChild ); - if( !EQUAL(poNode->GetValue(),"PARAMETER") - || poNode->GetChildCount() != 2 ) + if( !EQUAL(poNode->GetValue(),"PARAMETER") + || poNode->GetChildCount() != 2 ) continue; -/* -------------------------------------------------------------------- */ -/* Look for projection attributes */ -/* -------------------------------------------------------------------- */ - pszParamStr = poNode->GetChild(0)->GetValue(); - pszParamVal = poNode->GetChild(1)->GetValue(); - + /* -------------------------------------------------------------------- */ + /* Look for projection attributes */ + /* -------------------------------------------------------------------- */ + pszParamStr = poNode->GetChild(0)->GetValue(); + pszParamVal = poNode->GetChild(1)->GetValue(); - for(i=0; poNetcdfSRS[i].netCDFSRS != NULL; i++ ) { - if( EQUAL( poNetcdfSRS[i].SRS, pszParamStr ) ) { - CPLDebug( "GDAL_netCDF", "%s = %s", - poNetcdfSRS[i].netCDFSRS, - pszParamVal ); - break; - } - } -/* -------------------------------------------------------------------- */ -/* Write Projection attribute */ -/* -------------------------------------------------------------------- */ - sscanf( pszParamVal, "%f", &fValue ); - if( poNetcdfSRS[i].netCDFSRS != NULL ) { - nc_put_att_float( fpImage, - NCDFVarID, - poNetcdfSRS[i].netCDFSRS, - NC_FLOAT, - 1, - &fValue ); - } + for(i=0; poNetcdfSRS[i].netCDFSRS != NULL; i++ ) { + if( EQUAL( poNetcdfSRS[i].SRS, pszParamStr ) ) { + CPLDebug( "GDAL_netCDF", "%s = %s", + poNetcdfSRS[i].netCDFSRS, + pszParamVal ); + break; + } + } + /* -------------------------------------------------------------------- */ + /* Write Projection attribute */ + /* -------------------------------------------------------------------- */ + sscanf( pszParamVal, "%f", &fValue ); + if( poNetcdfSRS[i].netCDFSRS != NULL ) { + nc_put_att_float( fpImage, + NCDFVarID, + poNetcdfSRS[i].netCDFSRS, + NC_FLOAT, + 1, + &fValue ); + + } } } -/* -------------------------------------------------------------------- */ -/* Initialize Band Map */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Initialize Band Map */ + /* -------------------------------------------------------------------- */ for(int j=1; j <= nBands; j++ ) { - anBandMap[j-1]=j; + anBandMap[j-1]=j; } - -/* -------------------------------------------------------------------- */ -/* Create netCDF variable */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Create netCDF variable */ + /* -------------------------------------------------------------------- */ + for( int i=1; i <= nBands; i++ ) { - char szBandName[ NC_MAX_NAME ]; - GByte *pabScanline = NULL; - GInt16 *pasScanline = NULL; - GInt32 *panScanline = NULL; - float *pafScanline = NULL; - double *padScanline = NULL; - int NCDFVarID; - size_t start[ GDALNBDIM ]; - size_t count[ GDALNBDIM ]; - size_t nBandNameLen; - double dfNoDataValue; - unsigned char cNoDataValue; - float fNoDataValue; - int nlNoDataValue; - short nsNoDataValue; - GDALRasterBandH hBand; + char szBandName[ NC_MAX_NAME ]; + GByte *pabScanline = NULL; + GInt16 *pasScanline = NULL; + GInt32 *panScanline = NULL; + float *pafScanline = NULL; + double *padScanline = NULL; + int NCDFVarID; + size_t start[ GDALNBDIM ]; + size_t count[ GDALNBDIM ]; + size_t nBandNameLen; + double dfNoDataValue; + unsigned char cNoDataValue; + float fNoDataValue; + int nlNoDataValue; + short nsNoDataValue; + GDALRasterBandH hBand; - GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand( i ); + GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand( i ); hBand = GDALGetRasterBand( poSrcDS, i ); - sprintf( szBandName, "Band%d", i ); + sprintf( szBandName, "Band%d", i ); - eDT = poSrcDS->GetRasterBand(i)->GetRasterDataType(); - anBandDims[0] = nYDimID; - anBandDims[1] = nXDimID; - CPLErr eErr = CE_None; + eDT = poSrcDS->GetRasterBand(i)->GetRasterDataType(); + anBandDims[0] = nYDimID; + anBandDims[1] = nXDimID; + CPLErr eErr = CE_None; - dfNoDataValue = poSrcBand->GetNoDataValue(0); + dfNoDataValue = poSrcBand->GetNoDataValue(0); - if( eDT == GDT_Byte ) { - CPLDebug( "GDAL_netCDF", "%s = GDT_Byte ", szBandName ); + if( eDT == GDT_Byte ) { + CPLDebug( "GDAL_netCDF", "%s = GDT_Byte ", szBandName ); - status = nc_def_var( fpImage, szBandName, NC_BYTE, - GDALNBDIM, anBandDims, &NCDFVarID ); + status = nc_def_var( fpImage, szBandName, NC_BYTE, + GDALNBDIM, anBandDims, &NCDFVarID ); -/* -------------------------------------------------------------------- */ -/* Write data line per line */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Write data line per line */ + /* -------------------------------------------------------------------- */ - pabScanline = (GByte *) CPLMalloc( nBands * nXSize * nYSize * - sizeof( GByte ) ); - for( int iLine = 0; iLine < nYSize && eErr == CE_None; iLine++ ) { + pabScanline = (GByte *) CPLMalloc( nBands * nXSize * nYSize * + sizeof( GByte ) ); + for( int iLine = 0; iLine < nYSize && eErr == CE_None; iLine++ ) { -/* -------------------------------------------------------------------- */ -/* Read data from band i */ -/* -------------------------------------------------------------------- */ - eErr = poSrcBand->RasterIO( GF_Read, 0, iLine, nXSize, 1, - pabScanline, nXSize, 1, GDT_Byte, - 0,0); + /* -------------------------------------------------------------------- */ + /* Read data from band i */ + /* -------------------------------------------------------------------- */ + eErr = poSrcBand->RasterIO( GF_Read, 0, iLine, nXSize, 1, + pabScanline, nXSize, 1, GDT_Byte, + 0,0); -/* -------------------------------------------------------------------- */ -/* Write Fill Value */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Write Fill Value */ + /* -------------------------------------------------------------------- */ - cNoDataValue=(unsigned char) dfNoDataValue; - nc_put_att_uchar( fpImage, - NCDFVarID, - _FillValue, - NC_CHAR, - 1, - &cNoDataValue ); - -/* -------------------------------------------------------------------- */ -/* Write Data from Band i */ -/* -------------------------------------------------------------------- */ - start[0]=iLine; - start[1]=0; - count[0]=1; - count[1]=nXSize; + cNoDataValue=(unsigned char) dfNoDataValue; + nc_put_att_uchar( fpImage, + NCDFVarID, + _FillValue, + NC_CHAR, + 1, + &cNoDataValue ); -/* -------------------------------------------------------------------- */ -/* Put NetCDF file in data mode. */ -/* -------------------------------------------------------------------- */ - status = nc_enddef( fpImage ); - status = nc_put_vara_uchar (fpImage, NCDFVarID, start, - count, pabScanline); + /* -------------------------------------------------------------------- */ + /* Write Data from Band i */ + /* -------------------------------------------------------------------- */ + start[0]=iLine; + start[1]=0; + count[0]=1; + count[1]=nXSize; + /* -------------------------------------------------------------------- */ + /* Put NetCDF file in data mode. */ + /* -------------------------------------------------------------------- */ + status = nc_enddef( fpImage ); + status = nc_put_vara_uchar (fpImage, NCDFVarID, start, + count, pabScanline); -/* -------------------------------------------------------------------- */ -/* Put NetCDF file back in define mode. */ -/* -------------------------------------------------------------------- */ - status = nc_redef( fpImage ); - - } - CPLFree( pabScanline ); -/* -------------------------------------------------------------------- */ -/* Int16 */ -/* -------------------------------------------------------------------- */ - } else if( ( eDT == GDT_UInt16 ) || ( eDT == GDT_Int16 ) ) { - CPLDebug( "GDAL_netCDF", "%s = GDT_Int16 ",szBandName ); - status = nc_def_var( fpImage, szBandName, NC_SHORT, - GDALNBDIM, anBandDims, &NCDFVarID ); + /* -------------------------------------------------------------------- */ + /* Put NetCDF file back in define mode. */ + /* -------------------------------------------------------------------- */ + status = nc_redef( fpImage ); - pasScanline = (GInt16 *) CPLMalloc( nBands * nXSize * nYSize * - sizeof( GInt16 ) ); -/* -------------------------------------------------------------------- */ -/* Write Fill Value */ -/* -------------------------------------------------------------------- */ - nsNoDataValue= (GInt16) dfNoDataValue; - nc_put_att_short( fpImage, - NCDFVarID, - _FillValue, - NC_SHORT, - 1, - &nsNoDataValue ); - for( int iLine = 0; iLine < nYSize && eErr == CE_None; iLine++ ) { + } + CPLFree( pabScanline ); + /* -------------------------------------------------------------------- */ + /* Int16 */ + /* -------------------------------------------------------------------- */ - eErr = poSrcBand->RasterIO( GF_Read, 0, iLine, nXSize, 1, - pasScanline, nXSize, 1, GDT_Int16, - 0,0); + } else if( ( eDT == GDT_UInt16 ) || ( eDT == GDT_Int16 ) ) { + CPLDebug( "GDAL_netCDF", "%s = GDT_Int16 ",szBandName ); + status = nc_def_var( fpImage, szBandName, NC_SHORT, + GDALNBDIM, anBandDims, &NCDFVarID ); - start[0]=iLine; - start[1]=0; - count[0]=1; - count[1]=nXSize; + pasScanline = (GInt16 *) CPLMalloc( nBands * nXSize * nYSize * + sizeof( GInt16 ) ); + /* -------------------------------------------------------------------- */ + /* Write Fill Value */ + /* -------------------------------------------------------------------- */ + nsNoDataValue= (GInt16) dfNoDataValue; + nc_put_att_short( fpImage, + NCDFVarID, + _FillValue, + NC_SHORT, + 1, + &nsNoDataValue ); + for( int iLine = 0; iLine < nYSize && eErr == CE_None; iLine++ ) { + eErr = poSrcBand->RasterIO( GF_Read, 0, iLine, nXSize, 1, + pasScanline, nXSize, 1, GDT_Int16, + 0,0); - status = nc_enddef( fpImage ); - status = nc_put_vara_short( fpImage, NCDFVarID, start, - count, pasScanline); - status = nc_redef( fpImage ); - } - CPLFree( pasScanline ); -/* -------------------------------------------------------------------- */ -/* Int32 */ -/* -------------------------------------------------------------------- */ + start[0]=iLine; + start[1]=0; + count[0]=1; + count[1]=nXSize; - } else if( (eDT == GDT_UInt32) || (eDT == GDT_Int32) ) { - CPLDebug( "GDAL_netCDF", "%s = GDT_Int32 ",szBandName ); - status = nc_def_var( fpImage, szBandName, NC_INT, - GDALNBDIM, anBandDims, &NCDFVarID ); - panScanline = (GInt32 *) CPLMalloc( nBands * nXSize * nYSize * - sizeof( GInt32 ) ); -/* -------------------------------------------------------------------- */ -/* Write Fill Value */ -/* -------------------------------------------------------------------- */ - nlNoDataValue= (GInt32) dfNoDataValue; + status = nc_enddef( fpImage ); + status = nc_put_vara_short( fpImage, NCDFVarID, start, + count, pasScanline); + status = nc_redef( fpImage ); + } + CPLFree( pasScanline ); + /* -------------------------------------------------------------------- */ + /* Int32 */ + /* -------------------------------------------------------------------- */ - nc_put_att_int( fpImage, - NCDFVarID, - _FillValue, - NC_INT, - 1, - &nlNoDataValue ); + } else if( (eDT == GDT_UInt32) || (eDT == GDT_Int32) ) { + CPLDebug( "GDAL_netCDF", "%s = GDT_Int32 ",szBandName ); + status = nc_def_var( fpImage, szBandName, NC_INT, + GDALNBDIM, anBandDims, &NCDFVarID ); + panScanline = (GInt32 *) CPLMalloc( nBands * nXSize * nYSize * + sizeof( GInt32 ) ); + /* -------------------------------------------------------------------- */ + /* Write Fill Value */ + /* -------------------------------------------------------------------- */ + nlNoDataValue= (GInt32) dfNoDataValue; - for( int iLine = 0; iLine < nYSize && eErr == CE_None; iLine++ ) { + nc_put_att_int( fpImage, + NCDFVarID, + _FillValue, + NC_INT, + 1, + &nlNoDataValue ); - eErr = poSrcBand->RasterIO( GF_Read, 0, iLine, nXSize, 1, - panScanline, nXSize, 1, GDT_Int32, - 0,0); - start[0]=iLine; - start[1]=0; - count[0]=1; - count[1]=nXSize; + for( int iLine = 0; iLine < nYSize && eErr == CE_None; iLine++ ) { + eErr = poSrcBand->RasterIO( GF_Read, 0, iLine, nXSize, 1, + panScanline, nXSize, 1, GDT_Int32, + 0,0); - status = nc_enddef( fpImage ); - status = nc_put_vara_int( fpImage, NCDFVarID, start, - count, panScanline); - status = nc_redef( fpImage ); - } - CPLFree( panScanline ); -/* -------------------------------------------------------------------- */ -/* float */ -/* -------------------------------------------------------------------- */ - } else if( (eDT == GDT_Float32) ) { - CPLDebug( "GDAL_netCDF", "%s = GDT_Float32 ",szBandName ); - status = nc_def_var( fpImage, szBandName, NC_FLOAT, - GDALNBDIM, anBandDims, &NCDFVarID ); + start[0]=iLine; + start[1]=0; + count[0]=1; + count[1]=nXSize; - pafScanline = (float *) CPLMalloc( nBands * nXSize * nYSize * - sizeof( float ) ); -/* -------------------------------------------------------------------- */ -/* Write Fill Value */ -/* -------------------------------------------------------------------- */ - fNoDataValue= (float) dfNoDataValue; - nc_put_att_float( fpImage, - NCDFVarID, - _FillValue, - NC_FLOAT, - 1, - &fNoDataValue ); - - for( int iLine = 0; iLine < nYSize && eErr == CE_None; iLine++ ) { + status = nc_enddef( fpImage ); + status = nc_put_vara_int( fpImage, NCDFVarID, start, + count, panScanline); + status = nc_redef( fpImage ); + } + CPLFree( panScanline ); + /* -------------------------------------------------------------------- */ + /* float */ + /* -------------------------------------------------------------------- */ + } else if( (eDT == GDT_Float32) ) { + CPLDebug( "GDAL_netCDF", "%s = GDT_Float32 ",szBandName ); + status = nc_def_var( fpImage, szBandName, NC_FLOAT, + GDALNBDIM, anBandDims, &NCDFVarID ); - eErr = poSrcBand->RasterIO( GF_Read, 0, iLine, nXSize, 1, - pafScanline, nXSize, 1, - GDT_Float32, - 0,0); + pafScanline = (float *) CPLMalloc( nBands * nXSize * nYSize * + sizeof( float ) ); - start[0]=iLine; - start[1]=0; - count[0]=1; - count[1]=nXSize; + /* -------------------------------------------------------------------- */ + /* Write Fill Value */ + /* -------------------------------------------------------------------- */ + fNoDataValue= (float) dfNoDataValue; + nc_put_att_float( fpImage, + NCDFVarID, + _FillValue, + NC_FLOAT, + 1, + &fNoDataValue ); + for( int iLine = 0; iLine < nYSize && eErr == CE_None; iLine++ ) { - status = nc_enddef( fpImage ); - status = nc_put_vara_float( fpImage, NCDFVarID, start, - count, pafScanline); - status = nc_redef( fpImage ); - } - CPLFree( pafScanline ); -/* -------------------------------------------------------------------- */ -/* double */ -/* -------------------------------------------------------------------- */ - } else if( (eDT == GDT_Float64) ) { - CPLDebug( "GDAL_netCDF", "%s = GDT_Float64 ",szBandName ); - status = nc_def_var( fpImage, szBandName, NC_DOUBLE, - GDALNBDIM, anBandDims, &NCDFVarID ); + eErr = poSrcBand->RasterIO( GF_Read, 0, iLine, nXSize, 1, + pafScanline, nXSize, 1, + GDT_Float32, + 0,0); - padScanline = (double *) CPLMalloc( nBands * nXSize * nYSize * - sizeof( double ) ); -/* -------------------------------------------------------------------- */ -/* Write Fill Value */ -/* -------------------------------------------------------------------- */ - - nc_put_att_double( fpImage, - NCDFVarID, - _FillValue, - NC_DOUBLE, - 1, - &dfNoDataValue ); + start[0]=iLine; + start[1]=0; + count[0]=1; + count[1]=nXSize; - for( int iLine = 0; iLine < nYSize && eErr == CE_None; iLine++ ) { - eErr = poSrcBand->RasterIO( GF_Read, 0, iLine, nXSize, 1, - padScanline, nXSize, 1, - GDT_Float64, - 0,0); + status = nc_enddef( fpImage ); + status = nc_put_vara_float( fpImage, NCDFVarID, start, + count, pafScanline); + status = nc_redef( fpImage ); + } + CPLFree( pafScanline ); + /* -------------------------------------------------------------------- */ + /* double */ + /* -------------------------------------------------------------------- */ + } else if( (eDT == GDT_Float64) ) { + CPLDebug( "GDAL_netCDF", "%s = GDT_Float64 ",szBandName ); + status = nc_def_var( fpImage, szBandName, NC_DOUBLE, + GDALNBDIM, anBandDims, &NCDFVarID ); - start[0]=iLine; - start[1]=0; - count[0]=1; - count[1]=nXSize; + padScanline = (double *) CPLMalloc( nBands * nXSize * nYSize * + sizeof( double ) ); + /* -------------------------------------------------------------------- */ + /* Write Fill Value */ + /* -------------------------------------------------------------------- */ + nc_put_att_double( fpImage, + NCDFVarID, + _FillValue, + NC_DOUBLE, + 1, + &dfNoDataValue ); - status = nc_enddef( fpImage ); - status = nc_put_vara_double( fpImage, NCDFVarID, start, - count, padScanline); - status = nc_redef( fpImage ); - } - CPLFree( padScanline ); - } - -/* -------------------------------------------------------------------- */ -/* Write Projection for band */ -/* -------------------------------------------------------------------- */ - if( bWriteGeoTransform == TRUE ) { - /* nc_put_att_text( fpImage, NCDFVarID, - COORDINATES, - 7, - LONLAT ); - */ - nc_put_att_text( fpImage, NCDFVarID, - GRD_MAPPING, - strlen( pszNetcdfProjection ), - pszNetcdfProjection ); - } + for( int iLine = 0; iLine < nYSize && eErr == CE_None; iLine++ ) { - sprintf( szBandName, "GDAL Band Number %d", i); - nBandNameLen = strlen( szBandName ); - nc_put_att_text( fpImage, - NCDFVarID, - "long_name", - nBandNameLen, - szBandName ); + eErr = poSrcBand->RasterIO( GF_Read, 0, iLine, nXSize, 1, + padScanline, nXSize, 1, + GDT_Float64, + 0,0); - CopyMetadata( (void *) hBand, fpImage, NCDFVarID ); + start[0]=iLine; + start[1]=0; + count[0]=1; + count[1]=nXSize; - + status = nc_enddef( fpImage ); + status = nc_put_vara_double( fpImage, NCDFVarID, start, + count, padScanline); + status = nc_redef( fpImage ); + } + CPLFree( padScanline ); + } + + /* -------------------------------------------------------------------- */ + /* Write Projection for band */ + /* -------------------------------------------------------------------- */ + if( bWriteGeoTransform == TRUE ) { + /* nc_put_att_text( fpImage, NCDFVarID, + COORDINATES, + 7, + LONLAT ); + */ + nc_put_att_text( fpImage, NCDFVarID, + GRD_MAPPING, + strlen( pszNetcdfProjection ), + pszNetcdfProjection ); + } + + sprintf( szBandName, "GDAL Band Number %d", i); + nBandNameLen = strlen( szBandName ); + nc_put_att_text( fpImage, + NCDFVarID, + "long_name", + nBandNameLen, + szBandName ); + + CopyMetadata( (void *) hBand, fpImage, NCDFVarID ); + + + } @@ -2251,15 +2563,15 @@ // poDstDS->SetGeoTransform( adfGeoTransform ); -/* -------------------------------------------------------------------- */ -/* Cleanup and close. */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Cleanup and close. */ + /* -------------------------------------------------------------------- */ // CPLFree( pabScanline ); nc_close( fpImage ); -/* -------------------------------------------------------------------- */ -/* Re-open dataset, and copy any auxilary pam information. */ -/* -------------------------------------------------------------------- */ + /* -------------------------------------------------------------------- */ + /* Re-open dataset, and copy any auxilary pam information. */ + /* -------------------------------------------------------------------- */ netCDFDataset *poDS = (netCDFDataset *) GDALOpen( pszFilename, GA_ReadOnly ); if( poDS ) @@ -2268,7 +2580,7 @@ return poDS; - } +} /************************************************************************/ /* GDALRegister_netCDF() */ @@ -2279,17 +2591,14 @@ { GDALDriver *poDriver; - if (! GDAL_CHECK_VERSION("netCDF driver")) - return; - if( GDALGetDriverByName( "netCDF" ) == NULL ) { poDriver = new GDALDriver( ); - + poDriver->SetDescription( "netCDF" ); - poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, + poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, "Network Common Data Format" ); - poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, + poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, "frmt_various.html#netCDF" ); poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "nc" ); @@ -2308,7 +2617,7 @@ /* -------------------------------------------------------------------- */ - + //Albers equal area // //grid_mapping_name = albers_conical_equal_area @@ -2388,18 +2697,18 @@ // //Grid mapping attributes // -//false_easting -//false_northing -//grid_mapping_name +//false_easting +//false_northing +//grid_mapping_name //grid_north_pole_latitude //grid_north_pole_longitude -//latitude_of_projection_origin -//longitude_of_central_meridian +//latitude_of_projection_origin +//longitude_of_central_meridian //longitude_of_projection_origin -//north_pole_grid_longitude -//scale_factor_at_central_meridian -//scale_factor_at_projection_origin -//standard_parallel -//straight_vertical_longitude_from_pole +//north_pole_grid_longitude +//scale_factor_at_central_meridian +//scale_factor_at_projection_origin +//standard_parallel +//straight_vertical_longitude_from_pole