add -ms_shift option to gdal_pansharpen.py

[tbonfort]

No description provided.

[tbonfort]

@rouault is this the correct wording to describe msshift ?

[jratike80]

I am not sure if the name is logical. Shift makes me think that -ms_shift 0 0 wouldn't do anything, but is is so that actually ms_shift 1 1 does not do anything?

[tbonfort]

@jratike80 I'm not sure to understand your comment. For the name of the flag, ms_shift is just replicating the MSShift tags that are already present in the VRT. For the actual values, 0 0 does nothing, 1 1 will offset the ms bands by 1 pixel in the bottom-right direction.

[jratike80]

I read the comment e.g. -ms_shift 1 2 specifies that the top left multispectral pixel corresponds to the top left of the panchromatic pixel of row 2 and column 1. So I thought that values 0 0 would correspond to row 0 and col 0, thus worth one pixel outside the panchromatic image. Maybe worth trying to clarify that the unit of the shift is in panchromatic pixels.

[rouault]

@tbonfort How you describe the option would be the logical thing, but looking at the code in GDALPansharpenOperation::ProcessRegion(), my impression is that unfortunately it has currently the reverse convention. And I can't find anything in the doc & tests about those options.
I believe:

• the code should be adjusted to reflect the semantics you propose (hopefully that shouldn't break a lot of people...)
• autotest/gdrivers/vrtpansharpen.py should test the MSShiftX and MSShiftY elements
• doc/source/drivers/raster/vrt.rst should mention them in the section related to pansharpening towards the end
• autotest/pyscripts/test_gdal_pansharpen.py should test the new switch

[tbonfort]

Hijacking my own PR to understand a bit more how the pansharpening windows are computed...

Looking at the code, it seems to me that the relative geotransforms of the pan and ms bands are not used at all in order to compute the target MS window for a given PAN window. i.e. (supposing ms_shift is 0 for now), the code assumes that the PAN has a 1,0,0,0,-1,0 geotransform, and the MS has a pan_cols/ms_cols,0,0,0,-pan_rows/ms_rows,0 geotransform. Is that correct?

Taking ms_shift into account will shift the offset of the ms window, and as Even pointed out is effectively reversed as to what would be the logical thinking, i.e. an ms_shift of 1,1 currently implies that that the ms pixels are shifted 1 pixel in the top-left direction compared to the pan band.

If my previous assumptions are correct, there is I believe an issue with the current implementation. Taking for example SPOT imagery (PNEO is similar, but not Pleiades), the relative positioning of the pan and ms planes can be found on page 57 of https://earth.esa.int/eogateway/documents/20142/37627/SPOT-6-7-imagery-user-guide.pdf . In order to obtain correct pansharpening with this layout, we should set ms_shift to -0.5,-0.5, but around lines 1107 of gdalpansharpen.cpp we should also replace

double dfRatioX = static_cast<double>(poPanchroBand->GetXSize()) /
                      aMSBands[0]->GetXSize();
double dfRatioY = static_cast<double>(poPanchroBand->GetYSize()) /
                      aMSBands[0]->GetYSize();

by

double dfRatioX = static_cast<double>(poPanchroBand->GetXSize()-1) /
                      aMSBands[0]->GetXSize();
double dfRatioY = static_cast<double>(poPanchroBand->GetYSize()-1) /
                      aMSBands[0]->GetYSize();

or in other words, the pan/ms pixel size ratio is defined by the sensor itself, and not by the ratio of the image dimensions.

I see two possible fixes for this:

1. compute the current ms_shifts and x/y ratios from the input image geotransforms if they are present, and fallback/keep the current ones if no geotransforms are present
2. and/or explicitely add a msgeotransform tag to the pansharpening vrt that specifies the MS band's relative position/scaling compared to the PAN one.

[rouault]

Looking at the code, it seems to me that the relative geotransforms of the pan and ms bands are not used at all in order to compute the target MS window for a given PAN window

no, they are not. The current state of the code related to shifts between the pan and ms bands should be mostly considered as experimental, because at the time where I developed this, I was highly confused by what I saw in test products.

e.g with the test product I had during the dev:

$ gdalinfo IMG_SPOT6_MS_201407081617172_SEN_1278748101_R1C1.JP2
Size is 3072, 3584
Origin = (1.000000000000000,-1.000000000000000)
Pixel Size = (4.000000000000000,-4.000000000000000)

$ gdalinfo IMG_SPOT6_P_201407081617172_SEN_1278748101_R1C1.JP2 
Size is 12288, 14336
Origin = (0.500000000000000,-0.500000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)

The geotransforms here come from .J2W files. I now see that seems consistent with the figure at page 57 of the SPOT manual.

we should also replace by

double dfRatioX = static_cast<double>(poPanchroBand->GetXSize()-1) /
                      aMSBands[0]->GetXSize();
double dfRatioY = static_cast<double>(poPanchroBand->GetYSize()-1) /
                      aMSBands[0]->GetYSize();

I can't make sense of that. Can you explain why we should subtract 1 to the size of the pan ? SPOT products typically seem to have exactly 4 times more pixels in width/height for PAN compared to MS, and I'd expect the ratio to be exactly 4, not 3.9999something as the above modified formula would return.

I see two possible fixes for this:

• compute the current ms_shifts and x/y ratios from the input image geotransforms if they are present, and fallback/keep the current ones if no geotransforms are present

The pixel size of the geotransform, when present, could probably be used to compute the ratio, if using the raster dimension isn't reliable. The shifts are more tricky as shown in the SPOT case.

• and/or explicitely add a msgeotransform tag to the pansharpening vrt that specifies the MS band's relative position/scaling compared to the PAN one.

could be an idea. There would be an interesting decision to decide if this geotransform is to transform the top-left corner of the top-left pixel of the MS to the top-left corner of top-left pixel of the PAN, or the center of the top-left pixel of the MS to the center of the top-left pixel of the PAN. To be consistent with GDAL's general conventions, the first hypothesis would probably be better.

[tbonfort]

$ gdalinfo IMG_SPOT6_MS_201407081617172_SEN_1278748101_R1C1.JP2
Size is 3072, 3584
Origin = (1.000000000000000,-1.000000000000000)
Pixel Size = (4.000000000000000,-4.000000000000000)

$ gdalinfo IMG_SPOT6_P_201407081617172_SEN_1278748101_R1C1.JP2 
Size is 12288, 14336
Origin = (0.500000000000000,-0.500000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)

The geotransforms here come from .J2W files. I now see that seems consistent with the figure at page 57 of the SPOT manual.

Yes, the geotransforms provided should consistently position the MS band w.r.t. the PAN band.

I'm not sure where your test product came from, and/or how it was produced. I thought the RxCx JP2 tiles were supposed to be of a fixed size of 8192x8192 (except for the right most column and bottom most row).

we should also replace by

double dfRatioX = static_cast<double>(poPanchroBand->GetXSize()-1) /
                      aMSBands[0]->GetXSize();
double dfRatioY = static_cast<double>(poPanchroBand->GetYSize()-1) /
                      aMSBands[0]->GetYSize();

I can't make sense of that. Can you explain why we should subtract 1 to the size of the pan ? SPOT products typically seem to have exactly 4 times more pixels in width/height for PAN compared to MS, and I'd expect the ratio to be exactly 4, not 3.9999something as the above modified formula would return.

I understand your perplexity given your initial test case. The -1 comes from the schematics on page 57, and/or can be recomputed by using the size of the whole dimap rather than the size of an individual jp2 tile. For example my current test image has sizes:

Pan: Size is 14809, 16797
MS: Size is 3702, 4199

which is inline with my -1 example, but in essence just means that ratios to use should be computed from the geotransforms rather than from the image dimensions.

The pixel size of the geotransform, when present, could probably be used to compute the ratio, if using the raster dimension isn't reliable. The shifts are more tricky as shown in the SPOT case.

I would be in favor of correcting this behavior now, and leaving it up to the advanced user to select the correct ms_shifts.

could be an idea. There would be an interesting decision to decide if this geotransform is to transform the top-left corner of the top-left pixel of the MS to the top-left corner of top-left pixel of the PAN, or the center of the top-left pixel of the MS to the center of the top-left pixel of the PAN. To be consistent with GDAL's general conventions, the first hypothesis would probably be better.

This would be a welcome enhancement as I believe it is less error prone and more explicit than the current ms_shift hack. Basically, ms_window = pan_window * pan_gt * ms_gt^-1 . I would also be very much in favor of sticking with gdal's corner convention in order to avoid another source of confusion, it would be up to the user to compute the adjusted geotransform if the specification is pixel-centric. (FYI, the pleiades layout is pixel centric, i.e. the center of the top-left pan pixel is aligned with the center of the top-left ms pixel)

[rouault]

superseded per PR #7373