--- netcdfdataset.cpp.bak 2011-08-24 15:08:33.000000000 -0300 +++ netcdfdataset.cpp 2011-08-24 16:19:42.000000000 -0300 @@ -483,9 +483,11 @@ } else { switch( vartype ) { case NC_BYTE: - /* don't do default fill-values for bytes, too risky */ - dfNoData = 0.0; - break; + /* 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; break; @@ -1486,6 +1488,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] ) { @@ -1569,6 +1573,15 @@ 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); @@ -2057,14 +2070,25 @@ 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 ) { 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. */ /* -------------------------------------------------------------------- */ @@ -2338,7 +2362,6 @@ /* CreateCopy() */ /************************************************************************/ - static GDALDataset* NCDFCreateCopy( const char * pszFilename, GDALDataset *poSrcDS, int bStrict, char ** papszOptions, @@ -2357,6 +2380,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 +2409,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 +2438,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; @@ -2441,83 +2476,184 @@ char szGeoTransform[ MAX_STR_LEN ]; char szTemp[ 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++ ) { @@ -2678,15 +2815,14 @@ size_t start[ GDALNBDIM ]; size_t count[ GDALNBDIM ]; size_t nBandNameLen; - double dfNoDataValue; + double dfNoDataValue; unsigned char cNoDataValue; float fNoDataValue; int nlNoDataValue; short nsNoDataValue; - GDALRasterBandH hBand; GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand( i ); - hBand = GDALGetRasterBand( poSrcDS, i ); + GDALRasterBandH hBand = GDALGetRasterBand( poSrcDS, i ); sprintf( szBandName, "Band%d", i ); @@ -2733,8 +2869,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 +2882,6 @@ status = nc_enddef( fpImage ); status = nc_put_vara_uchar (fpImage, NCDFVarID, start, count, pabScanline); - - /* -------------------------------------------------------------------- */ /* Put NetCDF file back in define mode. */ /* -------------------------------------------------------------------- */ @@ -2780,8 +2916,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 +2959,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 +3002,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 +3045,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,11 +3070,11 @@ 7, LONLAT ); */ - nc_put_att_text( fpImage, NCDFVarID, - GRD_MAPPING, - strlen( pszNetcdfProjection ), - pszNetcdfProjection ); - } + nc_put_att_text( fpImage, NCDFVarID, + GRD_MAPPING, + strlen( pszNetcdfProjection ), + pszNetcdfProjection ); + } sprintf( szBandName, "GDAL Band Number %d", i); nBandNameLen = strlen( szBandName ); @@ -2942,13 +3085,8 @@ szBandName ); CopyMetadata( (void *) hBand, fpImage, NCDFVarID ); - - - } - - // poDstDS->SetGeoTransform( adfGeoTransform ); @@ -2969,6 +3107,7 @@ return poDS; } + /************************************************************************/ /* GDALRegister_netCDF() */ /************************************************************************/