Merge branch 'master' of github.com:mapbox/rio-pansharpen

README.rst

@@ -17,7 +17,7 @@ Pansharpening is a process of using the spatial information in the high-resoluti
  
     P pan-pixel cluster + M single multispectral pixel = M pan-sharpened pixel
 
-Find more examples and information in the [Mapbox pansharpening blog post](https://www.mapbox.com/blog/l8-pansharpening/).
+Find more examples and information in the `Mapbox pansharpening blog post <https://www.mapbox.com/blog/l8-pansharpening/>`_.
 
 Install
 =======
@@ -38,13 +38,13 @@ Or install from source
     pip install -e .
 
 
-**********
+
 Python API
-**********
+==========
 
 pansharpen.worker
 -----------------
-The `worker` module pansharpens Landsat 8. Visit the `USGS Landsat page <http://landsat.usgs.gov/band_designations_landsat_satellites.php>`__ page for more information on Landsat 8 band designations.
+The ``worker`` module pansharpens Landsat 8. Visit the `USGS Landsat page <http://landsat.usgs.gov/band_designations_landsat_satellites.php>`_ page for more information on Landsat 8 band designations.
 
 1. ``worker.pansharpen``
 ------------------------
@@ -71,8 +71,8 @@ and outputs:
 
 
 
-2.`worker.calculate_landsat_pansharpen``
-----------------------------------------
+2. ``worker.calculate_landsat_pansharpen``
+------------------------------------------
 ::
 
     >>> from pansharpen import worker
@@ -84,9 +84,9 @@ and outputs:
             customwindow)
 
 
-***
+
 CLI
-***
+===
 
 
 pansharpen
@@ -121,9 +121,8 @@ pansharpen
 
 
 
-*********************************************
 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](https://github.com/mapbox/pansharpening/blob/master/docs/pansharpening_methods.md).
+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 <https://github.com/mapbox/pansharpening/blob/master/docs/pansharpening_methods.rst>`_.

docs/pansharpening_methods.md → docs/pansharpening_methods.rst

@@ -1,48 +1,60 @@
-## Comparison of Different Pansharpening Methods
+=============================================
+Comparison of Different Pansharpening Methods
+=============================================
 
-#### Brovey
+Brovey
+------
 
 The Brovey transformation is a sharpening method that uses a mathematical combination of the color image and high resolution data. Each resampled, multispectral pixel is multiplied by the ratio of the corresponding panchromatic pixel intensity to the sum of all the multispectral intensities. It assumes that the spectral range spanned by the panchromatic image is the same as that covered by the multispectral channels. This is done essentially by increasing the resolution of the color information in the data set to match that of the panchromatic band. Therefore, the output RGB images will have the pixel size of the input high-resolution panchromatic data. Various resampling methods include bilinear, lanczos, cubic, average, mode, min, and max.
 
-#### Weighted Brovey
+Weighted Brovey
+---------------
 
 Particularly for Landsat 8 imagery data, we know that the pan-band does not include the full blue band, so we take a fraction of blue (optimal weight computed in this sprint) in the pan-band and use this weight to compute the sudo_pan_band, which is a weighted average of the three bands. We then compute the ratio between the pan-band and the sudo-band and adjust each of the three bands by this ratio.
 
-```
-sudo_pan = (R + B + B * weight)/(2 + weight)
-ratio = pan / sudo_pan
-R_out = R * ratio
-G_out = G * ratio
-B_out = B * ratio
-```
-![screen shot 2015-04-13 at 10 14 29 pm](https://cloud.githubusercontent.com/assets/4450007/7141761/7a277a88-e288-11e4-9dd7-39e3f970603f.png)
+::
 
-#### IHS
+	sudo_pan = (R + B + B * weight)/(2 + weight)
+	ratio = pan / sudo_pan
+	R_out = R * ratio
+	G_out = G * ratio
+	B_out = B * ratio
+
+
+.. image:: https://cloud.githubusercontent.com/assets/4450007/7141761/7a277a88-e288-11e4-9dd7-39e3f970603f.png
+
+IHS
+---
 The IHS transformation first converts the color image to IHS color space. It then replaces the intensity band by the a weighted version of the panchromatic image. Finally, the fused image is converted back to RGB space. IHS fused images generally experience spectral distortion from the original multispectral image. 
 
-#### PCA:
+PCA:
+----
 The Principle Component Analysis (PCA) is a common statistical procedure that is used to reduce the dimensionality of multi-dimensional space. It is used for numerous applications in fields like statistics, machine learning, and signal processing.It is an orthogonal transformation that converts a set of correlated observations into a set of linearly uncorrelated values called principal components. This transformation leads to an interesting result - the first principal component accounts for the greatest proportion of variability in the data. 
 In this case, it can be used to convert intercorrelated multispectral bands into a set of uncorrelated components.  The first band, which has the highest variance, is then replaced by the Panchromatic image.  We can then obtain the high-resolution pansharpened image by applying an inverse PCA on the PCA.
 
 
-#### P+XS:
+P+XS:
+-----
 The P+XS is a variational method, which calculates the pansharpened image by minimizing an energy functional. It obtains the edge information of the panchromatic image by using the gradient. The spectral information is obtained by approximating the panchromatic image as a linear combination of the multispectral bands (Ballester, 2007). 
 
-#### Wavelet:
+Wavelet:
+--------
 The Wavelet method uses Discrete Wavelet Transforms (DWT) to decompose the original multispectral and panchromatic image into components. For each of the images, there is one component that contains low-resolution information, while the others contain more detailed local spatial information. The low-resolution component of the panchromatic image is replaced by the low-resolution multispectral component. The final image is created by performing an inverse wavelet transformation. 
 Runtime in Theory: O(n)
 
-![dwt1](https://cloud.githubusercontent.com/assets/4450007/7141344/da8b2cec-e285-11e4-9253-a3040b076bd2.jpg)
+.. image:: https://cloud.githubusercontent.com/assets/4450007/7141344/da8b2cec-e285-11e4-9253-a3040b076bd2.jpg
 
 where cJk are the scaling coefficients and djk are the wavelet coefficients. The first term in Eq. (8) gives the low-resolution approximation of the signal while the second term gives the detailed information at resolutions from the original down to the current resolution J. The process of applying the DWT can be represented as a bandk of filters, as in the figure below.
 
-![dwt_components](https://cloud.githubusercontent.com/assets/4450007/7141354/e990abcc-e285-11e4-984b-b35c22e18f6f.jpg)
+.. image:: https://cloud.githubusercontent.com/assets/4450007/7141354/e990abcc-e285-11e4-984b-b35c22e18f6f.jpg
 
 
 In case of a 2D image, a single level decomposition can be performed resulting in four different frequency bands namely LL, LH, HL and HH sub band and an N level decomposition can be performed resulting in 3N+1 different frequency bands and it is shown in figure 3. At each level of decomposition, the image is split into high frequency and low frequency components; the low frequency components can be further decomposed until the desired resolution is reached.
 
-#### VWP: 
+VWP: 
+----
 VWP combines the Wavelet method with P+XS. It uses the geometry matching term from P+XS and spectral information from wavelet decomposition. This method outperforms others by preserving the highest spectral quality (Moeller, 2008).
 
-#### Wavelet + Canny Edge Detector:
-Combining DWT with Canny Edge Detectors. More details [here](http://link.springer.com/chapter/10.1007%2F978-3-642-21783-8_6).
+Wavelet + Canny Edge Detector:
+------------------------------
+Combining DWT with Canny Edge Detectors. More details `here <http://link.springer.com/chapter/10.1007%2F978-3-642-21783-8_6>`_.

rio_pansharpen/scripts/cli.py

@@ -2,6 +2,7 @@
 
 import click
 from rio_pansharpen.worker import calculate_landsat_pansharpen
+from rasterio.rio.options import creation_options
 
 
 @click.command('pansharpen')
@@ -24,9 +25,13 @@
               default=0,
               help="Specify blocksize for custom windows > 150"
               "[default=src_blockswindows]")
-def pansharpen(src_paths, dst_path, dst_dtype,
+@click.option('--out-alpha/--no-out-alpha', default=True, is_flag=True,
+              help="Output an alpha band along with RGB")
+@creation_options
+def pansharpen(
+        src_paths, dst_path, dst_dtype,
         weight, verbosity, jobs,
-        half_window, customwindow):
+        half_window, customwindow, out_alpha, creation_options):
     """Pansharpens a landsat scene.
     Input is a panchromatic band, plus 3 color bands
 
@@ -43,5 +48,5 @@ def pansharpen(src_paths, dst_path, dst_dtype,
 
     return calculate_landsat_pansharpen(
         src_paths, dst_path, dst_dtype, weight, verbosity,
-        jobs, half_window, customwindow
-      )
+        jobs, half_window, customwindow, out_alpha, creation_options
+      )

rio_pansharpen/utils.py

@@ -135,7 +135,7 @@ def _calc_windows(pan_src, customwindow):
     return windows
 
 
-def _rescale(arr, ndv, dst_dtype):
+def _rescale(arr, ndv, dst_dtype, out_alpha=True):
     """Convert an array from output dtype, scaling up linearly
     """
     if dst_dtype == np.__dict__['uint16']:
@@ -144,12 +144,13 @@ def _rescale(arr, ndv, dst_dtype):
         # convert to 8bit value range in place
         scale = float(np.iinfo(np.uint16).max) / float(np.iinfo(np.uint8).max)
 
-    return np.concatenate(
-                [
-                    (arr / scale).astype(dst_dtype),
-                    _simple_mask(
-                        arr.astype(dst_dtype),
-                        (ndv, ndv, ndv)
-                    ).reshape(1, arr.shape[1], arr.shape[2])
-                ]
-            )
+    res = (arr / scale).astype(dst_dtype)
+
+    if out_alpha:
+        mask = _simple_mask(
+            arr.astype(dst_dtype),
+            (ndv, ndv, ndv)).reshape(
+                1, arr.shape[1], arr.shape[2])
+        return np.concatenate([res, mask])
+    else:
+        return res