--- netcdfdataset.cpp.orig 2011-04-15 14:49:17.000000000 -0300 +++ netcdfdataset.cpp 2011-04-19 17:28:51.000000000 -0300 @@ -331,7 +331,60 @@ Taken += result * Sum; } +/* etiennesky fix for netCDF import/export */ +/* -------------------------------------------------------------------- */ +/* Get all other metadata */ +/* -------------------------------------------------------------------- */ + + int nAtt=0; + //char attstr[NC_MAX_NAME]; + 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,"NETCDF_%s",szTemp); + 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,"%g",fval); + SetMetadataItem( szMetaTemp, szTemp ); + break; + case NC_DOUBLE: + status = nc_get_att_double( poDS->cdfid, nZId, + szTemp, &dval ); + sprintf( szTemp,"%g",dval); + SetMetadataItem( szMetaTemp, szTemp ); + break; + default: + break; + } + } return CE_None; } @@ -485,6 +538,9 @@ case NC_BYTE: /* don't do default fill-values for bytes, too risky */ dfNoData = 0.0; + /* from etiennesky: 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; @@ -1373,6 +1429,12 @@ //oSRS.SetWellKnownGeogCS( "WGS84" ); } } + /* etiennesky fix for netCDF import/export */ + /* no grid_mapping value - not a GDAL netcdf */ + else if( EQUAL( szDimNameX,"lon" ) ) { + oSRS.SetWellKnownGeogCS( "WGS84" ); + } + /* -------------------------------------------------------------------- */ /* Read projection coordinates */ /* -------------------------------------------------------------------- */ @@ -2402,10 +2464,10 @@ GDALDataType eDT; - status = nc_def_dim( fpImage, "x", nXSize, &nXDimID ); + status = nc_def_dim( fpImage, "lon", nXSize, &nXDimID ); CPLDebug( "GDAL_netCDF", "status nc_def_dim X = %d\n", status ); - status = nc_def_dim( fpImage, "y", nYSize, &nYDimID ); + status = nc_def_dim( fpImage, "lat", nYSize, &nYDimID ); CPLDebug( "GDAL_netCDF", "status nc_def_dim Y = %d\n", status ); CPLDebug( "GDAL_netCDF", "nYDimID = %d\n", nXDimID ); @@ -2440,13 +2502,22 @@ double adfGeoTransform[6]; char szGeoTransform[ MAX_STR_LEN ]; char szTemp[ MAX_STR_LEN ]; + char szGdalGeo[ MAX_STR_LEN ]; - -/* -------------------------------------------------------------------- */ -/* Copy GeoTransform array from source */ -/* -------------------------------------------------------------------- */ + double dfX0=0.0; + double dfDX=0.0; + double dfY0=0.0; + double dfDY=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; *szGeoTransform = '\0'; for( i=0; i<6; i++ ) { @@ -2455,69 +2526,212 @@ 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 */ + /* -------------------------------------------------------------------- */ - nc_put_att_text( fpImage, - NCDFVarID, - "GeoTransform", - strlen( szGeoTransform ), - szGeoTransform ); + if( bWriteGeoTransform == TRUE ) { + + strcpy( pszNetcdfProjection, "GDAL_Geographics" ); + nc_def_var( fpImage, + pszNetcdfProjection, + NC_CHAR, + 0, NULL, &NCDFVarID ); + + nc_put_att_double( fpImage, + NCDFVarID, + "Northernmost_Northing", + NC_DOUBLE, + 1, + &dfNN ); + nc_put_att_double( fpImage, + NCDFVarID, + "Southernmost_Northing", + NC_DOUBLE, + 1, + &dfSN ); + nc_put_att_double( fpImage, + NCDFVarID, + "Easternmost_Easting", + NC_DOUBLE, + 1, + &dfEE ); + nc_put_att_double( fpImage, + NCDFVarID, + "Westernmost_Easting", + NC_DOUBLE, + 1, + &dfWE ); + + 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, + 30, + "Geographics Coordinate System" ); + nc_put_att_text( fpImage, + NCDFVarID, + LNG_NAME, + 14, + "Grid_latitude" ); + } + + /* write geo info to global attribute GDAL_Geographics */ + *szGdalGeo = '\0'; + + sprintf( szTemp, "Northernmost_Northing = %.15g \n", + dfNN ); + strcat( szGdalGeo, szTemp ); + sprintf( szTemp, "Southernmost_Northing = %.15g \n", + dfSN ); + strcat( szGdalGeo, szTemp ); + sprintf( szTemp, "Easternmost_Easting = %.15g \n", + dfEE ); + strcat( szGdalGeo, szTemp ); + sprintf( szTemp, "Westernmost_Easting = %.15g \n", + dfWE ); + strcat( szGdalGeo, szTemp ); + + sprintf( szTemp, "spatial_ref = %s \n", + pszWKT ); + strcat( szGdalGeo, szTemp ); + sprintf( szTemp, "GeoTransform = %s \n", + szGeoTransform ); + strcat( szGdalGeo, szTemp ); + sprintf( szTemp, "%s = Geographics Coordinate System \n", + GRD_MAPPING_NAME ); + strcat( szGdalGeo, szTemp ); + sprintf( szTemp, "%s = Grid_latitude", + LNG_NAME ); + strcat( szGdalGeo, szTemp ); 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_GLOBAL, + "GDAL_Geographics", + strlen( szGdalGeo ), + szGdalGeo ); + + /* -------------------------------------------------------------------- */ + /* 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" ); + nc_put_att_text( fpImage, + NCDFVarID, + "axis", + 1, + "Y" ); + + /* -------------------------------------------------------------------- */ + /* 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; iGetRasterBand( i ); - hBand = GDALGetRasterBand( poSrcDS, i ); + GDALRasterBandH hBand = GDALGetRasterBand( poSrcDS, i ); + + +/* -------------------------------------------------------------------- */ +/* Set var name */ +/* -------------------------------------------------------------------- */ - sprintf( szBandName, "Band%d", i ); + /* etiennesky fix for netCDF import/export */ + + /* 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; @@ -2734,7 +2975,12 @@ /* Write Data from Band i */ /* -------------------------------------------------------------------- */ start[0]=iLine; - start[1]=0; + /* invert latitude values */ + /* must be done for projected also */ + /* NOTE: only tested for BYTE and Int16 data types */ + if( oSRS.IsGeographic() ) + start[0]=nYSize - iLine - 1; + start[1]=0; count[0]=1; count[1]=nXSize; @@ -2745,6 +2991,10 @@ status = nc_put_vara_uchar (fpImage, NCDFVarID, start, count, pabScanline); + if (status != NC_NOERR) + CPLError( CE_Failure, CPLE_AppDefined, + "netCDF write failed: %s", + nc_strerror( status ) ); /* -------------------------------------------------------------------- */ /* Put NetCDF file back in define mode. */ @@ -2781,6 +3031,11 @@ 0,0); start[0]=iLine; + /* invert latitude values */ + /* must be done for projected also */ + /* NOTE: only tested for BYTE and Int16 data types */ + if( oSRS.IsGeographic() ) + start[0]=nYSize - iLine - 1; start[1]=0; count[0]=1; count[1]=nXSize; @@ -2823,6 +3078,11 @@ 0,0); start[0]=iLine; + /* invert latitude values */ + /* must be done for projected also */ + /* NOTE: only tested for BYTE and Int16 data types */ + if( oSRS.IsGeographic() ) + start[0]=nYSize - iLine - 1; start[1]=0; count[0]=1; count[1]=nXSize; @@ -2864,6 +3124,11 @@ 0,0); start[0]=iLine; + /* invert latitude values */ + /* must be done for projected also */ + /* NOTE: only tested for BYTE and Int16 data types */ + if( oSRS.IsGeographic() ) + start[0]=nYSize - iLine - 1; start[1]=0; count[0]=1; count[1]=nXSize; @@ -2905,6 +3170,11 @@ 0,0); start[0]=iLine; + /* invert latitude values */ + /* must be done for projected also */ + /* NOTE: only tested for BYTE and Int16 data types */ + if( oSRS.IsGeographic() ) + start[0]=nYSize - iLine - 1; start[1]=0; count[0]=1; count[1]=nXSize; @@ -2932,23 +3202,21 @@ 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 );