10-band functionality in in place. Need to add/upload final 10-band s…

…ample imagery.

Alignment-RigRelatives.ipynb

@@ -465,7 +465,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.1"
+   "version": "3.7.3"
   },
   "toc": {
    "nav_menu": {},

Alignment.ipynb

@@ -179,8 +179,12 @@
     "# figsize=(30,23) # use this size for full-image-resolution display\n",
     "figsize=(16,13)   # use this size for export-sized display\n",
     "\n",
-    "rgb_band_indices = [2,1,0]\n",
-    "cir_band_indices = [3,2,1]\n",
+    "rgb_band_indices = [capture.band_names_lower().index('red'),\n",
+    "                    capture.band_names_lower().index('green'),\n",
+    "                    capture.band_names_lower().index('blue')]\n",
+    "cir_band_indices = [capture.band_names_lower().index('nir'),\n",
+    "                    capture.band_names_lower().index('red'),\n",
+    "                    capture.band_names_lower().index('green')]\n",
     "\n",
     "# Create a normalized stack for viewing\n",
     "im_display = np.zeros((im_aligned.shape[0],im_aligned.shape[1],im_aligned.shape[2]), dtype=np.float32 )\n",
@@ -360,17 +364,20 @@
     "from micasense import plotutils\n",
     "import matplotlib.pyplot as plt\n",
     "\n",
+    "nir_band = capture.band_names_lower().index('nir')\n",
+    "red_band = capture.band_names_lower().index('red')\n",
+    "\n",
     "np.seterr(divide='ignore', invalid='ignore') # ignore divide by zero errors in the index calculation\n",
     "\n",
     "# Compute Normalized Difference Vegetation Index (NDVI) from the NIR(3) and RED (2) bands\n",
-    "ndvi = (im_aligned[:,:,3] - im_aligned[:,:,2]) / (im_aligned[:,:,3] + im_aligned[:,:,2])\n",
+    "ndvi = (im_aligned[:,:,nir_band] - im_aligned[:,:,red_band]) / (im_aligned[:,:,nir_band] + im_aligned[:,:,red_band])\n",
     "\n",
     "# remove shadowed areas (mask pixels with NIR reflectance < 20%))\n",
     "if img_type == 'reflectance':\n",
-    "    ndvi = np.ma.masked_where(im_aligned[:,:,3] < 0.20, ndvi) \n",
+    "    ndvi = np.ma.masked_where(im_aligned[:,:,nir_band] < 0.20, ndvi) \n",
     "elif img_type == 'radiance':\n",
     "    lower_pct_radiance = np.percentile(im_aligned[:,:,3],  10.0)\n",
-    "    ndvi = np.ma.masked_where(im_aligned[:,:,3] < lower_pct_radiance, ndvi) \n",
+    "    ndvi = np.ma.masked_where(im_aligned[:,:,nir_band] < lower_pct_radiance, ndvi) \n",
     "    \n",
     "# Compute and display a histogram\n",
     "ndvi_hist_min = np.min(ndvi)\n",
@@ -416,7 +423,8 @@
    "outputs": [],
    "source": [
     "# Compute Normalized Difference Red Edge Index from the NIR(3) and RedEdge(4) bands\n",
-    "ndre = (im_aligned[:,:,3] - im_aligned[:,:,4]) / (im_aligned[:,:,3] + im_aligned[:,:,4])\n",
+    "rededge_band = capture.band_names_lower().index('red edge')\n",
+    "ndre = (im_aligned[:,:,nir_band] - im_aligned[:,:,rededge_band]) / (im_aligned[:,:,nir_band] + im_aligned[:,:,rededge_band])\n",
     "\n",
     "# Mask areas with shadows and low NDVI to remove soil\n",
     "masked_ndre = np.ma.masked_where(ndvi < min_display_ndvi, ndre)\n",
@@ -457,7 +465,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "if im_aligned.shape[2] >= 5:\n",
+    "if len(capture.lw_indices()) > 0:\n",
     "\n",
     "    # by default we don't mask the thermal, since it's native resolution is much lower than the MS\n",
     "    masked_thermal = im_aligned[:,:,5]\n",
@@ -506,9 +514,8 @@
    },
    "outputs": [],
    "source": [
-    "bands = ['blue','green','red','nir','rededge']\n",
-    "x_band = 2\n",
-    "y_band = 3\n",
+    "x_band = red_band\n",
+    "y_band = nir_band\n",
     "x_max = np.max(im_aligned[:,:,x_band])\n",
     "y_max = np.max(im_aligned[:,:,y_band])\n",
     "\n",
@@ -517,8 +524,8 @@
     "ax = fig.gca()\n",
     "ax.set_xlim([0,x_max])\n",
     "ax.set_ylim([0,y_max])\n",
-    "plt.xlabel(\"{} Reflectance\".format(bands[x_band]))\n",
-    "plt.ylabel(\"{} Reflectance\".format(bands[y_band]))\n",
+    "plt.xlabel(\"{} Reflectance\".format(capture.band_names()[x_band]))\n",
+    "plt.ylabel(\"{} Reflectance\".format(capture.band_names()[y_band]))\n",
     "plt.show()"
    ]
   },
@@ -566,7 +573,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.1"
+   "version": "3.7.3"
   },
   "toc": {
    "nav_menu": {},

Batch Processing.ipynb

@@ -51,7 +51,7 @@
     "outputPath = os.path.join(imagePath,'..','stacks')\n",
     "thumbnailPath = os.path.join(outputPath, '..', 'thumbnails')\n",
     "\n",
-    "overwrite = False # usefult to set to false to continue interrupted processing\n",
+    "overwrite = False # can be set to set to False to continue interrupted processing\n",
     "generateThumbnails = True\n",
     "\n",
     "# Allow this code to align both radiance and reflectance images; bu excluding\n",
@@ -305,14 +305,14 @@
    "source": [
     "import subprocess\n",
     "\n",
-    "old_dir = os.getcwd()\n",
-    "os.chdir(outputPath)\n",
-    "cmd = 'exiftool -csv=\"{}\" -overwrite_original .'.format(fullCsvPath)\n",
+    "if os.environ.get('exiftoolpath') is not None:\n",
+    "    exiftool_cmd = os.path.normpath(os.environ.get('exiftoolpath'))\n",
+    "else:\n",
+    "    exiftool_cmd = 'exiftool'\n",
+    "        \n",
+    "cmd = '{} -csv=\"{}\" -overwrite_original {}'.format(exiftool_cmd, fullCsvPath, outputPath)\n",
     "print(cmd)\n",
-    "try:\n",
-    "    subprocess.check_call(cmd)\n",
-    "finally:\n",
-    "    os.chdir(old_dir)"
+    "subprocess.check_call(cmd)"
    ]
   }
  ],
@@ -332,7 +332,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.1"
+   "version": "3.7.3"
   }
  },
  "nbformat": 4,

README.md

@@ -62,7 +62,13 @@ This code is a community effort and is not supported by MicaSense support. Pleas
 Tests for many library functions are included in the `tests` diretory. Install the `pytest` module through your package manager (e.g. `pip install pytest`) and then tests can be run from the main directory using the command:
 
 ```bash
-pytest .
+pytest
+```
+
+Test execution can be relatively slow (2-3 minutes) as there is a lot of image processing occuring in some of the tests, and quite a bit of re-used IO.  To speed up tests, install the `pytest-xdist` plugin using `conda` or `pip` and achieve a significant speed up by running tests in parallel.
+
+```bash
+pytest -n auto
 ```
 
 For faster testing we can use `pytest-xdist` to paralellize testing.

micasense/capture.py

@@ -44,7 +44,7 @@ class Capture(object):
     found in the same folder and also share the same filename prefix, such
     as IMG_0000_*.tif, but this is not required
     """
-    def __init__(self, images, panelCorners=[None]*5):
+    def __init__(self, images, panelCorners=None):
         if isinstance(images, image.Image):
             self.images = [images]
         elif isinstance(images, list):
@@ -59,7 +59,10 @@ def __init__(self, images, panelCorners=[None]*5):
         self.uuid = self.images[0].capture_id
         self.panels = None
         self.detected_panel_count = 0
-        self.panelCorners = panelCorners
+        if panelCorners is None:
+            self.panelCorners = [None]*len(self.eo_indices())
+        else:
+            self.panelCorners = panelCorners
 
         self.__aligned_capture = None
 
@@ -146,8 +149,12 @@ def center_wavelengths(self):
         return [img.center_wavelength for img in self.images]
 
     def band_names(self):
-        '''Returns a list of the image band names'''
+        '''Returns a list of the image band names as they are in the image metadata'''
         return [img.band_name for img in self.images]
+
+    def band_names_lower(self):
+        '''Returns a list of the image band names in all lower case for easier comparisons'''
+        return [img.band_name.lower() for img in self.images]
 
     def dls_present(self):
         '''Returns true if DLS metadata is present in the images'''
@@ -215,9 +222,16 @@ def compute_undistorted_reflectance(self, irradiance_list=None, force_recompute=
 
     def eo_images(self):
         return [img for img in self.images if img.band_name != 'LWIR']
+
     def lw_images(self):
         return [img for img in self.images if img.band_name == 'LWIR']
 
+    def eo_indices(self):
+        return [index for index,img in enumerate(self.images) if img.band_name != 'LWIR']
+
+    def lw_indices(self):
+        return [index for index,img in enumerate(self.images) if img.band_name == 'LWIR']
+
     def reflectance(self, irradiance_list):
         '''Comptute and return list of reflectance images for given irradiance'''
         eo_imgs = [img.reflectance(irradiance_list[i]) for i,img in enumerate(self.eo_images())]
@@ -352,32 +366,41 @@ def aligned_shape(self):
             raise RuntimeError("call Capture.create_aligned_capture prior to saving as stack")
         return self.__aligned_capture.shape
 
-    def save_capture_as_stack(self, outfilename):
+    def save_capture_as_stack(self, outfilename, sort_by_wavelength=False):
         from osgeo.gdal import GetDriverByName, GDT_UInt16
         if self.__aligned_capture is None:
             raise RuntimeError("call Capture.create_aligned_capture prior to saving as stack")
 
         rows, cols, bands = self.__aligned_capture.shape
         driver = GetDriverByName('GTiff')
+
         outRaster = driver.Create(outfilename, cols, rows, bands, GDT_UInt16, options = [ 'INTERLEAVE=BAND','COMPRESS=DEFLATE' ])
-        if outRaster is None:
-            raise IOError("could not load gdal GeoTiff driver")
-        for i in range(0,5):
-            outband = outRaster.GetRasterBand(i+1)
-            outdata = self.__aligned_capture[:,:,i]
-            outdata[outdata<0] = 0
-            outdata[outdata>2] = 2   #limit reflectance data to 200% to allow some specular reflections
-            outband.WriteArray(outdata*32768) # scale reflectance images so 100% = 32768
-            outband.FlushCache()
-
-        if bands == 6:
-            outband = outRaster.GetRasterBand(6)
-            outdata = (self.__aligned_capture[:,:,5]+273.15) * 100 # scale data from float degC to back to centi-Kelvin to fit into uint16
-            outdata[outdata<0] = 0
-            outdata[outdata>65535] = 65535
-            outband.WriteArray(outdata)
-            outband.FlushCache()
-        outRaster = None
+        try:
+            if outRaster is None:
+                raise IOError("could not load gdal GeoTiff driver")
+
+            if sort_by_wavelength:
+                eo_list = list(np.argsort(self.center_wavelengths()))
+            else:
+                eo_list = self.eo_indices()
+
+            for outband,inband in enumerate(eo_list):
+                outband = outRaster.GetRasterBand(outband+1)
+                outdata = self.__aligned_capture[:,:,inband]
+                outdata[outdata<0] = 0
+                outdata[outdata>2] = 2   #limit reflectance data to 200% to allow some specular reflections
+                outband.WriteArray(outdata*32768) # scale reflectance images so 100% = 32768
+                outband.FlushCache()
+
+            for outband,inband in enumerate(self.lw_indices()):
+                outband = outRaster.GetRasterBand(len(eo_list)+outband+1)
+                outdata = (self.__aligned_capture[:,:,inband]+273.15) * 100 # scale data from float degC to back to centi-Kelvin to fit into uint16
+                outdata[outdata<0] = 0
+                outdata[outdata>65535] = 65535
+                outband.WriteArray(outdata)
+                outband.FlushCache()
+        finally:
+            outRaster = None
 
     def save_capture_as_rgb(self, outfilename, gamma=1.4, downsample=1, white_balance='norm', hist_min_percent=0.5, hist_max_percent=99.5, sharpen=True, rgb_band_indices = [2,1,0]):
         if self.__aligned_capture is None:
@@ -414,32 +437,35 @@ def save_capture_as_rgb(self, outfilename, gamma=1.4, downsample=1, white_balanc
             imageio.imwrite(outfilename, (255*gamma_corr_rgb).astype('uint8'))
         else:
             imageio.imwrite(outfilename, (255*unsharp_rgb).astype('uint8'))
-
-    def save_thermal_over_rgb(self, outfilename, figsize=(30,23)):
+    
+    def save_thermal_over_rgb(self, outfilename, figsize=(30,23), lw_index = None, hist_min_percent = 0.2, hist_max_percent = 99.8):
         if self.__aligned_capture is None:
             raise RuntimeError("call Capture.create_aligned_capture prior to saving as RGB")
 
         # by default we don't mask the thermal, since it's native resolution is much lower than the MS
-        masked_thermal = self.__aligned_capture[:,:,5]
-
+        if lw_index is None:
+            lw_index = self.lw_indices()[0]
+        masked_thermal = self.__aligned_capture[:,:,lw_index]
+
         im_display = np.zeros((self.__aligned_capture.shape[0],self.__aligned_capture.shape[1],3), dtype=np.float32 )
-        rgb_band_indices = [2,1,0]
-        hist_min_percent = 0.2
-        hist_max_percent = 99.8
-        # for rgb true color, we usually want to use the same min and max scaling across the 3 bands to
-        # maintain the "white balance" of the calibrated image
+        rgb_band_indices = [self.band_names_lower().index('red'),
+                            self.band_names_lower().index('green'),
+                            self.band_names_lower().index('blue')]
+
+        # for rgb true color, we usually want to use the same min and max scaling across the 3 bands to 
+        # maintain the "white balance" of the calibrated image  
         im_min = np.percentile(self.__aligned_capture[:,:,rgb_band_indices].flatten(), hist_min_percent)  # modify these percentiles to adjust contrast
         im_max = np.percentile(self.__aligned_capture[:,:,rgb_band_indices].flatten(), hist_max_percent)  # for many images, 0.5 and 99.5 are good values
         for dst_band,src_band in enumerate(rgb_band_indices):
             im_display[:,:,dst_band] =  imageutils.normalize(self.__aligned_capture[:,:,src_band], im_min, im_max)
 
         # Compute a histogram
-        min_display_therm = np.percentile(masked_thermal, 1)
-        max_display_therm = np.percentile(masked_thermal, 99)
+        min_display_therm = np.percentile(masked_thermal, hist_min_percent)
+        max_display_therm = np.percentile(masked_thermal, hist_max_percent)
 
-        fig, axis = plotutils.plot_overlay_withcolorbar(im_display,
-                                            masked_thermal,
-                                            figsize=figsize,
+        fig, _ = plotutils.plot_overlay_withcolorbar(im_display,
+                                            masked_thermal, 
+                                            figsize=figsize, 
                                             title='Temperature over True Color',
                                             vmin=min_display_therm,vmax=max_display_therm,
                                             overlay_alpha=0.25,

micasense/image.py

@@ -90,9 +90,6 @@ def __init__(self, image_path, exiftool_obj=None):
         self.gain = self.meta.gain()
         self.bits_per_pixel = self.meta.bits_per_pixel()
 
-        if self.bits_per_pixel != 16:
-            NotImplemented("Unsupported pixel bit depth: {} bits".format(self.bits_per_pixel))
-
         self.vignette_center = self.meta.vignette_center()
         self.vignette_polynomial = self.meta.vignette_polynomial()
         self.distortion_parameters = self.meta.distortion_parameters()

micasense/imageutils.py

@@ -29,6 +29,7 @@
 from skimage import exposure
 from skimage.morphology import disk
 from skimage.filters import rank, gaussian
+from skimage.util import img_as_ubyte
 
 def normalize(im, min=None, max=None):
     width, height = im.shape
@@ -42,8 +43,8 @@ def normalize(im, min=None, max=None):
     return norm
 
 def local_normalize(im):
-    norm = normalize(im) # TODO: mainly using this as a type conversion, but it's expensive
-    width, height = im.shape
+    norm = img_as_ubyte(normalize(im)) # TODO: mainly using this as a type conversion, but it's expensive
+    width, _ = im.shape
     disksize = int(width/5)
     if disksize % 2 == 0:
         disksize = disksize + 1

micasense/metadata.py

@@ -100,7 +100,7 @@ def print_all(self):
     def dls_present(self):
         return self.get_item("XMP:Irradiance") is not None \
                or self.get_item("XMP:HorizontalIrradiance") is not None \
-               or self.get_item("XMP:DirectIrradiance is not None") is not None
+               or self.get_item("XMP:DirectIrradiance") is not None
 
     def supports_radiometric_calibration(self):
         if(self.get_item('XMP:RadiometricCalibration')) is None:

micasense/plotutils.py

@@ -59,7 +59,7 @@ def subplotwithcolorbar(rows, cols, images, titles=None, figsize=None):
     plt.show()
     return fig, axes
 
-def plot_overlay_withcolorbar(imgbase, imgcolor, title=None, figsize=None, vmin=None, vmax=None, overlay_alpha=1.0, overlay_colormap='viridis', overlay_steps=None, display_contours=False, contour_fmt=None, contour_steps=None, contour_alpha=None):
+def plot_overlay_withcolorbar(imgbase, imgcolor, title=None, figsize=None, vmin=None, vmax=None, overlay_alpha=1.0, overlay_colormap='viridis', overlay_steps=None, display_contours=False, contour_fmt=None, contour_steps=None, contour_alpha=None, show=True):
     ''' Plot an image with a colorbar '''
     fig, axis = plt.subplots(1, 1, figsize=figsize, squeeze=False)
     base = axis[0][0].imshow(imgbase)
@@ -81,7 +81,8 @@ def plot_overlay_withcolorbar(imgbase, imgcolor, title=None, figsize=None, vmin=
     cax = divider.append_axes("right", size="3%", pad=0.05)
     fig.colorbar(rad2, cax=cax)
     plt.tight_layout()
-    plt.show()
+    if show:
+        plt.show()
     return fig, axis[0][0]
 
 def subplot(rows, cols, images, titles=None, figsize=None):