Skip to content

Commit a4a9c26

Browse files
rouaultrouault
authored
Merge pull request #13911 from rouault/fix_13885
netCDF: fix writing datasets with non-axis aligned geotransform
2 parents a7ec195 + a68a845 commit a4a9c26

File tree

2 files changed

+170
-61
lines changed

2 files changed

+170
-61
lines changed

‎autotest/gdrivers/netcdf.py‎

Lines changed: 75 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6893,3 +6893,78 @@ def test_netcdf_open_geotransform_gt5_positive():
68936893
[440720.0, 60.0, 0.0, 3751320.0, 0.0, -60.0]
68946894
)
68956895
assert ds.GetRasterBand(1).Checksum() == 4672
6896+
6897+
6898+
###############################################################################
6899+
#
6900+
6901+
6902+
@pytest.mark.parametrize(
6903+
"options,expected_gt,expected_cs,expected_warp_gt,expected_warp_cs",
6904+
[
6905+
(
6906+
[],
6907+
(1841001.75, 1.5, -5.0, 1144003.25, -5.0, -1.5),
6908+
4672,
6909+
(1840901.75, 5.2201532544552744, 0.0, 1144003.25, 0.0, -5.2201532544552744),
6910+
4454,
6911+
),
6912+
(
6913+
["WRITE_LONLAT=YES"],
6914+
(1841001.75, 1.5, -5.0, 1144003.25, -5.0, -1.5),
6915+
4672,
6916+
(
6917+
1840898.6923189342,
6918+
5.245503735406323,
6919+
0.0,
6920+
1144003.2801153609,
6921+
0.0,
6922+
-5.245503735406323,
6923+
),
6924+
4512,
6925+
),
6926+
(
6927+
["WRITE_BOTTOMUP=NO"],
6928+
(1840901.75, 1.5, 5.0, 1143973.25, -5.0, 1.5),
6929+
4855,
6930+
(1840901.75, 5.2201532544552744, 0.0, 1144003.25, 0.0, -5.2201532544552744),
6931+
4454,
6932+
),
6933+
(
6934+
["WRITE_LONLAT=YES", "WRITE_BOTTOMUP=NO"],
6935+
(1840901.75, 1.5, 5.0, 1143973.25, -5.0, 1.5),
6936+
4855,
6937+
(
6938+
1840901.2295569642,
6939+
5.249091373447578,
6940+
0.0,
6941+
1144004.0234124723,
6942+
0.0,
6943+
-5.249091373447578,
6944+
),
6945+
4443,
6946+
),
6947+
],
6948+
)
6949+
def test_netcdf_write_non_axis_aligned_geotransform(
6950+
tmp_path, options, expected_gt, expected_cs, expected_warp_gt, expected_warp_cs
6951+
):
6952+
6953+
src_ds = gdal.Open("../gcore/data/geomatrix.tif")
6954+
gdal.GetDriverByName("netCDF").CreateCopy(
6955+
tmp_path / "out.nc", src_ds, options=options
6956+
)
6957+
ds = gdal.Open(tmp_path / "out.nc")
6958+
assert ds.GetGeoTransform() == pytest.approx(expected_gt)
6959+
assert ds.GetRasterBand(1).Checksum() == expected_cs
6960+
6961+
if "WRITE_LONLAT=YES" in options:
6962+
warped_ds = gdal.Warp(
6963+
"", ds, options="-of MEM -to METHOD=GEOLOC_ARRAY -t_srs EPSG:32611"
6964+
)
6965+
else:
6966+
warped_ds = gdal.Warp(
6967+
"", ds, options="-of MEM -to METHOD=GEOTRANSFORM -t_srs EPSG:32611"
6968+
)
6969+
assert warped_ds.GetGeoTransform() == pytest.approx(expected_warp_gt)
6970+
assert warped_ds.GetRasterBand(1).Checksum() == expected_warp_cs

‎frmts/netcdf/netcdfdataset.cpp‎

Lines changed: 95 additions & 61 deletions
Original file line numberDiff line numberDiff line change
@@ -5620,31 +5620,26 @@ CPLErr netCDFDataset::AddProjectionVars(bool bDefsOnly,
56205620
return CE_Failure;
56215621

56225622
// Optional GDAL custom projection tags.
5623-
if (bWriteGDALTags)
5623+
if (bWriteGDALTags && bWriteGeoTransform && m_bHasGeoTransform)
56245624
{
5625-
std::string osGeoTransform = m_gt.ToString(" ");
5625+
GDALGeoTransform gt(m_gt);
5626+
if (!bBottomUp)
5627+
{
5628+
// Change origin from top to bottom and sign of coefficients
5629+
// indexed by row
5630+
gt.yorig += nRasterYSize * gt.yscale;
5631+
gt.xorig += nRasterYSize * gt.xrot;
5632+
gt.xrot = -gt.xrot;
5633+
gt.yscale = -gt.yscale;
5634+
}
5635+
std::string osGeoTransform = gt.ToString(" ");
56265636
CPLDebug("GDAL_netCDF", "szGeoTransform = %s",
56275637
osGeoTransform.c_str());
56285638

5629-
// if( strlen(pszProj4Defn) > 0 ) {
5630-
// nc_put_att_text(cdfid, NCDFVarID, "proj4",
5631-
// strlen(pszProj4Defn), pszProj4Defn);
5632-
// }
5633-
5634-
// For now, write the geotransform for back-compat or else
5635-
// the old (1.8.1) driver overrides the CF geotransform with
5636-
// empty values from dfNN, dfSN, dfEE, dfWE;
5637-
5638-
// TODO: fix this in 1.8 branch, and then remove this here.
5639-
if (bWriteGeoTransform && m_bHasGeoTransform)
5640-
{
5641-
{
5642-
const int status = nc_put_att_text(
5643-
cdfid, NCDFVarID, NCDF_GEOTRANSFORM,
5644-
osGeoTransform.size(), osGeoTransform.c_str());
5645-
NCDF_ERR(status);
5646-
}
5647-
}
5639+
const int status = nc_put_att_text(
5640+
cdfid, NCDFVarID, NCDF_GEOTRANSFORM, osGeoTransform.size(),
5641+
osGeoTransform.c_str());
5642+
NCDF_ERR(status);
56485643
}
56495644

56505645
// Write projection variable to band variable.
@@ -5657,7 +5652,13 @@ CPLErr netCDFDataset::AddProjectionVars(bool bDefsOnly,
56575652
const bool bIsRotatedPole =
56585653
pszCFProjection != nullptr &&
56595654
EQUAL(pszCFProjection, ROTATED_POLE_VAR_NAME);
5660-
if (bIsRotatedPole)
5655+
5656+
if (m_bHasGeoTransform && !m_gt.IsAxisAligned())
5657+
{
5658+
// Do not write X/Y coordinate arrays
5659+
}
5660+
5661+
else if (bIsRotatedPole)
56615662
{
56625663
// Rename dims to rlat/rlon.
56635664
papszDimName
@@ -5873,49 +5874,55 @@ CPLErr netCDFDataset::AddProjectionVars(bool bDefsOnly,
58735874
return CE_Failure;
58745875
}
58755876

5876-
// Get Y values.
5877-
const double dfY0 = (!bBottomUp) ? m_gt.yorig :
5878-
// Invert latitude values.
5879-
m_gt.yorig + (m_gt.yscale * nRasterYSize);
5880-
const double dfDY = m_gt.yscale;
5877+
// Make sure we are in data mode.
5878+
SetDefineMode(false);
58815879

5882-
for (int j = 0; j < nRasterYSize; j++)
5880+
int status = NC_NOERR;
5881+
5882+
if (m_gt.IsAxisAligned())
58835883
{
5884-
// The data point is centered inside the pixel.
5885-
if (!bBottomUp)
5886-
padYVal[j] = dfY0 + (j + 0.5) * dfDY;
5887-
else // Invert latitude values.
5888-
padYVal[j] = dfY0 - (j + 0.5) * dfDY;
5889-
}
5890-
startX[0] = 0;
5891-
countX[0] = nRasterXSize;
5884+
// Get Y values.
5885+
const double dfY0 =
5886+
(!bBottomUp) ? m_gt.yorig :
5887+
// Invert latitude values.
5888+
m_gt.yorig + (m_gt.yscale * nRasterYSize);
5889+
const double dfDY = m_gt.yscale;
58925890

5893-
// Get X values.
5894-
const double dfX0 = m_gt.xorig;
5895-
const double dfDX = m_gt.xscale;
5891+
for (int j = 0; j < nRasterYSize; j++)
5892+
{
5893+
// The data point is centered inside the pixel.
5894+
if (!bBottomUp)
5895+
padYVal[j] = dfY0 + (j + 0.5) * dfDY;
5896+
else // Invert latitude values.
5897+
padYVal[j] = dfY0 - (j + 0.5) * dfDY;
5898+
}
5899+
startX[0] = 0;
5900+
countX[0] = nRasterXSize;
58965901

5897-
for (int i = 0; i < nRasterXSize; i++)
5898-
{
5899-
// The data point is centered inside the pixel.
5900-
padXVal[i] = dfX0 + (i + 0.5) * dfDX;
5901-
}
5902-
startY[0] = 0;
5903-
countY[0] = nRasterYSize;
5902+
// Get X values.
5903+
const double dfX0 = m_gt.xorig;
5904+
const double dfDX = m_gt.xscale;
59045905

5905-
// Write X/Y values.
5906+
for (int i = 0; i < nRasterXSize; i++)
5907+
{
5908+
// The data point is centered inside the pixel.
5909+
padXVal[i] = dfX0 + (i + 0.5) * dfDX;
5910+
}
5911+
startY[0] = 0;
5912+
countY[0] = nRasterYSize;
59065913

5907-
// Make sure we are in data mode.
5908-
SetDefineMode(false);
5914+
// Write X/Y values.
59095915

5910-
CPLDebug("GDAL_netCDF", "Writing X values");
5911-
int status =
5912-
nc_put_vara_double(cdfid, nVarXID, startX, countX, padXVal);
5913-
NCDF_ERR(status);
5916+
CPLDebug("GDAL_netCDF", "Writing X values");
5917+
status =
5918+
nc_put_vara_double(cdfid, nVarXID, startX, countX, padXVal);
5919+
NCDF_ERR(status);
59145920

5915-
CPLDebug("GDAL_netCDF", "Writing Y values");
5916-
status =
5917-
nc_put_vara_double(cdfid, nVarYID, startY, countY, padYVal);
5918-
NCDF_ERR(status);
5921+
CPLDebug("GDAL_netCDF", "Writing Y values");
5922+
status =
5923+
nc_put_vara_double(cdfid, nVarYID, startY, countY, padYVal);
5924+
NCDF_ERR(status);
5925+
}
59195926

59205927
if (pfnProgress)
59215928
pfnProgress(0.20, nullptr, pProgressData);
@@ -5979,10 +5986,37 @@ CPLErr netCDFDataset::AddProjectionVars(bool bDefsOnly,
59795986
if (!bHasGeoloc)
59805987
{
59815988
// Fill values to transform.
5982-
for (int i = 0; i < nRasterXSize; i++)
5989+
if (m_gt.IsAxisAligned())
5990+
{
5991+
for (int i = 0; i < nRasterXSize; i++)
5992+
{
5993+
padLatVal[i] = padYVal[j];
5994+
padLonVal[i] = padXVal[i];
5995+
}
5996+
}
5997+
else
59835998
{
5984-
padLatVal[i] = padYVal[j];
5985-
padLonVal[i] = padXVal[i];
5999+
for (int i = 0; i < nRasterXSize; i++)
6000+
{
6001+
if (!bBottomUp)
6002+
{
6003+
padLatVal[i] = m_gt.yorig +
6004+
(i + 0.5) * m_gt.yrot +
6005+
(j + 0.5) * m_gt.yscale;
6006+
padLonVal[i] = m_gt.xorig +
6007+
(i + 0.5) * m_gt.xscale +
6008+
(j + 0.5) * m_gt.xrot;
6009+
}
6010+
else
6011+
{
6012+
padLatVal[i] =
6013+
m_gt.yorig + (i + 0.5) * m_gt.yrot +
6014+
(nRasterYSize - j - 0.5) * m_gt.yscale;
6015+
padLonVal[i] =
6016+
m_gt.xorig + (i + 0.5) * m_gt.xscale +
6017+
(nRasterYSize - j - 0.5) * m_gt.xrot;
6018+
}
6019+
}
59866020
}
59876021

59886022
// Do the transform.
@@ -6041,7 +6075,7 @@ CPLErr netCDFDataset::AddProjectionVars(bool bDefsOnly,
60416075
} // Projected
60426076

60436077
// If not projected/geographic and has geoloc
6044-
else if (!bIsGeographic && bHasGeoloc)
6078+
else if (!bIsGeographic && bHasGeoloc && m_gt.IsAxisAligned())
60456079
{
60466080
// Use
60476081
// https://cfconventions.org/Data/cf-conventions/cf-conventions-1.9/cf-conventions.html#_two_dimensional_latitude_longitude_coordinate_variables

0 commit comments

Comments
 (0)