Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Netcdf 6195 assume longlat #6910

Merged
merged 16 commits into from
Dec 15, 2022
Merged
Show file tree
Hide file tree
Changes from 9 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
32 changes: 32 additions & 0 deletions autotest/gdrivers/netcdf.py
Original file line number Diff line number Diff line change
Expand Up @@ -1107,6 +1107,38 @@ def test_netcdf_27():
gdal.SetConfigOption("GDAL_NETCDF_BOTTOMUP", config_bak)


###############################################################################
# check support for GDAL_NETCDF_ASSUME_LONGLAT configuration option


def test_netcdf_assume_longlat():

# test default config
with gdaltest.config_option("GDAL_NETCDF_ASSUME_LONGLAT", "YES"):
ds = gdal.Open("data/netcdf/trmm-nc2.nc")
srs = ds.GetSpatialRef()
assert srs is not None
assert (
srs.ExportToWkt()
== 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Longitude",EAST],AXIS["Latitude",NORTH]]'
)
mdsumner marked this conversation as resolved.
Show resolved Hide resolved
with gdaltest.config_option("GDAL_NETCDF_ASSUME_LONGLAT", "NO"):
ds = gdal.Open("data/netcdf/trmm-nc2.nc")
assert ds.GetSpatialRef() is None

# test open option and config overrides
with gdaltest.config_option("GDAL_NETCDF_ASSUME_LONGLAT", "YES"):
ds = gdal.OpenEx("data/netcdf/trmm-nc2.nc", open_options=["ASSUME_LONGLAT=NO"])
srs = ds.GetSpatialRef()
assert srs is None
with gdaltest.config_option("GDAL_NETCDF_ASSUME_LONGLAT", "NO"):
ds = gdal.OpenEx("data/netcdf/trmm-nc2.nc", open_options=["ASSUME_LONGLAT=YES"])
assert ds.GetSpatialRef() is not None
assert (
srs.ExportToWkt()
== 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Longitude",EAST],AXIS["Latitude",NORTH]]'
)
mdsumner marked this conversation as resolved.
Show resolved Hide resolved

rouault marked this conversation as resolved.
Show resolved Hide resolved
###############################################################################
# check support for writing multi-dimensional files (helper function)

Expand Down
11 changes: 11 additions & 0 deletions doc/source/drivers/raster/netcdf.rst
Original file line number Diff line number Diff line change
Expand Up @@ -283,6 +283,12 @@ The following open options are available:
same dimensions, then they should be reported as multiple bands of a same dataset.
Default is NO (that is each variable will be reported as a separate
subdataset)

- **ASSUME_LONGLAT=[YES/NO]** : Whether a Geographic CRS should
mdsumner marked this conversation as resolved.
Show resolved Hide resolved
be assumed and applied when, none has otherwise been found, a meaningful
geotransform has been found, and that geotransform is within the bounds
-180,360 -90,90, if YES assume OGC:CRS84. Default is NO.


Creation Issues
---------------
Expand Down Expand Up @@ -510,6 +516,11 @@ Configuration Options
should be always considered as geospatial axis, even if the lack
conventional attributes confirming it. Default is NO.

- **GDAL_NETCDF_ASSUME_LONGLAT=[YES/NO]** : Whether a Geographic CRS should
mdsumner marked this conversation as resolved.
Show resolved Hide resolved
be assumed and applied when, none has otherwise been found, a meaningful
geotransform has been found, and that geotransform is within the bounds
-180,360 -90,90, if YES assume OGC:CRS84. Default is NO.

VSI Virtual File System API support
-----------------------------------

Expand Down
49 changes: 48 additions & 1 deletion frmts/netcdf/netcdfdataset.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -2873,7 +2873,7 @@ netCDFDataset::netCDFDataset() :
nYDimID(-1),
bIsProjected(false),
bIsGeographic(false), // Can be not projected, and also not geographic

mdsumner marked this conversation as resolved.
Show resolved Hide resolved
// State vars.
bDefineMode(true),
bAddedGridMappingRef(false),
Expand Down Expand Up @@ -4777,6 +4777,49 @@ void netCDFDataset::SetProjectionFromVar( int nGroupId, int nVarId,
if( !bGotGeogCS && !bGotCfSRS && !bGotGdalSRS && !bGotCfGT && !bGotCfWktSRS)
CPLDebug("GDAL_netCDF", "did not get projection from CF nor GDAL!");

// wish of 6195
// we don't have CS/SRS, but we do have GT, and we live in -180,360 -90,90
if (!bGotGeogCS && !bGotCfSRS && !bGotGdalSRS && !bGotCfWktSRS)
{
if (bGotCfGT || bGotGdalGT)
{
bool bAssumedLongLat = CPLTestBool(
CSLFetchNameValueDef(papszOpenOptions, "ASSUME_LONGLAT",
CPLGetConfigOption("GDAL_NETCDF_ASSUME_LONGLAT", "NO")));

if (bAssumedLongLat &&
adfTempGeoTransform[0] >= -180 && adfTempGeoTransform[0] < 360 &&
(adfTempGeoTransform[0] + adfTempGeoTransform[1] * poDS->GetRasterXSize()) <= 360 &&
adfTempGeoTransform[3] <= 90 && adfTempGeoTransform[3] > -90 &&
(adfTempGeoTransform[3] + adfTempGeoTransform[5] * poDS->GetRasterYSize()) >= -90)
{

poDS->bIsGeographic = true;
char *pszTempProjection = nullptr;
// seems odd to use 4326 so OGC:CRS84
oSRS.SetFromUserInput("OGC:CRS84");
oSRS.exportToWkt(&pszTempProjection);
if(returnProjStr != nullptr)
{
(*returnProjStr) = std::string(pszTempProjection);
}
else
{
m_bAddedProjectionVarsDefs = true;
m_bAddedProjectionVarsData = true;
SetSpatialRefNoUpdate(&oSRS);
}
CPLFree(pszTempProjection);

CPLDebug("netCDF",
"Assummed Longitude Latitude CRS 'OGC:CRS84' because "
"none otherwise available and geotransform within suitable bounds. "
"Set GDAL_NETCDF_ASSUME_LONGLAT=NO as configuration option or "
" ASSUME_LONGLAT=NO as open option to bypass this assumption.");
}
}
}

// Search for Well-known GeogCS if got only CF WKT
// Disabled for now, as a named datum also include control points
// (see mailing list and bug#4281
Expand Down Expand Up @@ -10567,6 +10610,10 @@ void GDALRegister_netCDF()
"description='Whether 2D variables that share the same indexing dimensions "
"should be exposed as several bands of a same dataset instead of several "
"subdatasets.' default='NO'/>"
" <Option name='ASSUME_LONGLAT' type='boolean' scope='raster' "
"description='Whether when all else has failed for determining a CRS, a "
"meaningful geotransform has been found, and is within the "
"bounds -180,360 -90,90, assume OGC:CRS84.' default='NO'/>"
"</OpenOptionList>" );


Expand Down
5 changes: 5 additions & 0 deletions frmts/netcdf/netcdfdataset.h
Original file line number Diff line number Diff line change
Expand Up @@ -90,6 +90,11 @@
/* Config Options

GDAL_NETCDF_BOTTOMUP=YES/NO overrides bottom-up value on import
GDAL_NETCDF_VERIFY_DIMS=[YES/STRICT] : Try to guess which dimensions represent the latitude and longitude only by their attributes (STRICT) or also by guessing the name (YES), default is YES.
GDAL_NETCDF_IGNORE_XY_AXIS_NAME_CHECKS=[YES/NO] Whether X/Y dimensions should be always considered as geospatial axis, even if the lack conventional attributes confirming it. Default is NO.
GDAL_NETCDF_ASSUME_LONGLAT=[YES/NO] Whether when all else has failed for determining a CRS, a meaningful geotransform has been found, and is within the bounds -180,360 -90,90, if YES assume OGC:CRS84. Default is NO.

// TODO: this unusued and a few others occur in the source that are not documented, flush out unused opts and document the rest mdsumner@gmail.com
GDAL_NETCDF_CONVERT_LAT_180=YES/NO convert longitude values from ]180,360] to [-180,180]
*/

Expand Down