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.

Sign up for GitHub

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

Jump to bottom

Speed-up geolocation-array based warping of Sentinel3 datasets #10812

Merged
merged 6 commits into from
Sep 27, 2024
Merged

Speed-up geolocation-array based warping of Sentinel3 datasets #10812

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
596f7d2
VRTSourcedRasterBand::IRasterIO(): initialize output buffer to nodata…
rouault Sep 16, 2024
47bc43b
VRTComplexSource::RasterIO(): use double working data type for more s…
rouault Sep 16, 2024
99a17e5
VRTComplexSource::RasterIO(): speed-up pixel copy in more cases
rouault Sep 16, 2024
b8e608f
GDALCreateGeoLocTransformer(): fix inverted logic to decide for a deb…
rouault Sep 16, 2024
a93a76c
GDALCreateGeoLocTransformer(): increase threshold to use GTiff geoloc…
rouault Sep 16, 2024
642f7cc
GDALGeoLocDatasetAccessors: use smaller, but more, cached tiles
rouault Sep 16, 2024
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
14 changes: 10 additions & 4 deletions alg/gdalgeoloc.cpp
View file
Open in desktop
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
Open in desktop
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
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2493,30 +2493,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 @@ -2553,7 +2591,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 @@ -2563,24 +2601,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 @@ -2591,18 +2621,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
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -420,7 +420,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 @@ -436,6 +438,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
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2986,10 +2986,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 @@ -3043,7 +3043,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 @@ -3077,6 +3080,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 @@ -3232,8 +3281,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 @@ -3253,8 +3302,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 @@ -3412,6 +3461,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 @@ -3548,18 +3603,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
Loading