mapbox / rio-pansharpen Goto Github PK

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

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`

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

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.

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

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

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`
    """

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.

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

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

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.

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.

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

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

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

Add tests for changes made in #9

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:

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

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

2. 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

3. 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.

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:

• Better test coverage (master ticket for improving all things relating to testing): #5
• Refactor Python API | #8
• CircleCI, Codecov integration
• Pep8
• Include Docstrings: #6
• Rename repo: #7
• Update readme:#12
• Make repo public

In Progress:

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

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.