mapbox / rio-pansharpen Goto Github PK

View Code? Open in Web Editor NEW

50.0 93.0 24.0 83.36 MB

pansharpening landsat scenes

License: MIT License

Python 100.00%

satellite imagery rasterio

rio-pansharpen's Introduction

rio-pansharpen

pansharpens Landsat 8 scenes.

What is pansharpening?

Pansharpening is a process of using the spatial information in the high-resolution grayscale band (panchromatic, or pan-band) and color information in the multispectral bands to create a single high-resolution color image.

P pan-pixel cluster + M single multispectral pixel = M pan-sharpened pixel

Find more examples and information in the Mapbox pansharpening blog post.

Install

We highly recommend installing in a virtualenv. Once activated, :

pip install -U pip
pip install rio-pansharpen

Or install from source :

git clone https://github.com/mapbox/rio-pansharpen.git
cd rio-pansharpen
pip install -U pip
pip install -r requirements.txt
pip install -e .

Python API

pansharpen.worker

The worker module pansharpens Landsat 8. Visit the USGS Landsat page page for more information on Landsat 8 band designations.

1. `worker.pansharpen`

The worker.pansharpen function accepts the following as inputs:

numpy 3D array with shape == (3, vis_height, vis_width)
affine transform defining the georeferencing of the vis array
numpy 2D array with shape == (pan_height, pan_width)
affine transform defining the georeferencing of the pan array
pansharpening method

and outputs:

numpy 3D array with shape == (3, pan_height, pan_width)

>>> from pansharpen import worker
>>> from pansharpen.methods import Brovey
...
>>> pansharpened = worker.pansharpen(vis, vis_transform, pan, pan_transform,
                       pan_dtype, r_crs, dst_crs, weight,
                       method="Brovey", src_nodata=0)

2. `worker.calculate_landsat_pansharpen`

>>> from pansharpen import worker
>>> from pansharpen.utils import _calc_windows
>>> import riomucho
...
>>> worker.calculate_landsat_pansharpen(src_paths, dst_path, dst_dtype,
        weight, verbosity, jobs, half_window,
        customwindow)

CLI

pansharpen

Usage: rio pansharpen [OPTIONS] [SRC_PATHS]... DST_PATH

  Pansharpens a landsat scene. Input is a panchromatic band (B8), plus 3 color
  bands (B4, B3, B2)

     rio pansharpen B8.tif B4.tif B3.tif B2.tif out.tif

  Or with shell expansion

     rio pansharpen LC80410332015283LGN00_B{8,4,3,2}.tif out.tif

Options:
  --dst-dtype [uint16|uint8]
  -w, --weight FLOAT          Weight of blue band [default = 0.2]
  -v, --verbosity
  -j, --jobs INTEGER          Number of processes [default = 1]
  --half-window               Use a half window assuming pan in aligned with
                              rgb bands, default: False
  -c, --customwindow INTEGER  Specify blocksize for custom windows >
                              150[default=src_blockswindows]
  --help                      Show this message and exit.
  --help                 Show this message and exit.

Comparison of Different Pansharpening Methods

We've implemented the Weighted Brovey Transform for pansharpening, which is appropriate for data like Landsat where the panchromatic band is relatively similar in resolution to the color bands.

For more information on other pansharpening methods such as IHS, PCA, P+XS, Wavelet, VWP, Wavelet with Canny Edge Detector etc, please read our notes here.

rio-pansharpen's People

Contributors

Stargazers

Watchers

rio-pansharpen's Issues

rio pansharpen command fails with AttributeError: bilinear

I am trying to get this to run, but every try I do ends up with this error. Stack trace follows:

Traceback (most recent call last): File "/home/ubuntu/anaconda2/envs/rasterio-env/bin/rio", line 11, in <module> sys.exit(main_group()) File "/home/ubuntu/anaconda2/envs/rasterio-env/lib/python2.7/site-packages/click/core.py", line 716, in __call__ return self.main(*args, **kwargs) File "/home/ubuntu/anaconda2/envs/rasterio-env/lib/python2.7/site-packages/click/core.py", line 696, in main rv = self.invoke(ctx) File "/home/ubuntu/anaconda2/envs/rasterio-env/lib/python2.7/site-packages/click/core.py", line 1060, in invoke return _process_result(sub_ctx.command.invoke(sub_ctx)) File "/home/ubuntu/anaconda2/envs/rasterio-env/lib/python2.7/site-packages/click/core.py", line 889, in invoke return ctx.invoke(self.callback, **ctx.params) File "/home/ubuntu/anaconda2/envs/rasterio-env/lib/python2.7/site-packages/click/core.py", line 534, in invoke return callback(*args, **kwargs) File "/home/ubuntu/anaconda2/envs/rasterio-env/lib/python2.7/site-packages/rio_pansharpen/scripts/cli.py", line 46, in pansharpen jobs, half_window, customwindow File "/home/ubuntu/anaconda2/envs/rasterio-env/lib/python2.7/site-packages/rio_pansharpen/worker.py", line 167, in calculate_landsat_pansharpen rm.run(jobs) File "/home/ubuntu/anaconda2/envs/rasterio-env/lib/python2.7/site-packages/riomucho/__init__.py", line 108, in run for data, window in self.pool.imap_unordered(reader_worker, self.windows): File "/home/ubuntu/anaconda2/envs/rasterio-env/lib/python2.7/site-packages/riomucho/single_process_pool.py", line 9, in imap_unordered yield reader_worker(i) File "/home/ubuntu/anaconda2/envs/rasterio-env/lib/python2.7/site-packages/riomucho/__init__.py", line 33, in manualRead raise e AttributeError: bilinear

I am running this in a conda-env. This is the output of conda env export (untruncated since the cause may be some conflict):

`Using Anaconda Cloud api site https://api.anaconda.org
name: rasterio-env
dependencies:

affine=2.0.0=py27_0
click=6.6=py27_0
click-plugins=1.0.3=py27_0
cligj=0.4.0=py27_0
curl=7.49.0=1
enum34=1.1.6=py27_0
geos=3.5.0=0
geotiff=1.4.1=0
hdf4=4.2.12=0
hdf5=1.8.17=1
jbig=2.1=0
jpeg=8d=1
kealib=1.4.6=0
libgdal=2.1.0=0
libnetcdf=4.4.1=0
libtiff=4.0.6=2
menpo::opencv3=3.1.0=py27_0
mkl=11.3.3=0
numpy=1.11.1=py27_0
openssl=1.0.2h=1
pip=8.1.2=py27_0
proj4=4.9.2=0
pyparsing=2.1.4=py27_0
python=2.7.12=1
rasterio=0.36.0=np111py27_0
readline=6.2=2
setuptools=25.1.6=py27_0
snuggs=1.4.0=py27_0
sqlite=3.13.0=0
tk=8.5.18=0
wheel=0.29.0=py27_0
xerces-c=3.1.4=0
xz=5.2.2=0
zlib=1.2.8=3
pip:
- boto==2.39.0
- cycler==0.10.0
- dask==0.11.0
- decorator==4.0.10
- futures==3.0.5
- geocoder==1.9.0
- geojson==1.3.3
- homura==0.1.3
- humanize==0.5.1
- landsat-util==0.13.0
- matplotlib==1.5.1
- networkx==1.11
- pandas==0.18.1
- pillow==3.3.1
- polyline==1.1
- pycurl==7.43.0
- pymongo==3.3.0
- pyproj==1.9.5.1
- python-dateutil==2.5.1
- pytz==2016.6.1
- ratelim==0.1.6
- requests==2.7.0
- requests-futures==0.9.7
- rio-mucho==0.2.1
- rio-pansharpen==0.2.0.dev1
- scikit-image==0.12.3
- scipy==0.17.0
- shapely==1.5.16
- six==1.8.0
- termcolor==1.1.0
- toolz==0.8.0
- usgs==0.1.9
  prefix: /home/ubuntu/anaconda2/envs/rasterio-env`

CLI as a plugin in rio namespace

given the rio naming of the repo, we should move the setup.py entry point to register this as a rasterio plugin instead of a standalone script.

Docstrings

We should define the parameters, return values, and a short description.

Update Readme

We should complete the readme with coverage status, examples for using the cli and api, instructions for installing the module etc.

pansharpen method - string or function?

We're currently mixing the type of the method argument to pansharpen.

It defaults to a string but the code expects a function

We should either pass around a string (e.g. "Brovey") instead of a function to be safe from any distributed/multiprocessing gotchas

We should pansharpen with padded windows

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:

Read a buffered window of pan data.
Read the corresponding RGB data, geo-aware (i.e., using the window’s affine).
Transform (upsample) the RGB to match the pan.
Clip off the window’s buffer.
Null out any pixel where RGB data was missing.
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.

Refactor python API

We should build a pure function (data-in, data-out, no IO) that we can unit test easily. This might take some refactoring to accomplish.

I'm thinking of a function with a signature like:

def pansharpen(vis, vis_transform, pan, pan_transform, method="Brovey"):
    """Pansharpen a lower-resolution visual band

    Parameters
    =========
    vis: ndarray, 3D with shape == (3, vh, vw)
        Visual band array with RGB bands
    vis_transform: Affine
        affine transform defining the georeferencing of the vis array
    pan: ndarray, 2D with shape == (ph, pw)
        Panchromatic band array
    pan_transform: Affine
        affine transform defining the georeferencing of the pan array
    method: string
        Algorithm for pansharpening; default Brovey

    Returns:
    ======
    pansharp: ndarray, 3D with shape == (3, ph, pw)
        pansharpened visual band
        affine transform is identical to `pan_transform`
    """

Option to output 16bit tifs

would allow us to keep color detail throughout the process for use with rio color, etc

Pan path to release [Master Ticket]

Last July, we integrated pansharpening into the l8sr pipeline and blogged about it. We'd like to clean up this current repo, modularize pansharpening and have it ready as a plugin for the new l8sr-live stack which handles all preprocessing work necessary for Landsat (https://github.com/mapbox/satellite/issues/1204#issue-167649314). Here's the checklist of things we need to complete before releasing pansharpening into the wild:

To Do:

In Progress:

Item	Dev	Test	PR ready for 👀	Merge PR
New window approach	✅ #2	✅ #1 (comment)	✅ #2	✅
Add option to output 16bit tifs	✅ #3	✅ #11	✅ #9	✅
Refactor python API	✅ #13	✅ #5	✅ #13	✅
Better test coverage	✅ #13	✅ #5	✅ #13	✅
Update Readme	✅ #14		✅	✅
Add docstrings	✅ #14		✅	✅

Completed:

New window approach | #1 (comment)
Add option to output 16bit tifs | #3
Refactor python API | #8
Better test coverage (current coverage = 96%)| #5
Updated Readme | #14
Added Docstrings | #14

pip install failed with release [0.1.0]

Problem:

Pip install failed with
IOError: [Errno 2] No such file or directory: 'README.md'

Command:

pip install rio-pansharpen

Full backtrace:

$pip install rio-pansharpen
Collecting rio-pansharpen
  Downloading rio-pansharpen-0.1.0.tar.gz
    Complete output from command python setup.py egg_info:
    Traceback (most recent call last):
      File "<string>", line 1, in <module>
      File "/private/var/folders/mp/nl3vzcd93979ktmr6rxtt48r0000gn/T/pip-build-fxrufF/rio-pansharpen/setup.py", line 6, in <module>
        with codecs_open('README.md', encoding='utf-8') as f:
      File "/Users/VirginiaMapbox/venv/lib/python2.7/codecs.py", line 896, in open
        file = __builtin__.open(filename, mode, buffering)
    IOError: [Errno 2] No such file or directory: 'README.md'

    ----------------------------------------
Command "python setup.py egg_info" failed with error code 1 in /private/var/folders/mp/nl3vzcd93979ktmr6rxtt48r0000gn/T/pip-build-fxrufF/rio-pansharpen/

Current fix implemented in `0.1.1`:

I should have used README.rst in the setup.py instead of the markdown file. For release[0.1.1], I removed setup.cfg which tells PyPI where my README file is:

[metadata]
description-file = README.md

I was able to run pip install rio-pansharpen successfully after this change. For the next release, I will include README.rst as well as try to add the README files to MANIFEST.in.

Rename to "pansharpen"

The cli command and python module are named pansharpen. Maybe the repo should be too?

Test coverage

We need to do more direct unit tests and make sure we're hitting all the code.

Add tests for uint16 bit output

Add tests for changes made in #9

make public, rio pansharpen

Is there any reason to keep this repo private? I propose we make it a proper rio pansharpen tool and spend a few hours cleaning it up for a release.

TIR pan?

Hi there,

Do you think I could use this package to pansharpening Landsat8 TIR images (Band 10), or it just work on the visible bands?
Thanks in advance,

LSR

Dealing with products other than Level1 L8 - dtype and nodata issues

Currently the src_nodata is hardcoded to zero which is problematic for pansharpening products derived from level 1 L8 Bands.

For example, using LaSRC, the surface reflectance and toa products are output as int16 (-32768 to 32768) with a nodata value of -9999. When we pansharpen these, the nodata mask is incorrect since it assumes zero.

Window TypeError

Trying to run this, I am getting the following error, probably because Window objects in rasterio no longer are tuples.

Rasterio version 1.1.2

MuchoChildError: 
Child process's traceback:
Traceback (most recent call last):
  File "/home/jovyan/.conda/envs/py3/lib/python3.6/site-packages/riomucho/__init__.py", line 63, in wrapper
    return func(*args, **kwds)
  File "/home/jovyan/.conda/envs/py3/lib/python3.6/site-packages/riomucho/__init__.py", line 119, in __call__
    return self.user_func(srcs, window, ij, global_args), window
  File "/home/jovyan/.conda/envs/py3/lib/python3.6/site-packages/rio_pansharpen/worker.py", line 77, in _pansharpen_worker
    rgb_window = _pad_window(rgb_base_window, padding)
  File "/home/jovyan/.conda/envs/py3/lib/python3.6/site-packages/rio_pansharpen/utils.py", line 116, in _pad_window
    (wnd[0][0] - pad, wnd[0][1] + pad),
TypeError: 'Window' object does not support indexing

ValueError: blockxsize exceeds raster width with

Not sure if this is a problem with rio-pansharpen itself or rio-mucho.

I tried running the pansharpening over cropped Landsat images, but a particular window size keeps giving me this error, no matter if I try running from Python or directly via rio pansharpen.

After some tests, I managed to create the following repeatable sequence of actions to reproduce the problem:

Download some arbitrary Landsat scene. In my tests, any scene will cause the same problem.

landsat download -b 8,4,3,2 LC80390282016240LGN00

Crop the scene using gdal_translate to this exact window:

gdal_translate -srcwin 3812 2460 366 740 /home/ubuntu/landsat/downloads/LC82220762016210LGN00/LC82220762016210LGN00_B8.TIF gdal_b8.tif

gdal_translate -srcwin 1906 1230 183 371 /home/ubuntu/landsat/downloads/LC82220762016210LGN00/LC82220762016210LGN00_B4.TIF gdal_b4.tif

gdal_translate -srcwin 1906 1230 183 371 /home/ubuntu/landsat/downloads/LC82220762016210LGN00/LC82220762016210LGN00_B3.TIF gdal_b3.tif

gdal_translate -srcwin 1906 1230 183 371 /home/ubuntu/landsat/downloads/LC82220762016210LGN00/LC82220762016210LGN00_B2.TIF gdal_b2.tif

Run the pansharpening:

rio pansharpen gdal_b{8,4,3,2}.tif pan.tif

The full output is as follows:

/home/ubuntu/anaconda2/envs/benchmark/lib/python2.7/site-packages/rio_pansharpen/scripts/cli.py:46: FutureWarning: The value of this property will change in version 1.0. Please see https://github.com/mapbox/rasterio/issues/86 for details.
  jobs, half_window, customwindow
/home/ubuntu/anaconda2/envs/benchmark/lib/python2.7/site-packages/rio_pansharpen/worker.py:134: FutureWarning: GDAL-style transforms are deprecated and will not be supported in Rasterio 1.0.
  transform=guard_transform(pan_src.transform),
/home/ubuntu/anaconda2/envs/benchmark/lib/python2.7/site-packages/rio_pansharpen/worker.py:160: FutureWarning: GDAL-style transforms are deprecated and will not be supported in Rasterio 1.0.
  "r_aff": guard_transform(r_meta['transform']),
Traceback (most recent call last):
  File "/home/ubuntu/anaconda2/envs/benchmark/lib/python2.7/site-packages/rio_pansharpen/worker.py", line 167, in calculate_landsat_pansharpen
    rm.run(jobs)
  File "/home/ubuntu/anaconda2/envs/benchmark/lib/python2.7/site-packages/riomucho/__init__.py", line 108, in run
    with rio.open(self.outpath, 'w', **self.options) as dst:
  File "/home/ubuntu/anaconda2/envs/benchmark/lib/python2.7/site-packages/rasterio/__init__.py", line 193, in open
    s.start()
  File "rasterio/_io.pyx", line 1401, in rasterio._io.RasterUpdater.start (rasterio/_io.c:19697)
ValueError: blockxsize exceeds raster width.
Traceback (most recent call last):
  File "/home/ubuntu/anaconda2/envs/benchmark/bin/rio", line 11, in <module>
    sys.exit(main_group())
  File "/home/ubuntu/anaconda2/envs/benchmark/lib/python2.7/site-packages/click/core.py", line 716, in __call__
    return self.main(*args, **kwargs)
  File "/home/ubuntu/anaconda2/envs/benchmark/lib/python2.7/site-packages/click/core.py", line 696, in main
    rv = self.invoke(ctx)
  File "/home/ubuntu/anaconda2/envs/benchmark/lib/python2.7/site-packages/click/core.py", line 1060, in invoke
    return _process_result(sub_ctx.command.invoke(sub_ctx))
  File "/home/ubuntu/anaconda2/envs/benchmark/lib/python2.7/site-packages/click/core.py", line 889, in invoke
    return ctx.invoke(self.callback, **ctx.params)
  File "/home/ubuntu/anaconda2/envs/benchmark/lib/python2.7/site-packages/click/core.py", line 534, in invoke
    return callback(*args, **kwargs)
  File "/home/ubuntu/anaconda2/envs/benchmark/lib/python2.7/site-packages/rio_pansharpen/scripts/cli.py", line 46, in pansharpen
    jobs, half_window, customwindow
  File "/home/ubuntu/anaconda2/envs/benchmark/lib/python2.7/site-packages/rio_pansharpen/worker.py", line 167, in calculate_landsat_pansharpen
    rm.run(jobs)
  File "/home/ubuntu/anaconda2/envs/benchmark/lib/python2.7/site-packages/riomucho/__init__.py", line 108, in run
    with rio.open(self.outpath, 'w', **self.options) as dst:
  File "/home/ubuntu/anaconda2/envs/benchmark/lib/python2.7/site-packages/rasterio/__init__.py", line 193, in open
    s.start()
  File "rasterio/_io.pyx", line 1401, in rasterio._io.RasterUpdater.start (rasterio/_io.c:19697)
ValueError: blockxsize exceeds raster width.

Alpha band optional

For some workflows, we should drop the alpha band. An --out-alpha CLI option and a corresponding out_alpha function argument should handle it.