Skip to content

Commit 446d5ee

Browse files
rouaultrouault
authored
Merge pull request OSGeo#12454 from rouault/hillshade_improvements
gdaldem/gdal raster hillshade/slope/aspect/roughness/tpi/tri streaming: expose overviews if source has overviews
2 parents 3eab681 + f64ed06 commit 446d5ee

11 files changed

+508
-153
lines changed

‎apps/gdal_translate_lib.cpp‎

Lines changed: 3 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1566,12 +1566,9 @@ GDALDatasetH GDALTranslate(const char *pszDest, GDALDatasetH hSrcDataset,
15661566
psOptions->srcWin.dfXOff * adfDstGeoTransform[4] +
15671567
psOptions->srcWin.dfYOff * adfDstGeoTransform[5];
15681568

1569-
const double dfX = static_cast<double>(nOXSize);
1570-
const double dfY = static_cast<double>(nOYSize);
1571-
adfDstGeoTransform[1] *= psOptions->srcWin.dfXSize / dfX;
1572-
adfDstGeoTransform[2] *= psOptions->srcWin.dfYSize / dfY;
1573-
adfDstGeoTransform[4] *= psOptions->srcWin.dfXSize / dfX;
1574-
adfDstGeoTransform[5] *= psOptions->srcWin.dfYSize / dfY;
1569+
const double dfXRatio = psOptions->srcWin.dfXSize / nOXSize;
1570+
const double dfYRatio = psOptions->srcWin.dfYSize / nOYSize;
1571+
GDALRescaleGeoTransform(adfDstGeoTransform, dfXRatio, dfYRatio);
15751572

15761573
if (psOptions->dfXRes != 0.0)
15771574
{

‎apps/gdaldem_lib.cpp‎

Lines changed: 291 additions & 138 deletions
Large diffs are not rendered by default.

‎autotest/utilities/test_gdalalg_raster_aspect.py‎

Lines changed: 45 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -123,3 +123,48 @@ def test_gdalalg_raster_aspect_vrt_output_pipeline_from_filename():
123123
match=r"aspect: VRT output is not supported. Consider using the GDALG driver instead \(files with \.gdalg\.json extension\)",
124124
):
125125
alg.Run()
126+
127+
128+
def test_gdalalg_raster_aspect_overview():
129+
130+
src_ds = gdal.Translate(
131+
"",
132+
"../gdrivers/data/n43.tif",
133+
format="MEM",
134+
width=61,
135+
resampleAlg=gdal.GRIORA_Bilinear,
136+
)
137+
138+
with gdal.Run(
139+
"raster",
140+
"aspect",
141+
input=src_ds,
142+
output="",
143+
output_format="stream",
144+
) as alg:
145+
out_ds = alg.Output()
146+
cs = out_ds.GetRasterBand(1).Checksum()
147+
stats = out_ds.GetRasterBand(1).ComputeStatistics(False)
148+
149+
src_ds = gdal.Translate("", "../gdrivers/data/n43.tif", format="MEM")
150+
src_ds.BuildOverviews("BILINEAR", [2, 4])
151+
152+
with gdal.Run(
153+
"raster",
154+
"aspect",
155+
input=src_ds,
156+
output="",
157+
output_format="stream",
158+
) as alg:
159+
out_ds = alg.Output()
160+
del src_ds
161+
162+
assert out_ds.GetRasterBand(1).GetOverviewCount() == 2
163+
assert out_ds.GetRasterBand(1).GetOverview(-1) is None
164+
assert out_ds.GetRasterBand(1).GetOverview(2) is None
165+
assert out_ds.GetRasterBand(1).GetOverview(0).XSize == 61
166+
assert out_ds.GetRasterBand(1).GetOverview(0).YSize == 61
167+
assert out_ds.GetRasterBand(1).GetOverview(1).XSize == 31
168+
assert out_ds.GetRasterBand(1).GetOverview(1).YSize == 31
169+
assert out_ds.GetRasterBand(1).GetOverview(0).Checksum() == cs
170+
assert out_ds.GetRasterBand(1).GetOverview(0).ComputeStatistics(False) == stats

‎autotest/utilities/test_gdalalg_raster_hillshade.py‎

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -154,3 +154,40 @@ def test_gdalalg_raster_hillshade_vrt_output_pipeline_from_filename():
154154
match=r"hillshade: VRT output is not supported. Consider using the GDALG driver instead \(files with \.gdalg\.json extension\)",
155155
):
156156
alg.Run()
157+
158+
159+
def test_gdalalg_raster_hillshade_overview():
160+
161+
src_ds = gdal.Translate(
162+
"",
163+
"../gdrivers/data/n43.tif",
164+
format="MEM",
165+
width=61,
166+
resampleAlg=gdal.GRIORA_Bilinear,
167+
)
168+
169+
with gdal.Run(
170+
"raster", "hillshade", input=src_ds, output="", output_format="stream", z=30
171+
) as alg:
172+
out_ds = alg.Output()
173+
cs = out_ds.GetRasterBand(1).Checksum()
174+
stats = out_ds.GetRasterBand(1).ComputeStatistics(False)
175+
176+
src_ds = gdal.Translate("", "../gdrivers/data/n43.tif", format="MEM")
177+
src_ds.BuildOverviews("BILINEAR", [2, 4])
178+
179+
with gdal.Run(
180+
"raster", "hillshade", input=src_ds, output="", output_format="stream", z=30
181+
) as alg:
182+
out_ds = alg.Output()
183+
del src_ds
184+
185+
assert out_ds.GetRasterBand(1).GetOverviewCount() == 2
186+
assert out_ds.GetRasterBand(1).GetOverview(-1) is None
187+
assert out_ds.GetRasterBand(1).GetOverview(2) is None
188+
assert out_ds.GetRasterBand(1).GetOverview(0).XSize == 61
189+
assert out_ds.GetRasterBand(1).GetOverview(0).YSize == 61
190+
assert out_ds.GetRasterBand(1).GetOverview(1).XSize == 31
191+
assert out_ds.GetRasterBand(1).GetOverview(1).YSize == 31
192+
assert out_ds.GetRasterBand(1).GetOverview(0).Checksum() == cs
193+
assert out_ds.GetRasterBand(1).GetOverview(0).ComputeStatistics(False) == stats

‎autotest/utilities/test_gdalalg_raster_roughness.py‎

Lines changed: 45 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -119,3 +119,48 @@ def test_gdalalg_raster_roughness_vrt_output_pipeline_from_filename():
119119
match=r"roughness: VRT output is not supported. Consider using the GDALG driver instead \(files with \.gdalg\.json extension\)",
120120
):
121121
alg.Run()
122+
123+
124+
def test_gdalalg_raster_roughness_overview():
125+
126+
src_ds = gdal.Translate(
127+
"",
128+
"../gdrivers/data/n43.tif",
129+
format="MEM",
130+
width=61,
131+
resampleAlg=gdal.GRIORA_Bilinear,
132+
)
133+
134+
with gdal.Run(
135+
"raster",
136+
"roughness",
137+
input=src_ds,
138+
output="",
139+
output_format="stream",
140+
) as alg:
141+
out_ds = alg.Output()
142+
cs = out_ds.GetRasterBand(1).Checksum()
143+
stats = out_ds.GetRasterBand(1).ComputeStatistics(False)
144+
145+
src_ds = gdal.Translate("", "../gdrivers/data/n43.tif", format="MEM")
146+
src_ds.BuildOverviews("BILINEAR", [2, 4])
147+
148+
with gdal.Run(
149+
"raster",
150+
"roughness",
151+
input=src_ds,
152+
output="",
153+
output_format="stream",
154+
) as alg:
155+
out_ds = alg.Output()
156+
del src_ds
157+
158+
assert out_ds.GetRasterBand(1).GetOverviewCount() == 2
159+
assert out_ds.GetRasterBand(1).GetOverview(-1) is None
160+
assert out_ds.GetRasterBand(1).GetOverview(2) is None
161+
assert out_ds.GetRasterBand(1).GetOverview(0).XSize == 61
162+
assert out_ds.GetRasterBand(1).GetOverview(0).YSize == 61
163+
assert out_ds.GetRasterBand(1).GetOverview(1).XSize == 31
164+
assert out_ds.GetRasterBand(1).GetOverview(1).YSize == 31
165+
assert out_ds.GetRasterBand(1).GetOverview(0).Checksum() == cs
166+
assert out_ds.GetRasterBand(1).GetOverview(0).ComputeStatistics(False) == stats

‎autotest/utilities/test_gdalalg_raster_slope.py‎

Lines changed: 45 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -124,3 +124,48 @@ def test_gdalalg_raster_slope_vrt_output_pipeline_from_filename():
124124
match=r"slope: VRT output is not supported. Consider using the GDALG driver instead \(files with \.gdalg\.json extension\)",
125125
):
126126
alg.Run()
127+
128+
129+
def test_gdalalg_raster_slope_overview():
130+
131+
src_ds = gdal.Translate(
132+
"",
133+
"../gdrivers/data/n43.tif",
134+
format="MEM",
135+
width=61,
136+
resampleAlg=gdal.GRIORA_Bilinear,
137+
)
138+
139+
with gdal.Run(
140+
"raster",
141+
"slope",
142+
input=src_ds,
143+
output="",
144+
output_format="stream",
145+
) as alg:
146+
out_ds = alg.Output()
147+
cs = out_ds.GetRasterBand(1).Checksum()
148+
stats = out_ds.GetRasterBand(1).ComputeStatistics(False)
149+
150+
src_ds = gdal.Translate("", "../gdrivers/data/n43.tif", format="MEM")
151+
src_ds.BuildOverviews("BILINEAR", [2, 4])
152+
153+
with gdal.Run(
154+
"raster",
155+
"slope",
156+
input=src_ds,
157+
output="",
158+
output_format="stream",
159+
) as alg:
160+
out_ds = alg.Output()
161+
del src_ds
162+
163+
assert out_ds.GetRasterBand(1).GetOverviewCount() == 2
164+
assert out_ds.GetRasterBand(1).GetOverview(-1) is None
165+
assert out_ds.GetRasterBand(1).GetOverview(2) is None
166+
assert out_ds.GetRasterBand(1).GetOverview(0).XSize == 61
167+
assert out_ds.GetRasterBand(1).GetOverview(0).YSize == 61
168+
assert out_ds.GetRasterBand(1).GetOverview(1).XSize == 31
169+
assert out_ds.GetRasterBand(1).GetOverview(1).YSize == 31
170+
assert out_ds.GetRasterBand(1).GetOverview(0).Checksum() == cs
171+
assert out_ds.GetRasterBand(1).GetOverview(0).ComputeStatistics(False) == stats

‎frmts/mem/memdataset.cpp‎

Lines changed: 12 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -888,6 +888,16 @@ CPLErr MEMDataset::IBuildOverviews(const char *pszResampling, int nOverviews,
888888
(nRasterXSize + panOverviewList[i] - 1) / panOverviewList[i];
889889
poOvrDS->nRasterYSize =
890890
(nRasterYSize + panOverviewList[i] - 1) / panOverviewList[i];
891+
poOvrDS->bGeoTransformSet = bGeoTransformSet;
892+
memcpy(poOvrDS->adfGeoTransform, adfGeoTransform,
893+
6 * sizeof(double));
894+
const double dfOvrXRatio =
895+
static_cast<double>(nRasterXSize) / poOvrDS->nRasterXSize;
896+
const double dfOvrYRatio =
897+
static_cast<double>(nRasterYSize) / poOvrDS->nRasterYSize;
898+
GDALRescaleGeoTransform(poOvrDS->adfGeoTransform, dfOvrXRatio,
899+
dfOvrYRatio);
900+
poOvrDS->m_oSRS = m_oSRS;
891901
for (int iBand = 0; iBand < nBands; iBand++)
892902
{
893903
const GDALDataType eDT =
@@ -897,7 +907,7 @@ CPLErr MEMDataset::IBuildOverviews(const char *pszResampling, int nOverviews,
897907
return CE_Failure;
898908
}
899909
}
900-
m_apoOverviewDS.emplace_back(std::move(poOvrDS));
910+
m_apoOverviewDS.emplace_back(poOvrDS.release());
901911
}
902912
}
903913

@@ -1079,7 +1089,7 @@ std::unique_ptr<GDALDataset> MEMDataset::Clone(int nScopeFlags,
10791089
for (const auto &poOvrDS : m_apoOverviewDS)
10801090
{
10811091
poNewDS->m_apoOverviewDS.emplace_back(
1082-
poOvrDS->Clone(nScopeFlags, bCanShareState));
1092+
poOvrDS->Clone(nScopeFlags, bCanShareState).release());
10831093
}
10841094

10851095
poNewDS->SetDescription(GetDescription());

‎frmts/mem/memdataset.h‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -52,7 +52,10 @@ class CPL_DLL MEMDataset CPL_NON_FINAL : public GDALDataset
5252
std::vector<gdal::GCP> m_aoGCPs{};
5353
OGRSpatialReference m_oGCPSRS{};
5454

55-
std::vector<std::unique_ptr<GDALDataset>> m_apoOverviewDS{};
55+
using GDALDatasetRefCountedPtr =
56+
std::unique_ptr<GDALDataset, GDALDatasetUniquePtrReleaser>;
57+
58+
std::vector<GDALDatasetRefCountedPtr> m_apoOverviewDS{};
5659

5760
struct Private;
5861
std::unique_ptr<Private> m_poPrivate;

‎gcore/gdal_misc.cpp‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5842,3 +5842,23 @@ bool GDALDoesFileOrDatasetExist(const char *pszName, const char **ppszType,
58425842

58435843
return false;
58445844
}
5845+
5846+
/************************************************************************/
5847+
/* GDALRescaleGeoTransform() */
5848+
/************************************************************************/
5849+
5850+
/** Rescale a geotransform by multiplying its scale and rotation terms by
5851+
* the provided ratios.
5852+
*
5853+
* This is typically used to compute the geotransform matrix of an overview
5854+
* dataset from the full resolution dataset, where the ratios are the size
5855+
* of the full resolution dataset divided by the size of the overview.
5856+
*/
5857+
void GDALRescaleGeoTransform(double adfGeoTransform[6], double dfXRatio,
5858+
double dfYRatio)
5859+
{
5860+
adfGeoTransform[1] *= dfXRatio;
5861+
adfGeoTransform[2] *= dfYRatio;
5862+
adfGeoTransform[4] *= dfXRatio;
5863+
adfGeoTransform[5] *= dfYRatio;
5864+
}

‎gcore/gdal_priv.h‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4774,6 +4774,9 @@ struct GDALColorAssociation
47744774
std::vector<GDALColorAssociation> GDALLoadTextColorMap(const char *pszFilename,
47754775
GDALRasterBand *poBand);
47764776

4777+
void GDALRescaleGeoTransform(double adfGeoTransform[6], double dfXRatio,
4778+
double dfYRatio);
4779+
47774780
// Macro used so that Identify and driver metadata methods in drivers built
47784781
// as plugin can be duplicated in libgdal core and in the driver under different
47794782
// names

0 commit comments

Comments
 (0)