Skip to content

Commit

Permalink
Merge pull request #10884 from OSGeo/backport-10812-to-release/3.9
Browse files Browse the repository at this point in the history 
[Backport release/3.9] Speed-up geolocation-array based warping of Sentinel3 datasets
  • Loading branch information
@rouault
rouault committed Sep 27, 2024
2 parents ae7e834 + 5b42a3f commit 664f137
Show file tree
Hide file tree
Showing 5 changed files with 185 additions and 61 deletions.
14 changes: 10 additions & 4 deletions alg/gdalgeoloc.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2006,14 +2006,20 @@ void *GDALCreateGeoLocTransformerEx(GDALDatasetH hBaseDS,
psTransform->bUseArray = !CPLTestBool(pszUseTempDatasets);
else
{
psTransform->bUseArray = nXSize < 16 * 1000 * 1000 / nYSize;
if (psTransform->bUseArray)
constexpr int MEGAPIXEL_LIMIT = 24;
psTransform->bUseArray =
nXSize < MEGAPIXEL_LIMIT * 1000 * 1000 / nYSize;
if (!psTransform->bUseArray)
{
CPLDebug("GEOLOC",
"Using temporary GTiff backing to store backmap, because "
"geoloc arrays exceed 16 megapixels. You can set the "
"geoloc arrays require %d megapixels, exceeding the %d "
"megapixels limit. You can set the "
"GDAL_GEOLOC_USE_TEMP_DATASETS configuration option to "
"NO to force RAM storage of backmap");
"NO to force RAM storage of backmap",
static_cast<int>(static_cast<int64_t>(nXSize) * nYSize /
(1000 * 1000)),
MEGAPIXEL_LIMIT);
}
}

Expand Down
15 changes: 8 additions & 7 deletions alg/gdalgeoloc_dataset_accessor.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,13 +53,14 @@ class GDALGeoLocDatasetAccessors
bool LoadGeoloc(bool bIsRegularGrid);

public:
static constexpr int TILE_SIZE = 1024;

GDALCachedPixelAccessor<double, TILE_SIZE> geolocXAccessor;
GDALCachedPixelAccessor<double, TILE_SIZE> geolocYAccessor;
GDALCachedPixelAccessor<float, TILE_SIZE> backMapXAccessor;
GDALCachedPixelAccessor<float, TILE_SIZE> backMapYAccessor;
GDALCachedPixelAccessor<float, TILE_SIZE> backMapWeightAccessor;
static constexpr int TILE_SIZE = 256;
static constexpr int TILE_COUNT = 64;

GDALCachedPixelAccessor<double, TILE_SIZE, TILE_COUNT> geolocXAccessor;
GDALCachedPixelAccessor<double, TILE_SIZE, TILE_COUNT> geolocYAccessor;
GDALCachedPixelAccessor<float, TILE_SIZE, TILE_COUNT> backMapXAccessor;
GDALCachedPixelAccessor<float, TILE_SIZE, TILE_COUNT> backMapYAccessor;
GDALCachedPixelAccessor<float, TILE_SIZE, TILE_COUNT> backMapWeightAccessor;

explicit GDALGeoLocDatasetAccessors(GDALGeoLocTransformInfo *psTransform)
: m_psTransform(psTransform), geolocXAccessor(nullptr),
Expand Down
106 changes: 68 additions & 38 deletions autotest/gcore/vrt_read.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2451,30 +2451,68 @@ def test_vrt_read_top_and_bottom_strips_average():


@pytest.mark.parametrize(
"input_datatype", [gdal.GDT_Byte, gdal.GDT_UInt16, gdal.GDT_Int16]
)
@pytest.mark.parametrize("vrt_type", ["Byte", "UInt16", "Int16"])
@pytest.mark.parametrize("nodata", [0, 254])
@pytest.mark.parametrize(
"request_type", [gdal.GDT_Byte, gdal.GDT_UInt16, gdal.GDT_Int16]
"input_datatype,vrt_type,nodata,vrt_nodata,request_type",
[
(gdal.GDT_Byte, "Byte", 0, 255, gdal.GDT_Byte),
(gdal.GDT_Byte, "Byte", 254, 255, gdal.GDT_Byte),
(gdal.GDT_Byte, "Int8", 254, 255, gdal.GDT_Byte),
(gdal.GDT_Byte, "Byte", 254, 127, gdal.GDT_Int8),
(gdal.GDT_Byte, "UInt16", 254, 255, gdal.GDT_Byte),
(gdal.GDT_Byte, "Byte", 254, 255, gdal.GDT_UInt16),
(gdal.GDT_Int8, "Int8", 0, 127, gdal.GDT_Int8),
(gdal.GDT_Int8, "Int16", 0, 127, gdal.GDT_Int8),
(gdal.GDT_UInt16, "UInt16", 0, 65535, gdal.GDT_UInt16),
(gdal.GDT_Int16, "Int16", 0, 32767, gdal.GDT_Int16),
(gdal.GDT_UInt32, "UInt32", 0, (1 << 31) - 1, gdal.GDT_UInt32),
(gdal.GDT_Int32, "Int32", 0, (1 << 30) - 1, gdal.GDT_Int32),
(gdal.GDT_Int32, "Float32", 0, (1 << 30) - 1, gdal.GDT_Float64),
(gdal.GDT_UInt64, "UInt64", 0, (1 << 63) - 1, gdal.GDT_UInt64),
(gdal.GDT_Int64, "Int64", 0, (1 << 62) - 1, gdal.GDT_Int64),
(gdal.GDT_Int64, "Float32", 0, (1 << 62), gdal.GDT_Int64),
(gdal.GDT_Float32, "Float32", 0, 1.5, gdal.GDT_Float32),
(gdal.GDT_Float32, "Float32", 0, 1.5, gdal.GDT_Float64),
(gdal.GDT_Float32, "Float64", 0, 1.5, gdal.GDT_Float32),
(gdal.GDT_Float32, "Float64", 0, 1.5, gdal.GDT_Float64),
(gdal.GDT_Float64, "Float64", 0, 1.5, gdal.GDT_Float64),
(gdal.GDT_Float64, "Float32", 0, 1.5, gdal.GDT_Float64),
],
)
def test_vrt_read_complex_source_nodata(
tmp_vsimem, input_datatype, vrt_type, nodata, request_type
tmp_vsimem, input_datatype, vrt_type, nodata, vrt_nodata, request_type
):

if input_datatype == gdal.GDT_Byte:
array_type = "B"
elif input_datatype == gdal.GDT_UInt16:
array_type = "H"
elif input_datatype == gdal.GDT_Int16:
array_type = "h"
def get_array_type(dt):
m = {
gdal.GDT_Byte: "B",
gdal.GDT_Int8: "b",
gdal.GDT_UInt16: "H",
gdal.GDT_Int16: "h",
gdal.GDT_UInt32: "I",
gdal.GDT_Int32: "i",
gdal.GDT_UInt64: "Q",
gdal.GDT_Int64: "q",
gdal.GDT_Float32: "f",
gdal.GDT_Float64: "d",
}
return m[dt]

if input_datatype in (gdal.GDT_Float32, gdal.GDT_Float64):
input_val = 1.75
if vrt_type in ("Float32", "Float64") and request_type in (
gdal.GDT_Float32,
gdal.GDT_Float64,
):
expected_val = input_val
else:
expected_val = math.round(input_val)
else:
assert False
input_val = 1
expected_val = input_val

input_data = array.array(
array_type,
get_array_type(input_datatype),
[
nodata,
1,
input_val,
2,
3,
nodata, # EOL
Expand Down Expand Up @@ -2511,7 +2549,7 @@ def test_vrt_read_complex_source_nodata(
ds.Close()
complex_xml = f"""<VRTDataset rasterXSize="5" rasterYSize="6">
<VRTRasterBand dataType="{vrt_type}" band="1">
<NoDataValue>255</NoDataValue>
<NoDataValue>{vrt_nodata}</NoDataValue>
<ComplexSource>
<SourceFilename relativeToVRT="0">{input_filename}</SourceFilename>
<SourceBand>1</SourceBand>
Expand All @@ -2521,24 +2559,16 @@ def test_vrt_read_complex_source_nodata(
</VRTDataset>
"""
vrt_ds = gdal.Open(complex_xml)
if request_type == gdal.GDT_Byte:
array_request_type = "B"
elif request_type == gdal.GDT_UInt16:
array_request_type = "H"
elif request_type == gdal.GDT_Int16:
array_request_type = "h"
else:
assert False
got_data = vrt_ds.ReadRaster(buf_type=request_type)
got_data = struct.unpack(array_request_type * (5 * 6), got_data)
got_data = struct.unpack(get_array_type(request_type) * (5 * 6), got_data)
assert got_data == (
255,
1,
vrt_nodata,
expected_val,
2,
3,
255, # EOL
vrt_nodata, # EOL
4,
255,
vrt_nodata,
5,
6,
7, # EOL
Expand All @@ -2549,18 +2579,18 @@ def test_vrt_read_complex_source_nodata(
20, # EOL
8,
9,
255,
vrt_nodata,
10,
11, # EOL
255,
255,
255,
255,
255, # EOL
vrt_nodata,
vrt_nodata,
vrt_nodata,
vrt_nodata,
vrt_nodata, # EOL
12,
13,
14,
255,
vrt_nodata,
15, # EOL
)

Expand Down
24 changes: 23 additions & 1 deletion frmts/vrt/vrtsourcedrasterband.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -320,7 +320,9 @@ CPLErr VRTSourcedRasterBand::IRasterIO(
// Do nothing
}
else if (nPixelSpace == GDALGetDataTypeSizeBytes(eBufType) &&
(!m_bNoDataValueSet || m_dfNoDataValue == 0.0))
!(m_bNoDataValueSet && m_dfNoDataValue != 0.0) &&
!(m_bNoDataSetAsInt64 && m_nNoDataValueInt64 != 0) &&
!(m_bNoDataSetAsUInt64 && m_nNoDataValueUInt64 != 0))
{
if (nLineSpace == nBufXSize * nPixelSpace)
{
Expand All @@ -336,6 +338,26 @@ CPLErr VRTSourcedRasterBand::IRasterIO(
}
}
}
else if (m_bNoDataSetAsInt64)
{
for (int iLine = 0; iLine < nBufYSize; iLine++)
{
GDALCopyWords(&m_nNoDataValueInt64, GDT_Int64, 0,
static_cast<GByte *>(pData) +
static_cast<GIntBig>(nLineSpace) * iLine,
eBufType, static_cast<int>(nPixelSpace), nBufXSize);
}
}
else if (m_bNoDataSetAsUInt64)
{
for (int iLine = 0; iLine < nBufYSize; iLine++)
{
GDALCopyWords(&m_nNoDataValueUInt64, GDT_UInt64, 0,
static_cast<GByte *>(pData) +
static_cast<GIntBig>(nLineSpace) * iLine,
eBufType, static_cast<int>(nPixelSpace), nBufXSize);
}
}
else
{
double dfWriteValue = 0.0;
Expand Down
87 changes: 76 additions & 11 deletions frmts/vrt/vrtsources.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2860,10 +2860,10 @@ CPLErr VRTComplexSource::RasterIO(GDALDataType eVRTBandDataType, int nXOff,
GByte *const pabyOut = static_cast<GByte *>(pData) +
nPixelSpace * nOutXOff +
static_cast<GPtrDiff_t>(nLineSpace) * nOutYOff;
const auto eSourceType = poSourceBand->GetRasterDataType();
if (m_nProcessingFlags == PROCESSING_FLAG_NODATA)
{
// Optimization if doing only nodata processing
const auto eSourceType = poSourceBand->GetRasterDataType();
if (eSourceType == GDT_Byte)
{
if (!GDALIsValueInRange<GByte>(m_dfNoDataValue))
Expand Down Expand Up @@ -2917,7 +2917,10 @@ CPLErr VRTComplexSource::RasterIO(GDALDataType eVRTBandDataType, int nXOff,
// For Int32, float32 isn't sufficiently precise as working data type
if (eVRTBandDataType == GDT_CInt32 || eVRTBandDataType == GDT_CFloat64 ||
eVRTBandDataType == GDT_Int32 || eVRTBandDataType == GDT_UInt32 ||
eVRTBandDataType == GDT_Float64)
eVRTBandDataType == GDT_Int64 || eVRTBandDataType == GDT_UInt64 ||
eVRTBandDataType == GDT_Float64 || eSourceType == GDT_Int32 ||
eSourceType == GDT_UInt32 || eSourceType == GDT_Int64 ||
eSourceType == GDT_UInt64)
{
eErr = RasterIOInternal<double>(
poSourceBand, eVRTBandDataType, nReqXOff, nReqYOff, nReqXSize,
Expand Down Expand Up @@ -2948,6 +2951,52 @@ static inline bool hasZeroByte(uint32_t v)
return (((v)-0x01010101U) & ~(v)&0x80808080U) != 0;
}

/************************************************************************/
/* CopyWord() */
/************************************************************************/

template <class SrcType>
static void CopyWord(const SrcType *pSrcVal, GDALDataType eSrcType,
void *pDstVal, GDALDataType eDstType)
{
switch (eDstType)
{
case GDT_Byte:
GDALCopyWord(*pSrcVal, *static_cast<uint8_t *>(pDstVal));
break;
case GDT_Int8:
GDALCopyWord(*pSrcVal, *static_cast<int8_t *>(pDstVal));
break;
case GDT_UInt16:
GDALCopyWord(*pSrcVal, *static_cast<uint16_t *>(pDstVal));
break;
case GDT_Int16:
GDALCopyWord(*pSrcVal, *static_cast<int16_t *>(pDstVal));
break;
case GDT_UInt32:
GDALCopyWord(*pSrcVal, *static_cast<uint32_t *>(pDstVal));
break;
case GDT_Int32:
GDALCopyWord(*pSrcVal, *static_cast<int32_t *>(pDstVal));
break;
case GDT_UInt64:
GDALCopyWord(*pSrcVal, *static_cast<uint64_t *>(pDstVal));
break;
case GDT_Int64:
GDALCopyWord(*pSrcVal, *static_cast<int64_t *>(pDstVal));
break;
case GDT_Float32:
GDALCopyWord(*pSrcVal, *static_cast<float *>(pDstVal));
break;
case GDT_Float64:
GDALCopyWord(*pSrcVal, *static_cast<double *>(pDstVal));
break;
default:
GDALCopyWords(pSrcVal, eSrcType, 0, pDstVal, eDstType, 0, 1);
break;
}
}

/************************************************************************/
/* RasterIOProcessNoData() */
/************************************************************************/
Expand Down Expand Up @@ -3103,8 +3152,8 @@ CPLErr VRTComplexSource::RasterIOProcessNoData(
{
if (paSrcData[idxBuffer] != nNoDataValue)
{
GDALCopyWords(&paSrcData[idxBuffer], eSourceType, 0,
pDstLocation, eBufType, 0, 1);
CopyWord(&paSrcData[idxBuffer], eSourceType, pDstLocation,
eBufType);
}
}
}
Expand All @@ -3124,8 +3173,8 @@ CPLErr VRTComplexSource::RasterIOProcessNoData(
{
// Convert first to the VRTRasterBand data type
// to get its clamping, before outputting to buffer data type
GDALCopyWords(&paSrcData[idxBuffer], eSourceType, 0,
abyTemp, eVRTBandDataType, 0, 1);
CopyWord(&paSrcData[idxBuffer], eSourceType, abyTemp,
eVRTBandDataType);
GDALCopyWords(abyTemp, eVRTBandDataType, 0, pDstLocation,
eBufType, 0, 1);
}
Expand Down Expand Up @@ -3283,6 +3332,12 @@ CPLErr VRTComplexSource::RasterIOInternal(
/* Selectively copy into output buffer with nodata masking, */
/* and/or scaling. */
/* -------------------------------------------------------------------- */

const bool bTwoStepDataTypeConversion =
CPL_TO_BOOL(
GDALDataTypeIsConversionLossy(eWrkDataType, eVRTBandDataType)) &&
!CPL_TO_BOOL(GDALDataTypeIsConversionLossy(eVRTBandDataType, eBufType));

size_t idxBuffer = 0;
for (int iY = 0; iY < nOutYSize; iY++)
{
Expand Down Expand Up @@ -3419,18 +3474,28 @@ CPLErr VRTComplexSource::RasterIOInternal(
255.0f,
std::max(0.0f, static_cast<float>(afResult[0]) + 0.5f)));
}
else if (eBufType == eVRTBandDataType)
else if (!bTwoStepDataTypeConversion)
{
CopyWord<WorkingDT>(afResult, eWrkDataType, pDstLocation,
eBufType);
}
else if (eVRTBandDataType == GDT_Float32 && eBufType == GDT_Float64)
{
GDALCopyWords(afResult, eWrkDataType, 0, pDstLocation, eBufType,
0, 1);
// Particular case of the below 2-step conversion.
// Helps a bit for some geolocation based warping with Sentinel3
// data where the longitude/latitude arrays are Int32 bands,
// rescaled in VRT as Float32 and requested as Float64
float fVal;
GDALCopyWord(afResult[0], fVal);
*reinterpret_cast<double *>(pDstLocation) = fVal;
}
else
{
GByte abyTemp[2 * sizeof(double)];
// Convert first to the VRTRasterBand data type
// to get its clamping, before outputting to buffer data type
GDALCopyWords(afResult, eWrkDataType, 0, abyTemp,
eVRTBandDataType, 0, 1);
CopyWord<WorkingDT>(afResult, eWrkDataType, abyTemp,
eVRTBandDataType);
GDALCopyWords(abyTemp, eVRTBandDataType, 0, pDstLocation,
eBufType, 0, 1);
}
Expand Down

0 comments on commit 664f137

Please sign in to comment.