We should pansharpen with padded windows

[celoyd]

Right now, our manually doubled affine numbers give us a slight multispectral/panchromatic offset. This isn’t a big deal today, but it’s something we should know how to handle.

We’ve used other methods in the past, but they’re prone to seams in unpredictable cases depending on pixel alignment. Per voice w/ @dnomadb, @virginiayung, @perrygeo, the most bulletproof way would look something like:

1. Read a buffered window of pan data.
2. Read the corresponding RGB data, geo-aware (i.e., using the window’s affine).
3. Transform (upsample) the RGB to match the pan.
4. Clip off the window’s buffer.
5. Null out any pixel where RGB data was missing.
6. Write out the window.

As a bonus, buffered windows would also solve the color-seam problem that would appear if we applied this code to commercial imagery.

[virginiayung]

I tested out @perrygeo's changes in #2 with LC81170442016179LGN00. This is a Landsat 8 scene of Kaohsiung, Taiwan. In the Northern Hemisphere in summer (June, July, August), the sun rises in the northeast, peaks out nearly straight overhead (depending on latitude), and then sets in the northwest. We should not be seeing shadows beneath the ship as the first example below. This could be an example of the multispectral/panchromatic offset that resulted from manually doubling the affine numbers.

👎 Previous half-window approach: Shadow beneath ship

👍 New approach utilizing rasterio's internal logic for determining windows, bounds and affine transforms: Shadow evenly around ship

Here's a gif just for fun:

Next steps:

I'll look for a couple more examples to confirm this.
cc @celoyd @perrygeo

[virginiayung]

Here's a quick note on running pansharpening on LC81170442016179LGN00 with the --half-window option:

I noticed seams at the edge of each processing window, similar to what we observed back in July 2015:

Zooming in:

This can be fixed by make_affine function, which uses an arbitrary transformation instead of the whole-dataset affine transformation for individual windows.

@perrygeo: @celoyd and I would love to learn more about how rasterio does affine transformation before we decide on ditching the --half-windowapproach for good :)

[virginiayung]

To recap my conversation with @celoyd and @perrygeo over slack yesterday, this is the best example I can find that shows how well the new window approach implemented in #2 is able to remove the artifacts created by our current half-window approach. As seen below, the example checks out visually. I haven't found examples where the new approach introduced edge effects along the edges of each block. Let's merge this into master for now and write more unit tests that specially checks for pixel alignment after we refactor the python API:

Old Half-Window Approach (assuming the panchromatic band is exactly half the resolution but otherwise aligned with RGB bands)

Old Approach (red arrows) versus New Approach:

More notes from our slack conversation explaining the two approaches

• The upper left corner of the upper left cell of pan band is 7.5 meters to the NW

$ rio info LC80380302015166LGN00_B8.TIF | jq .transform
[15,0,382492.5,0,-15,4899607.5,0,0,1]

$ rio info LC80380302015166LGN00_B4.TIF | jq .transform
[30, 0, 382485,0,-30, 4899615,0,0,1]

When we use the half-window approach, cells are no longer aligned so we would need to resample or else we would get a 7.5m shift to the NW. By aligning the affine transform when doing upsampling, we are artificially shifting the pan bands NW when really they need to be resampled. The upper left red value gets assigned to the 4 upper left pan cells:

Artificially shifting the pan bands NW may be the main reason why we were seeing a slight multispectral/panchromatic offset. Instead, the new approach uses rasterio methods to determine the windows for the pan and RGB band.

General reference on Affine transforms: http://www.perrygeo.com/python-affine-transforms.html

cc @dnomadb