--- netcdfdataset.cpp.bak 2011-08-24 15:08:33.000000000 -0300 +++ netcdfdataset.cpp 2011-08-24 16:02:58.000000000 -0300 @@ -31,7 +31,6 @@ #include "cpl_error.h" CPL_CVSID("$Id: netcdfdataset.cpp 20997 2010-10-28 19:10:19Z kyle $"); - /************************************************************************/ /* ==================================================================== */ /* netCDFRasterBand */ @@ -185,7 +184,7 @@ { char szVarName[NC_MAX_NAME]; char szMetaName[NC_MAX_NAME]; - char szMetaTemp[8192]; + char szMetaTemp[MAX_STR_LEN]; int nd; int i,j; int Sum = 1; @@ -293,7 +292,7 @@ status = nc_get_vara_double( poDS->cdfid, nVarID, start, count, &dfData); - sprintf( szMetaTemp,"%g", dfData ); + sprintf( szMetaTemp,"%.16g", dfData ); break; default: break; @@ -331,7 +330,58 @@ Taken += result * Sum; } +/* etiennesky fix for netCDF import/export */ +/* -------------------------------------------------------------------- */ +/* Get all other metadata */ +/* -------------------------------------------------------------------- */ + + int nAtt=0; + nc_type atttype=NC_NAT; + size_t attlen; + float fval; + double dval; + int ival; + + nc_inq_varnatts( poDS->cdfid, nZId, &nAtt ); + for( i=0; i < nAtt ; i++ ) { + status = nc_inq_attname( poDS->cdfid, nZId, + i, szTemp); + status = nc_inq_att( poDS->cdfid, nZId, + szTemp, &atttype, &attlen); + if(strcmp(szTemp,_FillValue) ==0) continue; + sprintf( szMetaTemp,"%s",szTemp); + switch( atttype ) { + case NC_CHAR: + char *fillc; + fillc = (char *) CPLCalloc( attlen+1, sizeof(char) ); + status=nc_get_att_text( poDS->cdfid, nZId, + szTemp, fillc ); + SetMetadataItem( szMetaTemp, fillc ); + CPLFree(fillc); + break; + case NC_INT: + status = nc_get_att_int( poDS->cdfid, nZId, + szTemp, &ival ); + sprintf( szTemp,"%d",ival); + SetMetadataItem( szMetaTemp, szTemp ); + break; + case NC_FLOAT: + status = nc_get_att_float( poDS->cdfid, nZId, + szTemp, &fval ); + sprintf( szTemp,"%f",fval); + SetMetadataItem( szMetaTemp, szTemp ); + break; + case NC_DOUBLE: + status = nc_get_att_double( poDS->cdfid, nZId, + szTemp, &dval ); + sprintf( szTemp,"%.16g",dval); + SetMetadataItem( szMetaTemp, szTemp ); + break; + default: + break; + } + } return CE_None; } @@ -355,7 +405,7 @@ nc_type atttype=NC_NAT; size_t attlen; int status; - char szNoValueName[8192]; + char szNoValueName[MAX_STR_LEN]; this->panBandZPos = NULL; @@ -483,8 +533,10 @@ } else { switch( vartype ) { case NC_BYTE: - /* don't do default fill-values for bytes, too risky */ - dfNoData = 0.0; + /* dfNoData = 0.0; */ + /* why is this risky? use default value */ + /* if not values that are 0.0 will be treated as missing */ + dfNoData = NC_FILL_BYTE; break; case NC_CHAR: dfNoData = NC_FILL_CHAR; @@ -826,7 +878,8 @@ nc_inq_varname( cdfid, var, szVarName ); strcpy(szTemp,szVarName); - strcat(szTemp,"#grid_mapping"); + strcat(szTemp,"#"); + strcat(szTemp,GRD_MAPPING); pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); if( pszValue ) { strcpy(szGridMappingName,szTemp); @@ -852,7 +905,7 @@ szDimNameY[3] = '\0'; /* -------------------------------------------------------------------- */ -/* Read grid_mappinginformation and set projections */ +/* Read grid_mapping information and set projections */ /* -------------------------------------------------------------------- */ if( !( EQUAL(szGridMappingName,"" ) ) ) { @@ -1363,7 +1416,7 @@ /* -------------------------------------------------------------------- */ /* Is this Latitude/Longitude Grid, default */ /* -------------------------------------------------------------------- */ - + } else if( EQUAL( szDimNameX,"lon" ) ) { oSRS.SetWellKnownGeogCS( "WGS84" ); @@ -1373,13 +1426,13 @@ //oSRS.SetWellKnownGeogCS( "WGS84" ); } } + /* -------------------------------------------------------------------- */ /* 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) ); @@ -1431,7 +1484,6 @@ /* For Latitude we allow an error of 0.1 degrees for gaussion */ /* gridding */ /* -------------------------------------------------------------------- */ - if((( abs( abs(nSpacingBegin) - abs(nSpacingLast) ) ) < 100 ) && (( abs( abs(nSpacingBegin) - abs(nSpacingMiddle) ) ) < 100 ) && (( abs( abs(nSpacingMiddle) - abs(nSpacingLast) ) ) < 100) ) { @@ -1486,6 +1538,8 @@ if( ( EQUAL( szDimNameY, "lat" ) || EQUAL( szDimNameY, "y" ) ) && pdfYCoord[0] < pdfYCoord[1] ) poDS->bBottomUp = TRUE; + else + poDS->bBottomUp = FALSE; // Check for reverse order of y-coordinate if ( yMinMax[0] > yMinMax[1] ) { @@ -1568,7 +1622,16 @@ pszWKT = CSLFetchNameValue(poDS->papszMetadata, szTemp); if( pszWKT != NULL ) { - + + /* Check if bottom up was previously set by GDAL */ + poDS->bBottomUp = FALSE; + strcpy(szTemp,szGridMappingValue); + strcat( szTemp, "#" ); + strcat( szTemp, "GDAL_Bottom_Up"); + pszValue = CSLFetchNameValue(poDS->papszMetadata, szTemp); + if( pszValue != NULL && strncmp(pszValue,"1",1) == 0 ) + poDS->bBottomUp = TRUE; + pszProjection = CPLStrdup( pszWKT ); strcpy(szTemp,szGridMappingValue); @@ -1743,7 +1806,7 @@ nMetaTempSize = nAttrLen + 1; pszMetaTemp = (char *) CPLCalloc( nMetaTempSize, sizeof( char )); *pszMetaTemp = '\0'; - + switch (nAttrType) { case NC_CHAR: nc_get_att_text( cdfid, var, szAttrName, pszMetaTemp ); @@ -1778,10 +1841,10 @@ 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] ); + sprintf( szTemp, "%f, ", pfTemp[m] ); SafeStrcat(&pszMetaTemp, szTemp, &nMetaTempSize); } - sprintf( szTemp, "%e", pfTemp[m] ); + sprintf( szTemp, "%f", pfTemp[m] ); SafeStrcat(&pszMetaTemp,szTemp, &nMetaTempSize); CPLFree(pfTemp); break; @@ -1790,10 +1853,10 @@ 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] ); + sprintf( szTemp, "%.16g, ", pdfTemp[m] ); SafeStrcat(&pszMetaTemp, szTemp, &nMetaTempSize); } - sprintf( szTemp, "%g", pdfTemp[m] ); + sprintf( szTemp, "%.16g", pdfTemp[m] ); SafeStrcat(&pszMetaTemp, szTemp, &nMetaTempSize); CPLFree(pdfTemp); break; @@ -2057,13 +2120,24 @@ CPLDebug( "GDAL_netCDF", "dim_count = %d\n", dim_count ); +/* -------------------------------------------------------------------- */ +/* Check Conventions attribute */ +/* note that 'Conventions' is always capital 'C' in CF spec */ +/* For CF-1/COARDS assume bottom first */ +/* -------------------------------------------------------------------- */ + poDS->bBottomUp = FALSE; if( (status = nc_get_att_text( cdfid, NC_GLOBAL, "Conventions", - attname )) != NC_NOERR ) { + attname )) != NC_NOERR ) + { CPLError( CE_Warning, CPLE_AppDefined, "No UNIDATA NC_GLOBAL:Conventions attribute"); - /* note that 'Conventions' is always capital 'C' in CF spec*/ } - + else + { + if ( strncmp(attname,"CF-1",4) == 0 || + strncmp(attname,"COARDS",6) == 0 ) + poDS->bBottomUp = TRUE; + } /* -------------------------------------------------------------------- */ /* Create band information objects. */ @@ -2281,6 +2355,7 @@ const char *pszField; char szMetaName[ MAX_STR_LEN ]; char szMetaValue[ MAX_STR_LEN ]; + char szTemp[ MAX_STR_LEN ]; int nDataLength; int nItems; @@ -2312,14 +2387,33 @@ if( papszFieldData[1] != NULL ) { strcpy( szMetaName, papszFieldData[ 0 ] ); strcpy( szMetaValue, papszFieldData[ 1 ] ); + + if( CDFVarID == NC_GLOBAL ) { -/* -------------------------------------------------------------------- */ -/* netCDF attributes do not like the '#' character. */ -/* -------------------------------------------------------------------- */ - - for( unsigned int h=0; h < strlen( szMetaName ) -1 ; h++ ) { - if( szMetaName[h] == '#' ) szMetaName[h] = '-'; - } + /* Fix NC_GLOBAL, time and GDAL- */ + if( strncmp( szMetaName, "NC_GLOBAL#", 10 ) == 0 ) { + strcpy( szTemp, szMetaName+10 ); + strcpy( szMetaName, szTemp ); + } + else if( strncmp( szMetaName, "time#", 5 ) == 0 ) { + szMetaName[4] = '-'; + } + else if ( strstr( szMetaName, "#" ) == NULL ) { + strcpy( szTemp, "GDAL_" ); + strcat( szTemp, szMetaName ); + strcpy( szMetaName, szTemp ); + } + + } + +/* -------------------------------------------------------------------- */ +/* Only copy data without # */ +/* -------------------------------------------------------------------- */ + // for( unsigned int h=0; h < strlen( szMetaName ) -1 ; h++ ) { + // if( szMetaName[h] == '#' ) szMetaName[h] = '-'; + // } + + if ( strstr( szMetaName, "#" ) == NULL ) { nDataLength = strlen( szMetaValue ); nc_put_att_text( fpImage, @@ -2328,6 +2422,8 @@ nDataLength, szMetaValue ); + } + } CSLDestroy( papszFieldData ); @@ -2357,6 +2453,9 @@ int bWriteGeoTransform = FALSE; char pszNetcdfProjection[ NC_MAX_NAME ]; + char pszXDimName[ MAX_STR_LEN ]; + char pszYDimName[ MAX_STR_LEN ]; + if (nBands == 0) { CPLError( CE_Failure, CPLE_NotSupported, @@ -2383,12 +2482,24 @@ bProgressive = CSLFetchBoolean( papszOptions, "PROGRESSIVE", FALSE ); /* -------------------------------------------------------------------- */ +/* Get Projection ref for netCDF data CF-1 Convention */ +/* -------------------------------------------------------------------- */ + + OGRSpatialReference oSRS; + char *pszWKT = (char *) poSrcDS->GetProjectionRef(); + + if( pszWKT != NULL ) + oSRS.importFromWkt( &pszWKT ); + +/* -------------------------------------------------------------------- */ /* Create the dataset. */ /* -------------------------------------------------------------------- */ + int fpImage; int status; int nXDimID = 0; int nYDimID = 0; + GDALDataType eDT; status = nc_create( pszFilename, NC_CLOBBER, &fpImage ); @@ -2400,34 +2511,31 @@ return NULL; } - GDALDataType eDT; + if( oSRS.IsProjected() ) + { + strcpy( pszXDimName, "x" ); + strcpy( pszYDimName, "y" ); + } + else + { + strcpy( pszXDimName, "lon" ); + strcpy( pszYDimName, "lat" ); + } - status = nc_def_dim( fpImage, "x", nXSize, &nXDimID ); - CPLDebug( "GDAL_netCDF", "status nc_def_dim X = %d\n", status ); + status = nc_def_dim( fpImage, pszXDimName, nXSize, &nXDimID ); + CPLDebug( "GDAL_netCDF", "status nc_def_dim %s = %d\n", pszXDimName, status ); - status = nc_def_dim( fpImage, "y", nYSize, &nYDimID ); - CPLDebug( "GDAL_netCDF", "status nc_def_dim Y = %d\n", status ); + status = nc_def_dim( fpImage, pszYDimName, nYSize, &nYDimID ); + CPLDebug( "GDAL_netCDF", "status nc_def_dim %s = %d\n", pszYDimName, status ); CPLDebug( "GDAL_netCDF", "nYDimID = %d\n", nXDimID ); CPLDebug( "GDAL_netCDF", "nXDimID = %d\n", nYDimID ); CPLDebug( "GDAL_netCDF", "nXSize = %d\n", nXSize ); CPLDebug( "GDAL_netCDF", "nYSize = %d\n", nYSize ); - CopyMetadata((void *) poSrcDS, fpImage, NC_GLOBAL ); -/* -------------------------------------------------------------------- */ -/* Set Projection for netCDF data CF-1 Convention */ -/* -------------------------------------------------------------------- */ - - OGRSpatialReference oSRS; - char *pszWKT = (char *) poSrcDS->GetProjectionRef(); - - - if( pszWKT != NULL ) - oSRS.importFromWkt( &pszWKT ); - - if( oSRS.IsGeographic() ) { + if( ! oSRS.IsProjected() ) { /* assume geographic... */ int status; int i; @@ -2440,84 +2548,186 @@ double adfGeoTransform[6]; char szGeoTransform[ MAX_STR_LEN ]; char szTemp[ MAX_STR_LEN ]; + char szGdalGeo[ MAX_STR_LEN ]; + double dfX0=0.0; + double dfDX=0.0; + double dfY0=0.0; + double dfDY=0.0; + double dfTemp=0.0; + double *pafLonLat = NULL; + int anLatLonDims[1]; + size_t startLonLat[1]; + size_t countLonLat[1]; /* -------------------------------------------------------------------- */ /* Copy GeoTransform array from source */ /* -------------------------------------------------------------------- */ poSrcDS->GetGeoTransform( adfGeoTransform ); - bWriteGeoTransform = TRUE; + bWriteGeoTransform = FALSE; + if ( oSRS.IsGeographic() ) + bWriteGeoTransform = TRUE; *szGeoTransform = '\0'; for( i=0; i<6; i++ ) { - sprintf( szTemp, "%g ", + sprintf( szTemp, "%.16g ", 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_put_att_double( fpImage, - 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() ; - - nc_put_att_text( fpImage, - NCDFVarID, - "spatial_ref", - strlen( pszWKT ), - pszWKT ); + /* -------------------------------------------------------------------- */ + /* Write GDAL_Geographics attributes */ + /* -------------------------------------------------------------------- */ + + if( bWriteGeoTransform == TRUE ) + { - nc_put_att_text( fpImage, - NCDFVarID, - "GeoTransform", - strlen( szGeoTransform ), - szGeoTransform ); + strcpy( pszNetcdfProjection, "crs" ); + nc_def_var( fpImage, + pszNetcdfProjection, + NC_CHAR, + 0, NULL, &NCDFVarID ); + nc_put_att_text( fpImage, + NCDFVarID, + GRD_MAPPING_NAME, + strlen("latitude_longitude"), + "latitude_longitude" ); + dfTemp = oSRS.GetPrimeMeridian(); + nc_put_att_double( fpImage, + NCDFVarID, + "longitude_of_prime_meridian", + NC_DOUBLE, + 1, + &dfTemp ); + dfTemp = oSRS.GetSemiMajor(); + nc_put_att_double( fpImage, + NCDFVarID, + "semi_major_axis", + NC_DOUBLE, + 1, + &dfTemp ); + dfTemp = oSRS.GetInvFlattening(); + nc_put_att_double( fpImage, + NCDFVarID, + "inverse_flattening", + NC_DOUBLE, + 1, + &dfTemp ); + // pszWKT = (char *) poSrcDS->GetProjectionRef(); + // nc_put_att_text( fpImage, + // NCDFVarID, + // "spatial_ref", + // strlen( pszWKT ), + // pszWKT ); + // nc_put_att_text( fpImage, + // NCDFVarID, + // "GeoTransform", + // strlen( szGeoTransform ), + // szGeoTransform ); + } + + /* -------------------------------------------------------------------- */ + /* Write latitude attributes */ + /* -------------------------------------------------------------------- */ + anLatLonDims[0] = nYDimID; + status = nc_def_var( fpImage, "lat", NC_DOUBLE, 1, anLatLonDims, &NCDFVarID ); + nc_put_att_text( fpImage, + NCDFVarID, + "standard_name", + 8, + "latitude" ); + nc_put_att_text( fpImage, + NCDFVarID, + "long_name", + 8, + "latitude" ); + nc_put_att_text( fpImage, + NCDFVarID, + "units", + 13, + "degrees_north" ); + + /* -------------------------------------------------------------------- */ + /* Write latitude values */ + /* -------------------------------------------------------------------- */ + // dfY0 = adfGeoTransform[3]; + /* invert latitude values */ + dfY0 = ( adfGeoTransform[5] * nYSize ) + adfGeoTransform[3]; + dfDY = adfGeoTransform[5]; + + pafLonLat = (double *) CPLMalloc( nYSize * sizeof( double ) ); + for( i=0; iGetProjectionRef() ; - + pszWKT = (char *) poSrcDS->GetProjectionRef(); 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, GRD_MAPPING_NAME, strlen( pszNetcdfProjection ), pszNetcdfProjection ); - + nc_put_att_text( fpImage, + NCDFVarID, + "GDAL-Bottom_Up", + 1, + "1"); for( int iChild = 0; iChild < poPROJCS->GetChildCount(); iChild++ ) { @@ -2677,18 +2888,40 @@ 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; + const char *tmpMetadata; + char szLongName[ NC_MAX_NAME ]; GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand( i ); - hBand = GDALGetRasterBand( poSrcDS, i ); + GDALRasterBandH hBand = GDALGetRasterBand( poSrcDS, i ); + + /* etiennesky fix for netCDF import/export */ - sprintf( szBandName, "Band%d", i ); + /* Get var name from NETCDF_VARNAME */ + tmpMetadata = poSrcBand->GetMetadataItem("NETCDF_VARNAME"); + if( tmpMetadata != NULL) { + if(nBands>1) sprintf(szBandName,"%s%d",tmpMetadata,i); + else strcpy( szBandName, tmpMetadata ); + // poSrcBand->SetMetadataItem("NETCDF_VARNAME",""); + } + else + sprintf( szBandName, "Band%d", i ); + + /* Get long_name from #long_name */ + sprintf(szLongName,"%s#long_name",poSrcBand->GetMetadataItem("NETCDF_VARNAME")); + tmpMetadata = poSrcDS->GetMetadataItem(szLongName); + if( tmpMetadata != NULL) + strcpy( szLongName, tmpMetadata); + else + sprintf( szLongName, "GDAL Band Number %d", i); + +/* -------------------------------------------------------------------- */ +/* Get Data type */ +/* -------------------------------------------------------------------- */ eDT = poSrcDS->GetRasterBand(i)->GetRasterDataType(); anBandDims[0] = nYDimID; @@ -2733,8 +2966,10 @@ /* -------------------------------------------------------------------- */ /* Write Data from Band i */ /* -------------------------------------------------------------------- */ - start[0]=iLine; - start[1]=0; + // start[0]=iLine; + /* invert latitude values */ + start[0]=nYSize - iLine - 1; + start[1]=0; count[0]=1; count[1]=nXSize; @@ -2744,8 +2979,6 @@ status = nc_enddef( fpImage ); status = nc_put_vara_uchar (fpImage, NCDFVarID, start, count, pabScanline); - - /* -------------------------------------------------------------------- */ /* Put NetCDF file back in define mode. */ /* -------------------------------------------------------------------- */ @@ -2780,8 +3013,9 @@ pasScanline, nXSize, 1, GDT_Int16, 0,0); - start[0]=iLine; - start[1]=0; + // start[0]=iLine; + start[0]=nYSize - iLine - 1; + start[1]=0; count[0]=1; count[1]=nXSize; @@ -2822,7 +3056,9 @@ panScanline, nXSize, 1, GDT_Int32, 0,0); - start[0]=iLine; + // start[0]=iLine; + /* invert latitude values */ + start[0]=nYSize - iLine - 1; start[1]=0; count[0]=1; count[1]=nXSize; @@ -2863,7 +3099,9 @@ GDT_Float32, 0,0); - start[0]=iLine; + // start[0]=iLine; + /* invert latitude values */ + start[0]=nYSize - iLine - 1; start[1]=0; count[0]=1; count[1]=nXSize; @@ -2904,7 +3142,9 @@ GDT_Float64, 0,0); - start[0]=iLine; + // start[0]=iLine; + /* invert latitude values */ + start[0]=nYSize - iLine - 1; start[1]=0; count[0]=1; count[1]=nXSize; @@ -2927,28 +3167,24 @@ 7, LONLAT ); */ - nc_put_att_text( fpImage, NCDFVarID, - GRD_MAPPING, - strlen( pszNetcdfProjection ), - pszNetcdfProjection ); - } + nc_put_att_text( fpImage, NCDFVarID, + GRD_MAPPING, + strlen( pszNetcdfProjection ), + pszNetcdfProjection ); + } + +/* -------------------------------------------------------------------- */ +/* Copy Metadata for band */ +/* -------------------------------------------------------------------- */ - sprintf( szBandName, "GDAL Band Number %d", i); - nBandNameLen = strlen( szBandName ); nc_put_att_text( fpImage, NCDFVarID, "long_name", - nBandNameLen, - szBandName ); - + strlen( szLongName ), + szLongName ); CopyMetadata( (void *) hBand, fpImage, NCDFVarID ); - - - } - - // poDstDS->SetGeoTransform( adfGeoTransform );