diff --git a/.github/workflows/gh-pages.yml b/.github/workflows/gh-pages.yml index 8c4130f..90e0185 100644 --- a/.github/workflows/gh-pages.yml +++ b/.github/workflows/gh-pages.yml @@ -19,7 +19,7 @@ jobs: build: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v5 + - uses: actions/checkout@v7 with: lfs: true - uses: actions/setup-python@v6 @@ -36,8 +36,8 @@ jobs: - name: Install dependencies run: pip install -e .[docs] - name: Build GitHub Pages site - run: scripts/build_gh_pages.sh - - uses: actions/upload-pages-artifact@v3 + run: EXECUTE=1 scripts/build_gh_pages.sh + - uses: actions/upload-pages-artifact@v5 with: path: _site @@ -49,4 +49,4 @@ jobs: url: ${{ steps.deployment.outputs.page_url }} steps: - id: deployment - uses: actions/deploy-pages@v4 + uses: actions/deploy-pages@v5 diff --git a/.github/workflows/pr-agent.yml b/.github/workflows/pr-agent.yml index ab5f2cc..d357fb4 100644 --- a/.github/workflows/pr-agent.yml +++ b/.github/workflows/pr-agent.yml @@ -13,7 +13,7 @@ jobs: contents: write steps: - name: PR Agent action step - uses: the-pr-agent/pr-agent@v0.36.1 + uses: the-pr-agent/pr-agent@v0.38.0 env: OPENAI_KEY: ${{ secrets.OPENAI_KEY }} GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} diff --git a/.github/workflows/publish.yml b/.github/workflows/publish.yml index cf59f60..626a61d 100644 --- a/.github/workflows/publish.yml +++ b/.github/workflows/publish.yml @@ -8,7 +8,7 @@ jobs: lint: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v5 + - uses: actions/checkout@v7 - uses: actions/setup-python@v6 with: python-version: "3.12" @@ -18,7 +18,7 @@ jobs: needs: [lint] runs-on: ubuntu-latest steps: - - uses: actions/checkout@v5 + - uses: actions/checkout@v7 with: lfs: true - uses: actions/setup-python@v6 @@ -41,7 +41,7 @@ jobs: python -m pip install build python -m build - name: Store the distribution packages - uses: actions/upload-artifact@v6 + uses: actions/upload-artifact@v7 with: name: python-package-distributions path: dist/ @@ -57,7 +57,7 @@ jobs: id-token: write # IMPORTANT: mandatory for trusted publishing steps: - name: Download all the dists - uses: actions/download-artifact@v7 + uses: actions/download-artifact@v8 with: name: python-package-distributions path: dist/ diff --git a/.github/workflows/python-test.yml b/.github/workflows/python-test.yml index 79a7c55..8399423 100644 --- a/.github/workflows/python-test.yml +++ b/.github/workflows/python-test.yml @@ -11,7 +11,7 @@ jobs: lint: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v5 + - uses: actions/checkout@v7 - uses: actions/setup-python@v6 with: python-version: "3.12" @@ -25,7 +25,7 @@ jobs: matrix: python-version: ["3.10", "3.11", "3.12", "3.13"] steps: - - uses: actions/checkout@v5 + - uses: actions/checkout@v7 with: lfs: true - uses: actions/setup-python@v6 @@ -41,7 +41,7 @@ jobs: run: pip install -e ".[test]" - name: Test with pytest run: pytest - - uses: actions/upload-artifact@v6 + - uses: actions/upload-artifact@v7 if: ${{ matrix.python-version == '3.12' }} with: name: coverage-report diff --git a/.gitignore b/.gitignore index 7a6a429..f120624 100644 --- a/.gitignore +++ b/.gitignore @@ -8,3 +8,5 @@ __pycache__ /dist/ /*.egg-info /_site/ +.idea +docs/*.html diff --git a/Alignment v2.ipynb b/Alignment v2.ipynb index 93d18d9..7bfeb7e 100644 --- a/Alignment v2.ipynb +++ b/Alignment v2.ipynb @@ -697,9 +697,9 @@ "# Compute and display a histogram\n", "# ndvi_hist_min = np.min(ndvi)\n", "# ndvi_hist_max = np.max(ndvi)\n", - "ndvi_hist_min = np.min(np.percentile(ndvi, 0.5))\n", - "ndvi_hist_max = np.max(np.percentile(ndvi, 99.5))\n", - "fig, axis = plt.subplots(1, 1, figsize=(10, 4))\n", + "ndvi_hist_min = np.min(np.percentile(ndvi.compressed(), 0.5))\n", + "ndvi_hist_max = np.max(np.percentile(ndvi.compressed(), 99.5))\n", + "fig, axis = plt.subplots(1, 1, figsize=(10,4))\n", "axis.hist(ndvi.ravel(), bins=512, range=(ndvi_hist_min, ndvi_hist_max))\n", "plt.title(\"NDVI Histogram\")\n", "plt.show()\n", @@ -707,7 +707,7 @@ "min_display_ndvi = 0.45 # further mask soil by removing low-ndvi values\n", "# min_display_ndvi = np.percentile(ndvi.flatten(), 5.0) # modify with these percentilse to adjust contrast\n", "max_display_ndvi = np.percentile(\n", - " ndvi.flatten(), 99.5\n", + " ndvi.compressed(), 99.5\n", ") # for many images, 0.5 and 99.5 are good values\n", "masked_ndvi = np.ma.masked_where(ndvi < min_display_ndvi, ndvi)\n", "\n", @@ -752,15 +752,15 @@ "masked_ndre = np.ma.masked_where(ndvi < min_display_ndvi, ndre)\n", "\n", "# Compute a histogram\n", - "ndre_hist_min = np.min(np.percentile(masked_ndre, 0.5))\n", - "ndre_hist_max = np.max(np.percentile(masked_ndre, 99.5))\n", - "fig, axis = plt.subplots(1, 1, figsize=(10, 4))\n", + "ndre_hist_min = np.min(np.percentile(masked_ndre.compressed(), 0.5))\n", + "ndre_hist_max = np.max(np.percentile(masked_ndre.compressed(), 99.5))\n", + "fig, axis = plt.subplots(1, 1, figsize=(10,4))\n", "axis.hist(masked_ndre.ravel(), bins=512, range=(ndre_hist_min, ndre_hist_max))\n", "plt.title(\"NDRE Histogram (filtered to only plants)\")\n", "plt.show()\n", "\n", - "min_display_ndre = np.percentile(masked_ndre, 5)\n", - "max_display_ndre = np.percentile(masked_ndre, 99.5)\n", + "min_display_ndre = np.percentile(masked_ndre.compressed(), 5)\n", + "max_display_ndre = np.percentile(masked_ndre.compressed(), 99.5)\n", "\n", "fig, axis = plotutils.plot_overlay_withcolorbar(\n", " gamma_corr_rgb,\n", diff --git a/Alignment-10Band.ipynb b/Alignment-10Band.ipynb index ef22176..2fa43d6 100644 --- a/Alignment-10Band.ipynb +++ b/Alignment-10Band.ipynb @@ -403,9 +403,9 @@ "plt.show()\n", "\n", "min_display_ndvi = 0.45 # further mask soil by removing low-ndvi values\n", - "# min_display_ndvi = np.percentile(ndvi.flatten(), 5.0) # modify with these percentilse to adjust contrast\n", + "# min_display_ndvi = np.percentile(ndvi.compressed(), 5.0) # modify with these percentilse to adjust contrast\n", "max_display_ndvi = np.percentile(\n", - " ndvi.flatten(), 99.5\n", + " ndvi.compressed(), 99.5\n", ") # for many images, 0.5 and 99.5 are good values\n", "masked_ndvi = np.ma.masked_where(ndvi < min_display_ndvi, ndvi)\n", "\n", @@ -456,8 +456,8 @@ "plt.title(\"NDRE Histogram (filtered to only plants)\")\n", "plt.show()\n", "\n", - "min_display_ndre = np.percentile(masked_ndre, 5)\n", - "max_display_ndre = np.percentile(masked_ndre, 99.5)\n", + "min_display_ndre = np.percentile(masked_ndre.compressed(), 5)\n", + "max_display_ndre = np.percentile(masked_ndre.compressed(), 99.5)\n", "\n", "fig, axis = plotutils.plot_overlay_withcolorbar(\n", " gamma_corr_rgb,\n", diff --git a/Alignment-RigRelatives.ipynb b/Alignment-RigRelatives.ipynb index f775312..02db654 100644 --- a/Alignment-RigRelatives.ipynb +++ b/Alignment-RigRelatives.ipynb @@ -372,9 +372,9 @@ "plt.show()\n", "\n", "min_display_ndvi = 0.45 # further mask soil by removing low-ndvi values\n", - "# min_display_ndvi = np.percentile(ndvi.flatten(), 5.0) # modify with these percentilse to adjust contrast\n", + "# min_display_ndvi = np.percentile(ndvi.compressed(), 5.0) # modify with these percentilse to adjust contrast\n", "max_display_ndvi = np.percentile(\n", - " ndvi.flatten(), 99.5\n", + " ndvi.compressed(), 99.5\n", ") # for many images, 0.5 and 99.5 are good values\n", "masked_ndvi = np.ma.masked_where(ndvi < min_display_ndvi, ndvi)\n", "\n", @@ -426,8 +426,8 @@ "plt.title(\"NDRE Histogram (filtered to only plants)\")\n", "plt.show()\n", "\n", - "min_display_ndre = np.percentile(masked_ndre, 5)\n", - "max_display_ndre = np.percentile(masked_ndre, 99.5)\n", + "min_display_ndre = np.percentile(masked_ndre.compressed(), 5)\n", + "max_display_ndre = np.percentile(masked_ndre.compressed(), 99.5)\n", "\n", "fig, axis = plotutils.plot_overlay_withcolorbar(\n", " gamma_corr_rgb,\n", diff --git a/Alignment.ipynb b/Alignment.ipynb index 82b2dc7..35aa509 100644 --- a/Alignment.ipynb +++ b/Alignment.ipynb @@ -422,9 +422,9 @@ "plt.show()\n", "\n", "min_display_ndvi = 0.45 # further mask soil by removing low-ndvi values\n", - "# min_display_ndvi = np.percentile(ndvi.flatten(), 5.0) # modify with these percentilse to adjust contrast\n", + "# min_display_ndvi = np.percentile(ndvi.compressed(), 5.0) # modify with these percentilse to adjust contrast\n", "max_display_ndvi = np.percentile(\n", - " ndvi.flatten(), 99.5\n", + " ndvi.compressed(), 99.5\n", ") # for many images, 0.5 and 99.5 are good values\n", "masked_ndvi = np.ma.masked_where(ndvi < min_display_ndvi, ndvi)\n", "\n", @@ -475,8 +475,8 @@ "plt.title(\"NDRE Histogram (filtered to only plants)\")\n", "plt.show()\n", "\n", - "min_display_ndre = np.percentile(masked_ndre, 5)\n", - "max_display_ndre = np.percentile(masked_ndre, 99.5)\n", + "min_display_ndre = np.percentile(masked_ndre.compressed(), 5)\n", + "max_display_ndre = np.percentile(masked_ndre.compressed(), 99.5)\n", "\n", "fig, axis = plotutils.plot_overlay_withcolorbar(\n", " gamma_corr_rgb,\n", diff --git a/docs/Alignment v2.html b/docs/Alignment v2.html deleted file mode 100644 index 004c9a5..0000000 --- a/docs/Alignment v2.html +++ /dev/null @@ -1,16037 +0,0 @@ - - -
- - -For most use cases, each band of a multispectral capture must be aligned with the other bands in order to create meaningful data. In this tutorial, we show how to align the band to each other using open source OpenCV utilities.
-Image alignment allows the combination of images into true-color (RGB) and false color (such as CIR) composites, useful for scouting using single images as well as for display and management uses. In addition to composite images, alignment allows the calculation of pixel-accurate indices such as NDVI or NDRE at the single image level which can be very useful for applications like plant counting and coverage estimations, where mosaicing artifacts may otherwise skew analysis results.
-The image alignment method described below tends to work well on images with abundant image features, or areas of significant contrast. Cars, buildings, parking lots, and roads tend to provide the best results. This approach may not work well on images which contain few features or very repetitive features, such as full canopy row crops or fields of repetitive small crops such lettuce or strawberries. We will disscuss more about the advantages and disadvantages of these methods below.
-The functions behind this alignment process can work with most versions of RedEdge and Altum firmware. They will work best with RedEdge (3,M,MX) versions above 3.2.0 which include the "RigRelatives" tags, and all RedEdge-P/Altum/Altum-PT imagery. These tags provide a starting point for the image transformation and can help to ensure convergence of the algorithm.
-As we have done in previous examples, we use the micasense.capture class to open, radiometrically correct, and visualize all the bands of a MicaSense capture.
-First, we'll load the autoreload extension. This lets us change underlying code (such as library functions) without having to reload the entire workbook and kernel. This is useful in this workbook because the cell that runs the alignment can take a long time to run, so with the autoreload extension we can update the code after the alignment step for analysis and visualization without needing to re-compute the alignments each time.
%load_ext autoreload
-%autoreload 2
-import os, glob
-import micasense.capture as capture
-%matplotlib inline
-from pathlib import Path
-import matplotlib.pyplot as plt
-plt.rcParams["figure.facecolor"] = "w"
-
-panelNames = None
-
-# set your image paths here. See more here: https://docs.python.org/3/library/pathlib.html
-# if using Windows, you need to an an "r" to the path like this: Path(r"C:\Files")
-
-# imagePath = Path("./data/REDEDGE-MX-DUAL")
-
-# # these will return lists of image paths as strings
-# imageNames = list(imagePath.glob('IMG_0007_*.tif'))
-# imageNames = [x.as_posix() for x in imageNames]
-
-# panelNames = list(imagePath.glob('IMG_0001_*.tif'))
-# panelNames = [x.as_posix() for x in panelNames]
-
-# imagePath = Path("./data/REDEDGE-MX")
-
-# # these will return lists of image paths as strings
-# imageNames = list(imagePath.glob('IMG_0020_*.tif'))
-# imageNames = [x.as_posix() for x in imageNames]
-
-# panelNames = list(imagePath.glob('IMG_0001_*.tif'))
-# panelNames = [x.as_posix() for x in panelNames]
-
-imagePath = Path("./data/ALTUM-PT")
-
-# these will return lists of image paths as strings
-imageNames = list(imagePath.glob('IMG_0010_*.tif'))
-imageNames = [x.as_posix() for x in imageNames]
-
-panelNames = list(imagePath.glob('IMG_0000_*.tif'))
-panelNames = [x.as_posix() for x in panelNames]
-
-# imagePath = Path("./data/ALTUM")
-
-# # these will return lists of image paths as strings
-# imageNames = list(imagePath.glob('IMG_0021_*.tif'))
-# imageNames = [x.as_posix() for x in imageNames]
-
-# panelNames = list(imagePath.glob('IMG_0000_*.tif'))
-# panelNames = [x.as_posix() for x in panelNames]
-
-# imagePath = Path("./data/REDEDGE-P")
-
-# # these will return lists of image paths as strings
-# imageNames = list(imagePath.glob('IMG_0011_*.tif'))
-# imageNames = [x.as_posix() for x in imageNames]
-
-# panelNames = list(imagePath.glob('IMG_0000_*.tif'))
-# panelNames = [x.as_posix() for x in panelNames]
-
-
-if panelNames is not None:
- panelCap = capture.Capture.from_filelist(panelNames)
-else:
- panelCap = None
-
-thecapture = capture.Capture.from_filelist(imageNames)
-
-# get camera model for future use
-cam_model = thecapture.camera_model
-# if this is a multicamera system like the RedEdge-MX Dual,
-# we can combine the two serial numbers to help identify
-# this camera system later.
-if len(thecapture.camera_serials) > 1:
- cam_serial = "_".join(thecapture.camera_serials)
- print(cam_serial)
-else:
- cam_serial = thecapture.camera_serial
-
-print("Camera model:",cam_model)
-print("Bit depth:", thecapture.bits_per_pixel)
-print("Camera serial number:", cam_serial)
-print("Capture ID:",thecapture.uuid)
-
-# determine if this sensor has a panchromatic band
-if cam_model == 'RedEdge-P' or cam_model == 'Altum-PT':
- panchroCam = True
-else:
- panchroCam = False
- panSharpen = False
-
-if panelCap is not None:
- if panelCap.panel_albedo() is not None:
- panel_reflectance_by_band = panelCap.panel_albedo()
- else:
- panel_reflectance_by_band = [0.49]*len(thecapture.eo_band_names()) #RedEdge band_index order
- panel_irradiance = panelCap.panel_irradiance(panel_reflectance_by_band)
- img_type = "reflectance"
- thecapture.plot_undistorted_reflectance(panel_irradiance)
-else:
- if thecapture.dls_present():
- img_type='reflectance'
- thecapture.plot_undistorted_reflectance(thecapture.dls_irradiance())
- else:
- img_type = "radiance"
- thecapture.plot_undistorted_radiance()
-Camera model: Altum-PT -Bit depth: 12 -Camera serial number: PA01-2210071-MS -Capture ID: Z9WKG9xZ2rL9gFGQopS0 --
If we have already successfully aligned captures from this specific camera, we can typically save some time and use the alignment warp matrices for other captures from the same camera
- -from skimage.transform import ProjectiveTransform
-import numpy as np
-
-if panchroCam:
- warp_matrices_filename = cam_serial + "_warp_matrices_SIFT.npy"
-else:
- warp_matrices_filename = cam_serial + "_warp_matrices_opencv.npy"
-
-if Path('./' + warp_matrices_filename).is_file():
- print("Found existing warp matrices for camera", cam_serial)
- load_warp_matrices = np.load(warp_matrices_filename, allow_pickle=True)
- loaded_warp_matrices = []
- for matrix in load_warp_matrices:
- if panchroCam:
- transform = ProjectiveTransform(matrix=matrix.astype('float64'))
- loaded_warp_matrices.append(transform)
- else:
- loaded_warp_matrices.append(matrix.astype('float32'))
- print("Warp matrices successfully loaded.")
-
- if panchroCam:
- warp_matrices_SIFT = loaded_warp_matrices
- else:
- warp_matrices = loaded_warp_matrices
-else:
- print("No existing warp matrices found. Create them later in the notebook.")
- warp_matrices_SIFT = False
- warp_matrices = False
-Found existing warp matrices for camera PA01-2210071-MS -Warp matrices successfully loaded. --
Alignment is a three-step process:
-Images are unwarped using the built-in lens calibration -A transformation is found to align each band to a common band -The aligned images are combined and cropped, removing pixels which don't overlap in all bands. -We provide utilities to find the alignment transformations within a single capture. Our experience shows that once a good alignment transformation is found, it tends to be stable over a flight and, in most cases, over many flights. The transformation may change if the camera undergoes a shock event (such as a hard landing or drop) or if the temperature changes substantially between flights. In these events, a new transformation may need to be found.
-Further, since this approach finds a 2-dimensional (affine) transformation between images, it won't work when the parallax between bands results in a 3-dimensional depth field. This can happen if very close to the target or when targets are visible at significantly different ranges, such as a nearby tree or building against a background much farther way. In these cases, it will be necessary to use photogrammetry techniques to find a 3-dimensional mapping between images.
-For best alignment results it's good to select a capture which has features which visible in all bands. Man-made objects such as cars, roads, and buildings tend to work very well, while captures of only repeating crop rows tend to work poorly. Remember, once a good transformation has been found for flight, it can be generally be applied across all of the images.
-It's also good to use an image for alignment which is taken near the same level above ground as the rest of the flights. Above approximately 35m AGL, the alignment will be consistent. However, if images taken closer to the ground are used, such as panel images, the same alignment transformation will not work for the flight data.
- -import cv2
-import time
-import numpy as np
-import matplotlib.pyplot as plt
-import micasense.imageutils as imageutils
-import micasense.plotutils as plotutils
-
-# We use a different alignment method for RedEdge-P and Altum-PT,
-# so if the imagery is from this kind of camera, skip this step
-if not panchroCam:
- st = time.time()
- # set to True if you'd like to ignore existing warp matrices and create new ones
- regenerate = True
- pyramid_levels = 0 # for images with RigRelatives, setting this to 0 or 1 may improve alignment
- max_alignment_iterations = 10
-
- # match_index:
- # for non-panchromatic cameras we want to use band 1, which is green
- # the green band has zero rig relative offsets
- # NOTE: These band numbers are zero-indexed
- # special parameters for RedEdge-MX dual camera system
- if len(thecapture.eo_band_names()) == 10:
- print("is 10 band")
- pyramid_levels = 3
- match_index = 4
- max_alignment_iterations = 20
- else:
- match_index = 1
-
- warp_mode = cv2.MOTION_HOMOGRAPHY # MOTION_HOMOGRAPHY or MOTION_AFFINE. For Altum images only use HOMOGRAPHY
-
- print("Aligning images. Depending on settings this can take from a few seconds to many minutes")
- # Can potentially increase max_iterations for better results, but longer runtimes
- if warp_matrices and not regenerate:
- print("Using existing warp matrices...")
- try:
- irradiance = panel_irradiance+[0]
- except NameError:
- irradiance = None
- else:
- warp_matrices, alignment_pairs = imageutils.align_capture(thecapture,
- ref_index = match_index,
- max_iterations = max_alignment_iterations,
- warp_mode = warp_mode,
- pyramid_levels = pyramid_levels)
-
- print("Finished Aligning")
- et = time.time()
- elapsed_time = et - st
- print('Alignment time:', int(elapsed_time), 'seconds')
-After finding image alignments, we may need to remove pixels around the edges which aren't present in every image in the capture. To do this we use the affine transforms found above and the image distortions from the image metadata. OpenCV provides a couple of handy helpers for this task in the cv2.undistortPoints() and cv2.transform() methods. These methods take a set of pixel coordinates and apply our undistortion matrix and our affine transform, respectively. So, just as we did when registering the images, we first apply the undistortion process to the coordinates of the image borders, then we apply the affine transformation to that result. The resulting pixel coordinates tell us where the image borders end up after this pair of transformations, and we can then crop the resultant image to these coordinates.
- -if not panchroCam:
- cropped_dimensions, edges = imageutils.find_crop_bounds(thecapture, warp_matrices, warp_mode=warp_mode, reference_band=match_index)
- print("Cropped dimensions:",cropped_dimensions)
- im_aligned = thecapture.create_aligned_capture(warp_matrices=warp_matrices, motion_type=warp_mode, img_type=img_type, match_index=match_index)
-For older cameras we use OpenCV for capture alignment. For RedEdge-P and Altum-PT, we use SIFT (scale-invariant feature transform). We will use SIFT to create the warp matrices, then use these warp matrices to align the capture. See more here: -https://scikit-image.org/docs/stable/auto_examples/features_detection/plot_sift.html -https://en.wikipedia.org/wiki/Scale-invariant_feature_transform
-For sensors with a panchromatic band (Altum-PT or RedEdge-P), we may wish to create a pan-sharpened stack. This example uses a linear interpolation method for pan-sharpening.
-The radiometric_pan_sharpen function takes in the SIFT warp matrices and outputs an upsampled image stack (all bands are changed to the resolution of the panchromatic sensor) and a pan-sharpened stack.
Note: it is important to choose an initial alignment image that has a lot of straight, manmade features, such as roads or buildings. Trying to align a capture that is strictly vegetation will take a long time and may not be very good. This process can take a while, depending on how feature-rich the capture is. We have seen some Altum-PT captures align after 2 minutes, and some can take upwards of 30 minutes. Be patient! Once an alignment has been found, it can be used on other captures from the same camera, and the alignment process will be much faster.
- -# from micasense.imageutils import brovey_pan_sharpen,radiometric_pan_sharpen
-from skimage.transform import ProjectiveTransform
-
-import time
-
-if panchroCam:
- # set to True if you'd like to ignore existing warp matrices and create new ones
- regenerate = False
- st = time.time()
- if not warp_matrices_SIFT or regenerate:
- print("Generating new warp matrices...")
- warp_matrices_SIFT = thecapture.SIFT_align_capture(min_matches = 10)
-
- sharpened_stack, upsampled = thecapture.radiometric_pan_sharpened_aligned_capture(warp_matrices=warp_matrices_SIFT)
-
-# we can also use the Rig Relatives from the image metadata to do a quick, rudimentary alignment
-# warp_matrices0=thecapture.get_warp_matrices(ref_index=5)
-# sharpened_stack,upsampled = radiometric_pan_sharpen(thecapture,warp_matrices=warp_matrices0)
-
- print("Pansharpened shape:", sharpened_stack.shape)
- print("Upsampled shape:", upsampled.shape)
- # re-assign to im_aligned to match rest of code
- im_aligned = upsampled
- et = time.time()
- elapsed_time = et - st
- print('Alignment and pan-sharpening time:', int(elapsed_time), 'seconds')
-Pansharpened shape: (2990, 3980, 7) -Upsampled shape: (2990, 3980, 7) -Alignment and pan-sharpening time: 18 seconds --
Once an alignment for your camera has been found, it can be saved to a file for later use with this notebook. It can also be used on the Batch Alignment notebook for aligning all of the captures from an entire flight.
- -import numpy
-import skimage
-from skimage.transform import warp,matrix_transform,resize,FundamentalMatrixTransform,estimate_transform,ProjectiveTransform
-
-if panchroCam:
- working_wm = warp_matrices_SIFT
-else:
- working_wm = warp_matrices
-if not Path('./' + warp_matrices_filename).is_file() or regenerate:
- temp_matrices = []
- for x in working_wm:
- if isinstance(x, numpy.ndarray):
- temp_matrices.append(x)
- if isinstance(x, skimage.transform._geometric.ProjectiveTransform):
- temp_matrices.append(x.params)
- np.save(warp_matrices_filename, np.array(temp_matrices, dtype=object), allow_pickle=True)
- print("Saved to", Path('./' + warp_matrices_filename).resolve())
-else:
- print("Matrices already exist at",Path('./' + warp_matrices_filename).resolve())
-Matrices already exist at /Users/stephen/programs/gitlab/imageProcessingTutorial/PA01-2210071-MS_warp_matrices_SIFT.npy --
We can compare the radiance between multispectral bands, and between upsampled/pansharpened in the case of RedEdge-P and Altum-PT
- -theColors = {'Blue': 'blue', 'Green': 'green', 'Red': 'red', \
- 'Red edge': 'maroon', 'NIR': 'purple', 'Panchro': 'yellow', \
- 'Red edge-740': 'salmon', 'Red Edge': 'maroon', 'Blue-444': 'aqua', \
- 'Green-531': 'lime', 'Red-650': 'lightcoral', 'Red edge-705':'brown'}
-
-eo_count = len(thecapture.eo_indices())
-multispec_min = np.min(np.percentile(im_aligned[:,:,1:eo_count].flatten(),0.01))
-multispec_max = np.max(np.percentile(im_aligned[:,:,1:eo_count].flatten(), 99.99))
-
-theRange = (multispec_min,multispec_max)
-
-fig, axis = plt.subplots(1, 1, figsize=(10,4))
-for x,y in zip(thecapture.eo_indices(),thecapture.eo_band_names()):
- axis.hist(im_aligned[:,:,x].ravel(), bins=512, range=theRange, \
- histtype="step", label=y, color=theColors[y], linewidth=1.5)
-plt.title("Multispectral histogram (radiance)")
-axis.legend()
-plt.show()
-
-if panchroCam:
- eo_count = len(thecapture.eo_indices())
- multispec_min = np.min(np.percentile(sharpened_stack[:,:,1:eo_count].flatten(),0.01))
- multispec_max = np.max(np.percentile(sharpened_stack[:,:,1:eo_count].flatten(), 99.99))
-
- theRange = (multispec_min,multispec_max)
-
- fig, axis = plt.subplots(1, 1, figsize=(10,4))
- for x,y in zip(thecapture.eo_indices(),thecapture.eo_band_names()):
- axis.hist(sharpened_stack[:,:,x].ravel(), bins=512, range=theRange, \
- histtype="step", label=y, color=theColors[y], linewidth=1.5)
- plt.title("Pan-sharpened multispectral histogram (radiance)")
- axis.legend()
- plt.show()
-Once the transformation has been found, it can be verified by compositing the aligned images to check alignment. The image 'stack' containing all bands can also be exported to a multi-band TIFF file for viewing in external software such as QGIS. Useful composites are a naturally colored RGB as well as color infrared, or CIR.
- -# figsize=(30,23) # use this size for full-image-resolution display
-figsize=(16,13) # use this size for export-sized display
-
-rgb_band_indices = [thecapture.band_names_lower().index('red'),
- thecapture.band_names_lower().index('green'),
- thecapture.band_names_lower().index('blue')]
-cir_band_indices = [thecapture.band_names_lower().index('nir'),
- thecapture.band_names_lower().index('red'),
- thecapture.band_names_lower().index('green')]
-
-# Create normalized stacks for viewing
-im_display = np.zeros((im_aligned.shape[0],im_aligned.shape[1],im_aligned.shape[2]), dtype=np.float32)
-im_min = np.percentile(im_aligned[:,:,rgb_band_indices].flatten(), 0.5) # modify these percentiles to adjust contrast
-im_max = np.percentile(im_aligned[:,:,rgb_band_indices].flatten(), 99.5) # for many images, 0.5 and 99.5 are good values
-
-if panchroCam:
- im_display_sharp = np.zeros((sharpened_stack.shape[0],sharpened_stack.shape[1],sharpened_stack.shape[2]), dtype=np.float32 )
- im_min_sharp = np.percentile(sharpened_stack[:,:,rgb_band_indices].flatten(), 0.5) # modify these percentiles to adjust contrast
- im_max_sharp = np.percentile(sharpened_stack[:,:,rgb_band_indices].flatten(), 99.5) # for many images, 0.5 and 99.5 are good values
-
-
-# for rgb true color, we use the same min and max scaling across the 3 bands to
-# maintain the "white balance" of the calibrated image
-for i in rgb_band_indices:
- im_display[:,:,i] = imageutils.normalize(im_aligned[:,:,i], im_min, im_max)
- if panchroCam:
- im_display_sharp[:,:,i] = imageutils.normalize(sharpened_stack[:,:,i], im_min_sharp, im_max_sharp)
-
-rgb = im_display[:,:,rgb_band_indices]
-
-if panchroCam:
- rgb_sharp = im_display_sharp[:,:,rgb_band_indices]
-
-nir_band = thecapture.band_names_lower().index('nir')
-red_band = thecapture.band_names_lower().index('red')
-
-ndvi = (im_aligned[:,:,nir_band] - im_aligned[:,:,red_band]) / (im_aligned[:,:,nir_band] + im_aligned[:,:,red_band])
-
-# for cir false color imagery, we normalize the NIR,R,G bands within themselves, which provides
-# the classical CIR rendering where plants are red and soil takes on a blue tint
-for i in cir_band_indices:
- im_display[:,:,i] = imageutils.normalize(im_aligned[:,:,i])
-
-cir = im_display[:,:,cir_band_indices]
-if panchroCam:
- fig, (ax1,ax2) = plt.subplots(1, 2, figsize=figsize)
-else:
- fig, ax1 = plt.subplots(1, 1, figsize=figsize)
-ax1.set_title("Red-Green-Blue Composite")
-ax1.imshow(rgb)
-if panchroCam:
- ax2.set_title("Red-Green-Blue Composite (pan-sharpened)")
- ax2.imshow(rgb_sharp)
-
-fig, (ax3,ax4) = plt.subplots(1, 2, figsize=figsize)
-ax3.set_title("NDVI")
-ax3.imshow(ndvi)
-ax4.set_title("Color Infrared (CIR) Composite")
-ax4.imshow(cir)
-
-# set custom lims if you want to zoom in to image to see more detail
-# this is useful for comparing upsampled and pan-sharpened stacks
-# custom_xlim=(1500,2000)
-# custom_ylim=(2000,1500)
-# plt.setp([ax1,ax2,ax3,ax4], xlim=custom_xlim, ylim=custom_ylim)
-
-plt.show()
-There are many techniques for image enhancement, but one which is commonly used to improve the visual sharpness of imagery is the unsharp mask. Here, we apply an unsharp mask to the RGB image to improve the visualization, and then apply a gamma curve to make the darkest areas brighter.
- -if panchroCam:
- rgb = rgb_sharp
-# Create an enhanced version of the RGB render using an unsharp mask
-gaussian_rgb = cv2.GaussianBlur(rgb, (9,9), 10.0)
-gaussian_rgb[gaussian_rgb<0] = 0
-gaussian_rgb[gaussian_rgb>1] = 1
-unsharp_rgb = cv2.addWeighted(rgb, 1.5, gaussian_rgb, -0.5, 0)
-unsharp_rgb[unsharp_rgb<0] = 0
-unsharp_rgb[unsharp_rgb>1] = 1
-
-# Apply a gamma correction to make the render appear closer to what our eyes would see
-gamma = 1.4
-gamma_corr_rgb = unsharp_rgb**(1.0/gamma)
-fig = plt.figure(figsize=figsize)
-plt.imshow(gamma_corr_rgb, aspect='equal')
-plt.axis('off')
-plt.show()
-We can easily export the stacks into an image using the GDAL library (http://www.glal.org). Once exported, these image stacks can be opened in software such as QGIS and raster operations such as NDVI or NDRE computation can be done in that software. The stacks include geographic information.
-If you prefer, you may set sort_by_wavelength to True in the save_capture_as_stack function.
Unless otherwise specified, this will save in your working folder, that is your imageprocessing directory.
#set output name to unique capture ID, e.g. FWoNSvgDNBX63Xv378qs
-outputName = thecapture.uuid
-
-try:
- irradiance = panel_irradiance+[0]
-except NameError:
- irradiance = None
-
-st = time.time()
-if panchroCam:
- # in this example, we can export both a pan-sharpened stack and an upsampled stack
- # so you can compare them in GIS. In practice, you would typically only output the pansharpened stack
- thecapture.save_capture_as_stack(outputName+"-pansharpened.tif", sort_by_wavelength=False,pansharpen=True,img_type=img_type)
- thecapture.save_capture_as_stack(outputName+"-upsampled.tif", sort_by_wavelength=False,pansharpen=False,img_type=img_type)
-else:
- thecapture.save_capture_as_stack(outputName+".tif", sort_by_wavelength=True,img_type=img_type)
-
-et = time.time()
-elapsed_time = et - st
-print("Time to save stacks:", int(elapsed_time), "seconds.")
-Time to save stacks: 34 seconds. --
For raw index computation on single images, the numpy package provides a simple way to do math and simple visualization on images. Below, we compute and visualize an image histogram, and then use that to pick a color map range for visualizing the NDVI of an image.
After computing the NDVI and prior to displaying it, we use a very rudimentary method for focusing on the plants and removing the soil and shadow information from our images and histograms. Below, we remove non-plant pixels by setting to zero any pixels in the image where the NIR reflectance is less than 20%. This helps to ensure that the NDVI and NDRE histograms aren't skewed substantially by soil noise.
- -from micasense import plotutils
-import matplotlib.pyplot as plt
-
-nir_band = thecapture.band_names_lower().index('nir')
-red_band = thecapture.band_names_lower().index('red')
-
-thelayer = im_aligned
-if panchroCam:
- thelayer = sharpened_stack
-np.seterr(divide='ignore', invalid='ignore') # ignore divide by zero errors in the index calculation
-
-# Compute Normalized Difference Vegetation Index (NDVI) from the NIR(3) and RED (2) bands
-ndvi = (thelayer[:,:,nir_band] - thelayer[:,:,red_band]) / (thelayer[:,:,nir_band] + thelayer[:,:,red_band])
-print("Image type:",img_type)
-
-# remove shadowed areas (mask pixels with NIR reflectance < 20%))
-# this does not seem to work on panchro stacks
-if img_type == 'reflectance':
- ndvi = np.ma.masked_where(thelayer[:,:,nir_band] < 0.20, ndvi)
-elif img_type == 'radiance':
- lower_pct_radiance = np.percentile(thelayer[:,:,nir_band], 10.0)
- ndvi = np.ma.masked_where(thelayer[:,:,nir_band] < lower_pct_radiance, ndvi)
-# Compute and display a histogram
-# ndvi_hist_min = np.min(ndvi)
-# ndvi_hist_max = np.max(ndvi)
-ndvi_hist_min = np.min(np.percentile(ndvi,0.5))
-ndvi_hist_max = np.max(np.percentile(ndvi,99.5))
-fig, axis = plt.subplots(1, 1, figsize=(10,4))
-axis.hist(ndvi.ravel(), bins=512, range=(ndvi_hist_min, ndvi_hist_max))
-plt.title("NDVI Histogram")
-plt.show()
-
-min_display_ndvi = 0.45 # further mask soil by removing low-ndvi values
-#min_display_ndvi = np.percentile(ndvi.flatten(), 5.0) # modify with these percentilse to adjust contrast
-max_display_ndvi = np.percentile(ndvi.flatten(), 99.5) # for many images, 0.5 and 99.5 are good values
-masked_ndvi = np.ma.masked_where(ndvi < min_display_ndvi, ndvi)
-
-#reduce the figure size to account for colorbar
-figsize=np.asarray(figsize) - np.array([3,2])
-
-#plot NDVI over an RGB basemap, with a colorbar showing the NDVI scale
-fig, axis = plotutils.plot_overlay_withcolorbar(gamma_corr_rgb,
- masked_ndvi,
- figsize = (14,7),
- title = 'NDVI filtered to only plants over RGB base layer',
- vmin = min_display_ndvi,
- vmax = max_display_ndvi)
-fig.savefig(thecapture.uuid+'_ndvi_over_rgb.png')
-Image type: reflectance --
In the same manner, we can compute, filter, and display another index useful for MicaSense cameras, the Normalized Difference Red Edge (NDRE) index. We also filter out shadows and soil to ensure our display focuses only on the plant health.
- -# Compute Normalized Difference Red Edge Index from the NIR(3) and RedEdge(4) bands
-rededge_band = thecapture.band_names_lower().index('red edge')
-ndre = (thelayer[:,:,nir_band] - thelayer[:,:,rededge_band]) / (thelayer[:,:,nir_band] + thelayer[:,:,rededge_band])
-
-# Mask areas with shadows and low NDVI to remove soil
-masked_ndre = np.ma.masked_where(ndvi < min_display_ndvi, ndre)
-
-# Compute a histogram
-ndre_hist_min = np.min(np.percentile(masked_ndre,0.5))
-ndre_hist_max = np.max(np.percentile(masked_ndre,99.5))
-fig, axis = plt.subplots(1, 1, figsize=(10,4))
-axis.hist(masked_ndre.ravel(), bins=512, range=(ndre_hist_min, ndre_hist_max))
-plt.title("NDRE Histogram (filtered to only plants)")
-plt.show()
-
-min_display_ndre = np.percentile(masked_ndre, 5)
-max_display_ndre = np.percentile(masked_ndre, 99.5)
-
-fig, axis = plotutils.plot_overlay_withcolorbar(gamma_corr_rgb,
- masked_ndre,
- figsize=(14,7),
- title='NDRE filtered to only plants over RGB base layer',
- vmin=min_display_ndre,vmax=max_display_ndre)
-fig.savefig(thecapture.uuid+'_ndre_over_rgb.png')
-If our image is from an Altum or Altum-PT and includes a thermal band, we can display the re-sampled and aligned thermal data over the RGB data to maintain the context of the thermal information.
- -if len(thecapture.lw_indices()) > 0:
- lwir_band = thecapture.band_names_lower().index('lwir')
-
- # by default we don't mask the thermal, since it's native resolution is much lower than the MS
- if panchroCam:
- masked_thermal = sharpened_stack[:,:,lwir_band]
- else:
- masked_thermal = im_aligned[:,:,lwir_band]
- # Alternatively we can mask the thermal only to plants here, which is useful for large contiguous areas
- # masked_thermal = np.ma.masked_where(ndvi < 0.45, im_aligned[:,:,5])
-
-
- # Compute a histogram
- fig, axis = plt.subplots(1, 1, figsize=(10,4))
- axis.hist(masked_thermal.ravel(), bins=512, range=(np.min(np.percentile(masked_thermal,1)), np.max(np.percentile(masked_thermal,99))))
- plt.title("Thermal Histogram")
- plt.show()
-
- min_display_therm = np.percentile(masked_thermal, 1)
- max_display_therm = np.percentile(masked_thermal, 99)
-
- fig, axis = plotutils.plot_overlay_withcolorbar(gamma_corr_rgb,
- masked_thermal,
- figsize=(14,7),
- title='Temperature over True Color',
- vmin=min_display_therm,vmax=max_display_therm,
- overlay_alpha=0.25,
- overlay_colormap='jet',
- overlay_steps=16,
- display_contours=True,
- contour_steps=16,
- contour_alpha=.4,
- contour_fmt="%.0fC")
- fig.savefig(thecapture.uuid+'_thermal_over_rgb.png')
-Finally, we show a classic agricultural remote sensing output in the tassled cap plot. This plot can be useful for visualizing row crops and plots the Red Reflectance channel on the X-axis against the NIR reflectance channel on the Y-axis. This plot also clearly shows the line of the soil in that space. The tassled cap view isn't very useful for this arid data set; however, we can see the "badge of trees" of high NIR reflectance and relatively low red reflectance. This provides an example of one of the uses of aligned images for single capture analysis.
- -x_band = red_band
-y_band = nir_band
-x_max = np.max(np.percentile(im_aligned[:,:,x_band],99.99))
-y_max = np.max(np.percentile(im_aligned[:,:,y_band],99.99))
-
-fig = plt.figure(figsize=(12,12))
-plt.hexbin(im_aligned[:,:,x_band],im_aligned[:,:,y_band],gridsize=640,bins='log',extent=(0,x_max,0,y_max))
-ax = fig.gca()
-ax.set_xlim([0,x_max])
-ax.set_ylim([0,y_max])
-plt.xlabel("{} Reflectance".format(thecapture.band_names()[x_band]))
-plt.ylabel("{} Reflectance".format(thecapture.band_names()[y_band]))
-plt.show()
-
-if panchroCam:
- x_band = red_band
- y_band = nir_band
- x_max = np.max(np.percentile(sharpened_stack[:,:,x_band],99.99))
- y_max = np.max(np.percentile(sharpened_stack[:,:,y_band],99.99))
-
- fig = plt.figure(figsize=(12,12))
- plt.hexbin(sharpened_stack[:,:,x_band],sharpened_stack[:,:,y_band],gridsize=640,bins='log',extent=(0,x_max,0,y_max))
- ax = fig.gca()
- ax.set_xlim([0,x_max])
- ax.set_ylim([0,y_max])
- plt.xlabel("{} Reflectance (pan-sharpened)".format(thecapture.band_names()[x_band]))
- plt.ylabel("{} Reflectance (pan-sharpened)".format(thecapture.band_names()[y_band]))
- plt.show()
-Lastly, we output the warp_matrices that we got for this image stack for usage elsewhere. Currently, these can be used in the Batch Processing.ipynb notebook to save reflectance-compensated stacks of images to a directory.
if panchroCam:
- print(warp_matrices_SIFT)
-else:
- print(warp_matrices)
-[<ProjectiveTransform(matrix= - [[ 4.71051427e-01, 1.81635696e-03, 6.30150447e+01], - [-2.91181288e-03, 4.69196185e-01, 1.02276511e+02], - [ 5.38572502e-07, -1.51241146e-06, 1.00000000e+00]]) at 0x7f9109f44410>, <ProjectiveTransform(matrix= - [[ 4.68668790e-01, 9.34580148e-04, 5.84716582e+01], - [-2.23630030e-03, 4.67555112e-01, 9.71512698e+01], - [ 2.90282499e-07, -1.34188378e-06, 1.00000000e+00]]) at 0x7f9109f44110>, <ProjectiveTransform(matrix= - [[ 4.70308622e-01, 1.86821829e-03, 5.15342645e+01], - [-2.12534374e-03, 4.70189903e-01, 9.11463269e+01], - [-8.44058028e-08, -1.49050229e-07, 1.00000000e+00]]) at 0x7f9109f446d0>, <ProjectiveTransform(matrix= - [[ 4.67756541e-01, 2.58999879e-03, 6.22011674e+01], - [-3.30003520e-03, 4.66852334e-01, 9.37899796e+01], - [ 1.67165473e-07, -8.40910501e-07, 1.00000000e+00]]) at 0x7f9109f441d0>, <ProjectiveTransform(matrix= - [[ 4.68418989e-01, 1.79902364e-03, 5.07127933e+01], - [-3.26474577e-03, 4.68172796e-01, 9.85409805e+01], - [-6.14928956e-07, -9.35021962e-07, 1.00000000e+00]]) at 0x7f9109f44490>, <ProjectiveTransform(matrix= - [[1., 0., 0.], - [0., 1., 0.], - [0., 0., 1.]]) at 0x7f908f2e94d0>, <ProjectiveTransform(matrix= - [[ 6.63024562e-02, 3.46707043e-04, -9.35690463e-01], - [-6.41765998e-04, 6.52823865e-02, 2.98242993e+01], - [-2.09231083e-06, -3.01713640e-07, 1.00000000e+00]]) at 0x7f908f2fa810>] --
Note: This notebook specifically demonstrates alignment and stack creation for 10-band data. For alignment of 5-band (RedEdge) and 6-band (Altum) data, pr more details about alignment in general, see the "Active Image Alignment" notebook.
-As we have done in previous examples, we use the micasense.capture class to open, radiometrically correct, and visualize the 10 bands of a RedEdge-MX Dual (Red+Blue) capture.
-First, we'll load the autoreload extension. This lets us change underlying code (such as library functions) without having to reload the entire workbook and kernel. This is useful in this workbook because the cell that runs the alignment can take a long time to run, so with autoreload extension we can update the code after the alignment step for analysis and visualization without needing to re-compute the alignments each time.
%load_ext autoreload
-%autoreload 2
-import os, glob
-import micasense.capture as capture
-%matplotlib inline
-
-panelNames = None
-paneCap = None
-
-imagePath = os.path.join('.','data','REDEDGE-MX-DUAL')
-imageNames = glob.glob(os.path.join(imagePath,'IMG_0007_*.tif'))
-panelNames = glob.glob(os.path.join(imagePath,'IMG_0001_*.tif'))
-
-# Allow this code to align both radiance and reflectance images; bu excluding
-# a definition for panelNames above, radiance images will be used
-# For panel images, efforts will be made to automatically extract the panel information
-# but if the panel/firmware is before Altum 1.3.5, RedEdge 5.1.7 the panel reflectance
-# will need to be set in the panel_reflectance_by_band variable.
-# Note: radiance images will not be used to properly create NDVI/NDRE images below.
-if panelNames is not None:
- panelCap = capture.Capture.from_filelist(panelNames)
-else:
- panelCap = None
-
-capture = capture.Capture.from_filelist(imageNames)
-
-if panelCap is not None:
- if panelCap.panel_albedo() is not None:
- panel_reflectance_by_band = panelCap.panel_albedo()
- else:
- raise IOError("Comment this lne and set panel_reflectance_by_band here")
- panel_reflectance_by_band = [0.65]*len(imageNames)
- panel_irradiance = panelCap.panel_irradiance(panel_reflectance_by_band)
- img_type = "reflectance"
- capture.plot_undistorted_reflectance(panel_irradiance)
-else:
- if capture.dls_present():
- img_type='reflectance'
- capture.plot_undistorted_reflectance(capture.dls_irradiance())
- else:
- img_type = "radiance"
- capture.plot_undistorted_radiance()
-Alignment is a three step process:
-We provide utilities to find the alignement transformations within a single capture. Our experience shows that once a good alignmnet transformation is found, it tends to be stable over a flight and, in most cases, over many flights. The transformation may change if the camera undergoes a shock event (such as a hard landing or drop) or if the temperature changes substantially between flights. In these events a new transformation may need to be found.
-Further, since this approach finds a 2-dimensional (affine) transformation between images, it won't work when the parallax between bands results in a 3-dimensional depth field. This can happen if very close to the target or when targets are visible at significantly different ranges, such as a nearby tree or building against a background much farther way. In these cases it will be necessary to use photogrammetry techniques to find a 3-dimensional mapping between images.
-For best alignment results it's good to select a capture which has features which visible in all bands. Man-made objects such as cars, roads, and buildings tend to work very well, while captures of only repeating crop rows tend to work poorly. Remember, once a good transformation has been found for flight, it can be generally be applied across all of the images.
-It's also good to use an image for alignment which is taken near the same level above ground as the rest of the flights. Above approximately 35m AGL, the alignement will be consistent. However, if images taken closer to the ground are used, such as panel images, the same alignment transformation will not work for the flight data.
- -import cv2
-import numpy as np
-import matplotlib.pyplot as plt
-import micasense.imageutils as imageutils
-import micasense.plotutils as plotutils
-
-## Alignment settings
-match_index = 4 # Index of the band, here we use green
-max_alignment_iterations = 20
-warp_mode = cv2.MOTION_HOMOGRAPHY # MOTION_HOMOGRAPHY or MOTION_AFFINE. For Altum images only use HOMOGRAPHY
-pyramid_levels = 3 # for 10-band imagery we use a 3-level pyramid. In some cases
-
-print("Aligning images. Depending on settings this can take from a few seconds to many minutes")
-# Can potentially increase max_iterations for better results, but longer runtimes
-warp_matrices, alignment_pairs = imageutils.align_capture(capture,
- ref_index = match_index,
- max_iterations = max_alignment_iterations,
- warp_mode = warp_mode,
- pyramid_levels = pyramid_levels)
-
-print("Finished Aligning, warp matrices={}".format(warp_matrices))
-Aligning images. Depending on settings this can take from a few seconds to many minutes -Finished aligning band 4 -Finished aligning band 0 -Finished aligning band 3 -Finished aligning band 1 -Finished aligning band 2 -Finished aligning band 6 -Finished aligning band 5 -Finished aligning band 8 -Finished aligning band 7 -Finished aligning band 9 -Finished Aligning, warp matrices=[array([[ 9.7928989e-01, 1.3615261e-04, 1.9574374e-02], - [-3.7023663e-03, 9.9245304e-01, 2.8355631e+01], - [-1.0651237e-05, 2.6095395e-06, 1.0000000e+00]], dtype=float32), array([[ 9.9525303e-01, -7.0169556e-04, -6.3409243e+00], - [-7.1914408e-05, 9.9811482e-01, -6.2239196e-02], - [-3.6460631e-06, 1.8568957e-06, 1.0000000e+00]], dtype=float32), array([[ 9.9596852e-01, -2.5820474e-03, -1.2950540e+01], - [ 2.8751660e-03, 9.9919146e-01, 2.6645420e+01], - [-1.8280600e-06, 2.3902323e-06, 1.0000000e+00]], dtype=float32), array([[ 9.9592507e-01, 4.5777867e-03, 1.4982515e+01], - [-3.4031910e-03, 9.9773425e-01, 2.2611160e+01], - [-9.0995320e-07, 2.3299665e-06, 1.0000000e+00]], dtype=float32), array([[ 1.0000000e+00, 2.8202797e-20, 8.9091027e-15], - [-5.6938062e-19, 1.0000000e+00, 6.4872178e-15], - [-1.2599002e-21, 1.3424955e-23, 1.0000000e+00]], dtype=float32), array([[ 9.9348122e-01, -1.8129036e-02, 1.1003100e+00], - [ 1.2132729e-02, 9.9091125e-01, 2.8742294e+01], - [-2.9285566e-06, -6.2816152e-06, 1.0000000e+00]], dtype=float32), array([[ 9.8954433e-01, -1.8810771e-02, 1.1892141e+01], - [ 1.3862181e-02, 9.8975760e-01, 2.0044521e+01], - [-4.0683108e-06, -4.2664565e-06, 1.0000000e+00]], dtype=float32), array([[ 9.8794365e-01, -1.2153405e-02, 1.8298512e+01], - [ 8.4343022e-03, 9.8829275e-01, 2.8633160e+01], - [-2.5430195e-06, -3.1491891e-06, 1.0000000e+00]], dtype=float32), array([[ 9.9505156e-01, -1.7393423e-02, 3.6898847e+00], - [ 1.3367594e-02, 9.9287426e-01, 1.8703876e+01], - [-1.3000300e-06, -4.7085673e-06, 1.0000000e+00]], dtype=float32), array([[ 9.9182397e-01, -1.2042225e-02, 1.2686370e+01], - [ 7.6934961e-03, 9.8897797e-01, 2.7970253e+01], - [-9.2655773e-07, -6.0620891e-06, 1.0000000e+00]], dtype=float32)] --
After finding image alignments we may need to remove pixels around the edges which aren't present in every image in the capture. To do this we use the affine transforms found above and the image distortions from the image metadata. OpenCV provides a couple of handy helpers for this task in the cv2.undistortPoints() and cv2.transform() methods. These methods takes a set of pixel coordinates and apply our undistortion matrix and our affine transform, respectively. So, just as we did when registering the images, we first apply the undistortion process the coordinates of the image borders, then we apply the affine transformation to that result. The resulting pixel coordinates tell us where the image borders end up after this pair of transformations, and we can then crop the resultant image to these coordinates.
cropped_dimensions, edges = imageutils.find_crop_bounds(capture, warp_matrices, warp_mode=warp_mode)
-im_aligned = imageutils.aligned_capture(capture, warp_matrices, warp_mode, cropped_dimensions, match_index, img_type=img_type)
-display(im_aligned)
-print("fff")
-im_aligned = capture.create_aligned_capture(warp_matrices=warp_matrices, motion_type=warp_mode, img_type=img_type, match_index=match_index)
-display(im_aligned)
-array([[[0.04449825, 0.05834329, 0.05699025, ..., 0.0528319 , - 0.07618934, 0.13635507], - [0.04551493, 0.05871509, 0.0533512 , ..., 0.05375743, - 0.07973219, 0.13832156], - [0.04553824, 0.06390802, 0.05610467, ..., 0.05398859, - 0.08443005, 0.14852835], - ..., - [0.05893659, 0.07230809, 0.05917593, ..., 0.06073643, - 0.09398529, 0.21687469], - [0.05003324, 0.0664675 , 0.05700499, ..., 0.05660931, - 0.08722617, 0.20154284], - [0.04540764, 0.06866231, 0.05971338, ..., 0.05768275, - 0.08611467, 0.19954923]], - - [[0.043996 , 0.05818984, 0.05566213, ..., 0.0515653 , - 0.07534796, 0.12414938], - [0.04554788, 0.05788112, 0.05298578, ..., 0.05173871, - 0.076417 , 0.12747806], - [0.04729695, 0.06531452, 0.05632412, ..., 0.05696491, - 0.07990886, 0.15101017], - ..., - [0.06180993, 0.07528408, 0.06261701, ..., 0.06480689, - 0.1002305 , 0.21519387], - [0.05216246, 0.06675215, 0.05966157, ..., 0.05824959, - 0.09197214, 0.20033701], - [0.04686544, 0.0696869 , 0.06084135, ..., 0.05859195, - 0.08765141, 0.19360156]], - - [[0.04719513, 0.05703212, 0.05228514, ..., 0.05184768, - 0.08116278, 0.130595 ], - [0.04736948, 0.06186526, 0.05192235, ..., 0.05348465, - 0.07854285, 0.13192512], - [0.04892166, 0.06844062, 0.05518164, ..., 0.05929117, - 0.08475667, 0.15392922], - ..., - [0.06107106, 0.07277559, 0.06353831, ..., 0.06612159, - 0.09967496, 0.2089769 ], - [0.05610751, 0.06975867, 0.06273928, ..., 0.0632534 , - 0.09190915, 0.1867479 ], - [0.04991827, 0.07058094, 0.06168615, ..., 0.06201676, - 0.08668306, 0.17329696]], - - ..., - - [[0.03752367, 0.09495717, 0.06308448, ..., 0.06055268, - 0.11859115, 0.36167163], - [0.03716302, 0.09172834, 0.06215791, ..., 0.05819597, - 0.11055551, 0.36911854], - [0.03741182, 0.08998678, 0.05918581, ..., 0.05523164, - 0.10871691, 0.37251878], - ..., - [0.17053236, 0.20084937, 0.2656319 , ..., 0.23127253, - 0.21557146, 0.22828759], - [0.16936235, 0.19811206, 0.2181453 , ..., 0.19929384, - 0.19923317, 0.2175055 ], - [0.17199045, 0.19134107, 0.20697866, ..., 0.20128694, - 0.20539421, 0.21508098]], - - [[0.03794758, 0.09345208, 0.06537206, ..., 0.0648974 , - 0.11819842, 0.36217105], - [0.03913194, 0.09361731, 0.06283107, ..., 0.06258071, - 0.11491562, 0.37606257], - [0.0391482 , 0.09110615, 0.06024874, ..., 0.05788773, - 0.10959219, 0.36674654], - ..., - [0.16886954, 0.20589119, 0.282101 , ..., 0.24447776, - 0.22555993, 0.24048506], - [0.16768672, 0.20429438, 0.2303536 , ..., 0.2075679 , - 0.21025576, 0.22479147], - [0.16860068, 0.1963228 , 0.21333084, ..., 0.20610072, - 0.21020785, 0.2212299 ]], - - [[0.03912217, 0.09133695, 0.06471822, ..., 0.06551831, - 0.11582895, 0.3591986 ], - [0.03989396, 0.09222235, 0.06317808, ..., 0.06262669, - 0.11217536, 0.35820585], - [0.04026981, 0.08995611, 0.0575012 , ..., 0.05748435, - 0.1085691 , 0.35474142], - ..., - [0.17024775, 0.20644474, 0.29202968, ..., 0.2621528 , - 0.23570539, 0.24650206], - [0.16856626, 0.2006478 , 0.24013193, ..., 0.21821336, - 0.21899852, 0.22762969], - [0.16514897, 0.19955555, 0.22367412, ..., 0.21459754, - 0.21694328, 0.22302781]]], dtype=float32)-
fff --
array([[[0.04449825, 0.05834329, 0.05699025, ..., 0.0528319 , - 0.07618934, 0.13635507], - [0.04551493, 0.05871509, 0.0533512 , ..., 0.05375743, - 0.07973219, 0.13832156], - [0.04553824, 0.06390802, 0.05610467, ..., 0.05398859, - 0.08443005, 0.14852835], - ..., - [0.05893659, 0.07230809, 0.05917593, ..., 0.06073643, - 0.09398529, 0.21687469], - [0.05003324, 0.0664675 , 0.05700499, ..., 0.05660931, - 0.08722617, 0.20154284], - [0.04540764, 0.06866231, 0.05971338, ..., 0.05768275, - 0.08611467, 0.19954923]], - - [[0.043996 , 0.05818984, 0.05566213, ..., 0.0515653 , - 0.07534796, 0.12414938], - [0.04554788, 0.05788112, 0.05298578, ..., 0.05173871, - 0.076417 , 0.12747806], - [0.04729695, 0.06531452, 0.05632412, ..., 0.05696491, - 0.07990886, 0.15101017], - ..., - [0.06180993, 0.07528408, 0.06261701, ..., 0.06480689, - 0.1002305 , 0.21519387], - [0.05216246, 0.06675215, 0.05966157, ..., 0.05824959, - 0.09197214, 0.20033701], - [0.04686544, 0.0696869 , 0.06084135, ..., 0.05859195, - 0.08765141, 0.19360156]], - - [[0.04719513, 0.05703212, 0.05228514, ..., 0.05184768, - 0.08116278, 0.130595 ], - [0.04736948, 0.06186526, 0.05192235, ..., 0.05348465, - 0.07854285, 0.13192512], - [0.04892166, 0.06844062, 0.05518164, ..., 0.05929117, - 0.08475667, 0.15392922], - ..., - [0.06107106, 0.07277559, 0.06353831, ..., 0.06612159, - 0.09967496, 0.2089769 ], - [0.05610751, 0.06975867, 0.06273928, ..., 0.0632534 , - 0.09190915, 0.1867479 ], - [0.04991827, 0.07058094, 0.06168615, ..., 0.06201676, - 0.08668306, 0.17329696]], - - ..., - - [[0.03752367, 0.09495717, 0.06308448, ..., 0.06055268, - 0.11859115, 0.36167163], - [0.03716302, 0.09172834, 0.06215791, ..., 0.05819597, - 0.11055551, 0.36911854], - [0.03741182, 0.08998678, 0.05918581, ..., 0.05523164, - 0.10871691, 0.37251878], - ..., - [0.17053236, 0.20084937, 0.2656319 , ..., 0.23127253, - 0.21557146, 0.22828759], - [0.16936235, 0.19811206, 0.2181453 , ..., 0.19929384, - 0.19923317, 0.2175055 ], - [0.17199045, 0.19134107, 0.20697866, ..., 0.20128694, - 0.20539421, 0.21508098]], - - [[0.03794758, 0.09345208, 0.06537206, ..., 0.0648974 , - 0.11819842, 0.36217105], - [0.03913194, 0.09361731, 0.06283107, ..., 0.06258071, - 0.11491562, 0.37606257], - [0.0391482 , 0.09110615, 0.06024874, ..., 0.05788773, - 0.10959219, 0.36674654], - ..., - [0.16886954, 0.20589119, 0.282101 , ..., 0.24447776, - 0.22555993, 0.24048506], - [0.16768672, 0.20429438, 0.2303536 , ..., 0.2075679 , - 0.21025576, 0.22479147], - [0.16860068, 0.1963228 , 0.21333084, ..., 0.20610072, - 0.21020785, 0.2212299 ]], - - [[0.03912217, 0.09133695, 0.06471822, ..., 0.06551831, - 0.11582895, 0.3591986 ], - [0.03989396, 0.09222235, 0.06317808, ..., 0.06262669, - 0.11217536, 0.35820585], - [0.04026981, 0.08995611, 0.0575012 , ..., 0.05748435, - 0.1085691 , 0.35474142], - ..., - [0.17024775, 0.20644474, 0.29202968, ..., 0.2621528 , - 0.23570539, 0.24650206], - [0.16856626, 0.2006478 , 0.24013193, ..., 0.21821336, - 0.21899852, 0.22762969], - [0.16514897, 0.19955555, 0.22367412, ..., 0.21459754, - 0.21694328, 0.22302781]]], dtype=float32)-
Once the transformation has been found, it can be verified by composting the aligned images to check alignment. The image 'stack' containing all bands can also be exported to a multi-band TIFF file for viewing in extrernal software such as QGIS. Usef ul componsites are a naturally colored RGB as well as color infrared, or CIR.
- -# figsize=(30,23) # use this size for full-image-resolution display
-figsize=(16,13) # use this size for export-sized display
-
-rgb_band_indices = [capture.band_names().index('Red'),capture.band_names().index('Green'),capture.band_names().index('Blue-444')]
-cir_band_indices = [capture.band_names().index('NIR'),capture.band_names().index('Red'),capture.band_names().index('Green')]
-
-# Create a normalized stack for viewing
-im_display = np.zeros((im_aligned.shape[0],im_aligned.shape[1],im_aligned.shape[2]), dtype=np.float32 )
-
-im_min = np.percentile(im_aligned[:,:,rgb_band_indices].flatten(), 0.5) # modify these percentiles to adjust contrast
-im_max = np.percentile(im_aligned[:,:,rgb_band_indices].flatten(), 99.5) # for many images, 0.5 and 99.5 are good values
-
-# for rgb true color, we use the same min and max scaling across the 3 bands to
-# maintain the "white balance" of the calibrated image
-for i in rgb_band_indices:
- im_display[:,:,i] = imageutils.normalize(im_aligned[:,:,i], im_min, im_max)
-
-rgb = im_display[:,:,rgb_band_indices]
-
-# for cir false color imagery, we normalize the NIR,R,G bands within themselves, which provides
-# the classical CIR rendering where plants are red and soil takes on a blue tint
-for i in cir_band_indices:
- im_display[:,:,i] = imageutils.normalize(im_aligned[:,:,i])
-
-cir = im_display[:,:,cir_band_indices]
-fig, axes = plt.subplots(1, 2, figsize=figsize)
-axes[0].set_title("Red-Green-Blue Composite")
-axes[0].imshow(rgb)
-axes[1].set_title("Color Infrared (CIR) Composite")
-axes[1].imshow(cir)
-plt.show()
-There are many techniques for image enhancement, but one which is commonly used to improve the visual sharpness of imagery is the unsharp mask. Here we apply an unsharp mask to the RGB image to improve the visualization, and then apply a gamma curve to make the darkest areas brighter.
- -# Create an enhanced version of the RGB render using an unsharp mask
-gaussian_rgb = cv2.GaussianBlur(rgb, (9,9), 10.0)
-gaussian_rgb[gaussian_rgb<0] = 0
-gaussian_rgb[gaussian_rgb>1] = 1
-unsharp_rgb = cv2.addWeighted(rgb, 1.5, gaussian_rgb, -0.5, 0)
-unsharp_rgb[unsharp_rgb<0] = 0
-unsharp_rgb[unsharp_rgb>1] = 1
-
-# Apply a gamma correction to make the render appear closer to what our eyes would see
-gamma = 1.4
-gamma_corr_rgb = unsharp_rgb**(1.0/gamma)
-fig = plt.figure(figsize=figsize)
-plt.imshow(gamma_corr_rgb, aspect='equal')
-plt.axis('off')
-plt.show()
-Composite images can be exported to JPEG or PNG format using the imageio package. These images may be useful for visualization or thumbnailing, and creating RGB thumbnails of a set of images can provide a convenient way to browse the imagery in a more visually appealing way that browsing the raw imagery.
import imageio
-imtype = 'png' # or 'jpg'
-imageio.imwrite('rgb.'+imtype, (255*gamma_corr_rgb).astype('uint8'))
-imageio.imwrite('cir.'+imtype, (255*cir).astype('uint8'))
-We can export the image easily stacks using the gdal library (http://www.glal.org). Once exported, these image stacks can be opened in software such as QGIS and raster operations such as NDVI or NDRE computation can be done in that software. At this time the stacks don't include any geographic information.
from osgeo import gdal, gdal_array
-rows, cols, bands = im_display.shape
-driver = gdal.GetDriverByName('GTiff')
-filename = "bbggrreeen" #blue,blue,green,green,red,red,rededge,rededge,rededge,nir
-
-sort_by_wavelength = True # set to false if you want stacks in camera-band-index order
-
-outRaster = driver.Create(filename+"C.tiff", cols, rows, bands, gdal.GDT_UInt16, options = [ 'INTERLEAVE=BAND','COMPRESS=DEFLATE' ])
-try:
- if outRaster is None:
- raise IOError("could not load gdal GeoTiff driver")
-
- if sort_by_wavelength:
- eo_list = list(np.argsort(capture.center_wavelengths()))
- else:
- eo_list = capture.eo_indices()
-
- for outband,inband in enumerate(eo_list):
- outband = outRaster.GetRasterBand(outband+1)
- outdata = im_aligned[:,:,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(capture.lw_indices()):
- outband = outRaster.GetRasterBand(len(eo_list)+outband+1)
- outdata = (im_aligned[:,:,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
-"Stacks" as described above are useful in a number of processing cases. For example, at the time of this writing, many photogrammetry suites could import and process stack files without significantly impacting the radiometric processing which has already been accomplished.
-Running photogrammetry on stack files instead of raw image files has both advantages and drawbacks. The primary advantage has been found to be an increase in processing speed and a reduction in program memory usage. As the photogrammetric workflow generally operates on luminance images and may not use color information, stacked images may require similar resources and be processed at a similar speed as single-band images. This is because one band of the stack can be used to generate the matching feature space while the others are ignored for matching purposes. This reduces the feature space 5-fold over matching using all images separately.
-One disadvantage is that stacking images outside of the photogrammetric workflow may result in poor image matching. The RedEdge is known to have stable lens characteristics over the course of normal operation, but variations in temperature or impacts to the camera through handling or rough landings may change the image alignment parameters. For this reason, we recommend finding a matching transformation for each flight (each take-off and landing). Alignment transformations from multiple images within a flight can be compared to find the best transformation to apply to the set of the flight. While not described or supported in this generic implementation, some matching algorithms can use a "seed" value as a starting point to speed up matching. For most cases, this seed could be the transformation found in a previous flight, or another source of a known good transformation.
- -For raw index computation on single images, the numpy package provides a simple way to do math and simple visualizatoin on images. Below, we compute and visualize an image histogram and then use that to pick a colormap range for visualizing the NDVI of an image.
-After computing the NDVI and prior to displaying it, we use a very rudimentary method for focusing on the plants and removing the soil and shadow information from our images and histograms. Below we remove non-plant pixels by setting to zero any pixels in the image where the NIR reflectance is less than 20%. This helps to ensure that the NDVI and NDRE histograms aren't skewed substantially by soil noise.
- -from micasense import plotutils
-import matplotlib.pyplot as plt
-
-nir_band = [name.lower() for name in capture.band_names()].index('nir')
-red_band = [name.lower() for name in capture.band_names()].index('red')
-
-np.seterr(divide='ignore', invalid='ignore') # ignore divide by zero errors in the index calculation
-
-# Compute Normalized Difference Vegetation Index (NDVI) from the NIR(3) and RED (2) bands
-ndvi = (im_aligned[:,:,nir_band] - im_aligned[:,:,red_band]) / (im_aligned[:,:,nir_band] + im_aligned[:,:,red_band])
-
-# remove shadowed areas (mask pixels with NIR reflectance < 20%))
-if img_type == 'reflectance':
- ndvi = np.ma.masked_where(im_aligned[:,:,nir_band] < 0.20, ndvi)
-elif img_type == 'radiance':
- lower_pct_radiance = np.percentile(im_aligned[:,:,3], 10.0)
- ndvi = np.ma.masked_where(im_aligned[:,:,nir_band] < lower_pct_radiance, ndvi)
-
-# Compute and display a histogram
-ndvi_hist_min = np.min(ndvi)
-ndvi_hist_max = np.max(ndvi)
-fig, axis = plt.subplots(1, 1, figsize=(10,4))
-axis.hist(ndvi.ravel(), bins=512, range=(ndvi_hist_min, ndvi_hist_max))
-plt.title("NDVI Histogram")
-plt.show()
-
-min_display_ndvi = 0.45 # further mask soil by removing low-ndvi values
-#min_display_ndvi = np.percentile(ndvi.flatten(), 5.0) # modify with these percentilse to adjust contrast
-max_display_ndvi = np.percentile(ndvi.flatten(), 99.5) # for many images, 0.5 and 99.5 are good values
-masked_ndvi = np.ma.masked_where(ndvi < min_display_ndvi, ndvi)
-
-#reduce the figure size to account for colorbar
-figsize=np.asarray(figsize) - np.array([3,2])
-
-#plot NDVI over an RGB basemap, with a colorbar showing the NDVI scale
-fig, axis = plotutils.plot_overlay_withcolorbar(gamma_corr_rgb,
- masked_ndvi,
- figsize = figsize,
- title = 'NDVI filtered to only plants over RGB base layer',
- vmin = min_display_ndvi,
- vmax = max_display_ndvi)
-fig.savefig('ndvi_over_rgb.png')
-In the same manner, we can compute, filter, and display another index useful for the RedEdge camera, the Normalized Difference Red Edge (NDRE) index. We also filter out shadows and soil to ensure our display focuses only on the plant health.
- -# Compute Normalized Difference Red Edge Index from the NIR(3) and RedEdge(4) bands
-rededge_band = [name.lower() for name in capture.band_names()].index('red edge')
-ndre = (im_aligned[:,:,nir_band] - im_aligned[:,:,rededge_band]) / (im_aligned[:,:,nir_band] + im_aligned[:,:,rededge_band])
-
-# Mask areas with shadows and low NDVI to remove soil
-masked_ndre = np.ma.masked_where(ndvi < min_display_ndvi, ndre)
-
-# Compute a histogram
-ndre_hist_min = np.min(masked_ndre)
-ndre_hist_max = np.max(masked_ndre)
-fig, axis = plt.subplots(1, 1, figsize=(10,4))
-axis.hist(masked_ndre.ravel(), bins=512, range=(ndre_hist_min, ndre_hist_max))
-plt.title("NDRE Histogram (filtered to only plants)")
-plt.show()
-
-min_display_ndre = np.percentile(masked_ndre, 5)
-max_display_ndre = np.percentile(masked_ndre, 99.5)
-
-fig, axis = plotutils.plot_overlay_withcolorbar(gamma_corr_rgb,
- masked_ndre,
- figsize=figsize,
- title='NDRE filtered to only plants over RGB base layer',
- vmin=min_display_ndre,vmax=max_display_ndre)
-fig.savefig('ndre_over_rgb.png')
-Finally, we show a classic agricultural remote sensing output in the tassled cap plot. This plot can be useful for visualizing row crops and plots the Red Reflectance channel on the X-axis against the NIR reflectance channel on the Y-axis. This plot also clearly shows the line of the soil in that space. The tassled cap view isn't very useful for this arid data set; however, we can see the "badge of trees" of high NIR reflectance and relatively low red reflectance. This provides an example of one of the uses of aligned images for single capture analysis.
- -x_band = red_band
-y_band = nir_band
-x_max = np.max(im_aligned[:,:,x_band])
-y_max = np.max(im_aligned[:,:,y_band])
-
-fig = plt.figure(figsize=(12,12))
-plt.hexbin(im_aligned[:,:,x_band],im_aligned[:,:,y_band],gridsize=640,bins='log',extent=(0,x_max,0,y_max))
-ax = fig.gca()
-ax.set_xlim([0,x_max])
-ax.set_ylim([0,y_max])
-plt.xlabel("{} Reflectance".format(capture.band_names()[x_band]))
-plt.ylabel("{} Reflectance".format(capture.band_names()[y_band]))
-plt.show()
-Last, we output the warp_matrices that we got for this image stack for usage elsewhere. Currently these can be used in the Batch Processing.ipynb notebook to save reflectance-compensated stacks of images to a directory.
print(warp_matrices)
-[array([[ 9.7928989e-01, 1.3615261e-04, 1.9574374e-02], - [-3.7023663e-03, 9.9245304e-01, 2.8355631e+01], - [-1.0651237e-05, 2.6095395e-06, 1.0000000e+00]], dtype=float32), array([[ 9.9525303e-01, -7.0169556e-04, -6.3409243e+00], - [-7.1914408e-05, 9.9811482e-01, -6.2239196e-02], - [-3.6460631e-06, 1.8568957e-06, 1.0000000e+00]], dtype=float32), array([[ 9.9596852e-01, -2.5820474e-03, -1.2950540e+01], - [ 2.8751660e-03, 9.9919146e-01, 2.6645420e+01], - [-1.8280600e-06, 2.3902323e-06, 1.0000000e+00]], dtype=float32), array([[ 9.9592507e-01, 4.5777867e-03, 1.4982515e+01], - [-3.4031910e-03, 9.9773425e-01, 2.2611160e+01], - [-9.0995320e-07, 2.3299665e-06, 1.0000000e+00]], dtype=float32), array([[ 1.0000000e+00, 2.8202797e-20, 8.9091027e-15], - [-5.6938062e-19, 1.0000000e+00, 6.4872178e-15], - [-1.2599002e-21, 1.3424955e-23, 1.0000000e+00]], dtype=float32), array([[ 9.9348122e-01, -1.8129036e-02, 1.1003100e+00], - [ 1.2132729e-02, 9.9091125e-01, 2.8742294e+01], - [-2.9285566e-06, -6.2816152e-06, 1.0000000e+00]], dtype=float32), array([[ 9.8954433e-01, -1.8810771e-02, 1.1892141e+01], - [ 1.3862181e-02, 9.8975760e-01, 2.0044521e+01], - [-4.0683108e-06, -4.2664565e-06, 1.0000000e+00]], dtype=float32), array([[ 9.8794365e-01, -1.2153405e-02, 1.8298512e+01], - [ 8.4343022e-03, 9.8829275e-01, 2.8633160e+01], - [-2.5430195e-06, -3.1491891e-06, 1.0000000e+00]], dtype=float32), array([[ 9.9505156e-01, -1.7393423e-02, 3.6898847e+00], - [ 1.3367594e-02, 9.9287426e-01, 1.8703876e+01], - [-1.3000300e-06, -4.7085673e-06, 1.0000000e+00]], dtype=float32), array([[ 9.9182397e-01, -1.2042225e-02, 1.2686370e+01], - [ 7.6934961e-03, 9.8897797e-01, 2.7970253e+01], - [-9.2655773e-07, -6.0620891e-06, 1.0000000e+00]], dtype=float32)] --
Copyright (c) 2017-2018 MicaSense, Inc. For licensing information see the project git repository
- -This workbook shows a completely passive alignment function using only the RigRelatives present in image metadata. Older versions of firmware may not include these tags. If your images don't, check the MicaSense support site for some tips on how to update your camera firmware to have them, as well as how to add them to datasets prior to the update.
-While this method doesn't provide perfect alignment, it can be fast and very useful for visualization of images when processing power is limited or speed is more important than alignment quality.
-As we have done in previous examples, we use the micasense.capture class to open, radiometrically correct, and visualize the 5 bands of a RedEdge capture.
-First, we'll load the autoreload extension. This lets us change underlying code (such as library functions) without having to reload the entire workbook and kernel. This is useful in this workbook because the cell that runs the alignment can take a long time to run, so with autoreload extension we can change external code for analysis and visualization without needing to re-compute the alignments each time.
%load_ext autoreload
-%autoreload 2
-import os, glob
-import micasense.capture as capture
-%matplotlib inline
-from pathlib import Path
-
-panelNames = None
-
-# This is an Altum image with RigRelatives and a thermal band
-
-imagePath = Path("./data/ALTUM")
-
-# these will return lists of image paths as strings
-imageNames = list(imagePath.glob('IMG_0021_*.tif'))
-imageNames = [x.as_posix() for x in imageNames]
-
-panelNames = list(imagePath.glob('IMG_0000_*.tif'))
-panelNames = [x.as_posix() for x in panelNames]
-
-if panelNames is not None:
- panelCap = capture.Capture.from_filelist(panelNames)
-else:
- panelCap = None
-
-capture = capture.Capture.from_filelist(imageNames)
-
-for img in capture.images:
- if img.rig_relatives is None:
- raise ValueError("Images must have RigRelatives tags set which requires updated firmware and calibration. See the links in text above")
-
-if panelCap is not None:
- if panelCap.panel_albedo() is not None:
- panel_reflectance_by_band = panelCap.panel_albedo()
- else:
- panel_reflectance_by_band = [0.49, 0.49, 0.49, 0.49, 0.49] #RedEdge band_index order
- panel_irradiance = panelCap.panel_irradiance(panel_reflectance_by_band)
- img_type = "reflectance"
- capture.plot_undistorted_reflectance(panel_irradiance)
-else:
- if False: #capture.dls_present():
- img_type='reflectance'
- capture.plot_undistorted_reflectance(capture.dls_irradiance())
- else:
- img_type = "radiance"
- capture.plot_undistorted_radiance()
-For images with RigRelative tags present, we can find a rough alignment using only the built in relatives. These can be good for quick visualizatoins. For better results, use an algorithm like that in the other image alignment tutorial to sweeten the alignment.
- -import cv2
-import numpy as np
-import matplotlib.pyplot as plt
-import micasense.imageutils as imageutils
-import micasense.plotutils as plotutils
-
-reference_band = 5
-warp_mode = cv2.MOTION_HOMOGRAPHY
-warp_matrices = capture.get_warp_matrices(ref_index=reference_band)
-
-cropped_dimensions,edges = imageutils.find_crop_bounds(capture,warp_matrices,reference_band=reference_band)
-im_aligned = imageutils.aligned_capture(capture, warp_matrices, warp_mode, cropped_dimensions, reference_band, img_type=img_type)
-
-print(im_aligned.shape)
-
-print("warp_matrices={}".format(warp_matrices))
-(91, 123, 6) -warp_matrices=[array([[ 1.57521370e+01, 2.64721217e-01, -2.44963917e+02], - [-2.94310824e-01, 1.58442787e+01, -1.34414613e+02], - [-2.77775663e-05, -9.67369561e-06, 1.00000000e+00]]), array([[ 1.55968786e+01, 2.34242421e-01, -2.26810378e+02], - [-2.78826529e-01, 1.56989040e+01, -1.13524137e+02], - [-4.23356370e-05, -1.27338303e-05, 1.00000000e+00]]), array([[ 1.56259258e+01, 3.15964099e-01, -1.90938966e+02], - [-3.41730559e-01, 1.57251597e+01, -9.97306451e+01], - [-2.77625565e-05, -5.27840773e-06, 1.00000000e+00]]), array([[ 1.56569908e+01, 2.40978299e-01, -2.27165074e+02], - [-2.94899115e-01, 1.57570534e+01, -1.33111975e+02], - [-4.09540877e-05, -2.33779890e-05, 1.00000000e+00]]), array([[ 1.55950891e+01, 2.07864879e-01, -2.33521065e+02], - [-2.85249968e-01, 1.56841110e+01, -1.31908090e+02], - [-5.30525636e-05, -3.78432070e-05, 1.00000000e+00]]), array([[1.00000000e+00, 1.02405070e-18, 0.00000000e+00], - [1.09972007e-18, 1.00000000e+00, 0.00000000e+00], - [1.60052412e-21, 8.52545510e-23, 1.00000000e+00]])] --
Once the transformation has been found, it can be verified by composting the aligned images to check alignment. The image 'stack' containing all bands can also be exported to a multi-band TIFF file for viewing in extrernal software such as QGIS. Usef ul componsites are a naturally colored RGB as well as color infrared, or CIR.
- -# figsize=(30,23) # use this size for full-image-resolution display
-figsize=(16,13) # use this size for export-sized display
-
-rgb_band_indices = [2,1,0]
-cir_band_indices = [3,2,1]
-
-# Create an empty normalized stack for viewing
-im_display = np.zeros((im_aligned.shape[0],im_aligned.shape[1],im_aligned.shape[2]), dtype=np.float32 )
-
-im_min = np.percentile(im_aligned[:,:,0:2].flatten(), 0.1) # modify with these percentilse to adjust contrast
-im_max = np.percentile(im_aligned[:,:,0:2].flatten(), 99.9) # for many images, 0.5 and 99.5 are good values
-
-for i in range(0,im_aligned.shape[2]):
- if img_type == 'reflectance':
- # for reflectance images we maintain white-balance by applying the same display scaling to all bands
- im_display[:,:,i] = imageutils.normalize(im_aligned[:,:,i], im_min, im_max)
- elif img_type == 'radiance':
- # for radiance images we do an auto white balance since we don't know the input light spectrum by
- # stretching each display band histogram to it's own min and max
- im_display[:,:,i] = imageutils.normalize(im_aligned[:,:,i])
-
-rgb = im_display[:,:,rgb_band_indices]
-# for cir false color imagery, we normalize the NIR,R,G bands within themselves, which provides
-# the classical CIR rendering where plants are red and soil takes on a blue tint
-for i in cir_band_indices:
- im_display[:,:,i] = imageutils.normalize(im_aligned[:,:,i])
-
-cir = im_display[:,:,cir_band_indices]
-fig, axes = plt.subplots(1, 2, figsize=figsize)
-axes[0].set_title("Red-Green-Blue Composite")
-axes[0].imshow(rgb)
-axes[1].set_title("Color Infrared (CIR) Composite")
-axes[1].imshow(cir)
-plt.show()
-There are many techniques for image enhancement, but one which is commonly used to improve the visual sharpness of imagery is the unsharp mask. Here we apply an unsharp mask to the RGB image to improve the visualization, and then apply a gamma curve to make the darkest areas brighter.
- -# Create an enhanced version of the RGB render using an unsharp mask
-gaussian_rgb = cv2.GaussianBlur(rgb, (9,9), 10.0)
-gaussian_rgb[gaussian_rgb<0] = 0
-gaussian_rgb[gaussian_rgb>1] = 1
-unsharp_rgb = cv2.addWeighted(rgb, 1.5, gaussian_rgb, -0.5, 0)
-unsharp_rgb[unsharp_rgb<0] = 0
-unsharp_rgb[unsharp_rgb>1] = 1
-
-# Apply a gamma correction to make the render appear closer to what our eyes would see
-gamma = 1.4
-gamma_corr_rgb = unsharp_rgb**(1.0/gamma)
-fig = plt.figure(figsize=figsize)
-plt.imshow(gamma_corr_rgb, aspect='equal')
-plt.axis('off')
-plt.show()
-We can output the image to a PNG or JPEG file for viewing. This can also be useful in creating thumbnails of captures.
- -import imageio
-imtype = 'png' # or 'jpg'
-imageio.imwrite('rgb.'+imtype, (255*gamma_corr_rgb).astype('uint8'))
-We can export the image easily stacks using the gdal library (http://www.glal.org). Once exported, these image stacks can be opened in software such as QGIS and raster operations such as NDVI or NDRE computation can be done in that software. At this time the stacks don't include any geographic information.
from osgeo import gdal, gdal_array
-rows, cols, bands = im_display.shape
-driver = gdal.GetDriverByName('GTiff')
-filename = "bgrne" #blue,green,red,nir,redEdge
-
-if im_aligned.shape[2] == 6:
- filename = filename + "t" #thermal
-outRaster = driver.Create(filename+".tiff", cols, rows, im_aligned.shape[2], gdal.GDT_UInt16)
-
-normalize = (img_type == 'radiance') # normalize radiance images to fit with in UInt16
-
-# Output a 'stack' in the same band order as RedEdge/Alutm
-# Blue,Green,Red,NIR,RedEdge[,Thermal]
-# reflectance stacks are output with 32768=100% reflectance to provide some overhead for specular reflections
-# radiance stacks are output with 65535=100% radiance to provide some overhead for specular reflections
-
-# NOTE: NIR and RedEdge are not in wavelength order!
-
-multispec_min = np.min(im_aligned[:,:,1:5])
-multispec_max = np.max(im_aligned[:,:,1:5])
-
-for i in range(0,5):
- outband = outRaster.GetRasterBand(i+1)
- if normalize:
- outdata = imageutils.normalize(im_aligned[:,:,i],multispec_min,multispec_max)
- else:
- outdata = im_aligned[:,:,i]
- outdata[outdata<0] = 0
- outdata[outdata>2] = 2
-
- outdata = outdata*32767
- outdata[outdata<0] = 0
- outdata[outdata>65535] = 65535
- outband.WriteArray(outdata)
- outband.FlushCache()
-
-if im_aligned.shape[2] == 6:
- outband = outRaster.GetRasterBand(6)
- outdata = im_aligned[:,:,5] * 100 # scale to centi-C to fit into uint16
- outdata[outdata<0] = 0
- outdata[outdata>65535] = 65535
- outband.WriteArray(outdata)
- outband.FlushCache()
-outRaster = None
-"Stacks" as described above are useful in a number of processing cases. For example, at the time of this writing, many photogrammetry suites could import and process stack files without significantly impacting the radiometric processing which has already been accomplished.
-Running photogrammetry on stack files instead of raw image files has both advantages and drawbacks. The primary advantage has been found to be an increase in processing speed and a reduction in program memory usage. As the photogrammetric workflow generally operates on luminance images and may not use color information, stacked images may require similar resources and be processed at a similar speed as single-band images. This is because one band of the stack can be used to generate the matching feature space while the others are ignored for matching purposes. This reduces the feature space 5-fold over matching using all images separately.
-One disadvantage is that stacking images outside of the photogrammetric workflow may result in poor image matching. The RedEdge is known to have stable lens characteristics over the course of normal operation, but variations in temperature or impacts to the camera through handling or rough landings may change the image alignment parameters. For this reason, we recommend finding a matching transformation for each flight (each take-off and landing). Alignment transformations from multiple images within a flight can be compared to find the best transformation to apply to the set of the flight. While not described or supported in this generic implementation, some matching algorithms can use a "seed" value as a starting point to speed up matching. For most cases, this seed could be the transformation found in a previous flight, or another source of a known good transformation.
- -For raw index computation on single images, the numpy package provides a simple way to do math and simple visualizatoin on images. Below, we compute and visualize an image histogram and then use that to pick a colormap range for visualizing the NDVI of an image.
-After computing the NDVI and prior to displaying it, we use a very rudimentary method for focusing on the plants and removing the soil and shadow information from our images and histograms. Below we remove non-plant pixels by setting to zero any pixels in the image where the NIR reflectance is less than 20%. This helps to ensure that the NDVI and NDRE histograms aren't skewed substantially by soil noise.
- -from micasense import plotutils
-import matplotlib.pyplot as plt
-
-np.seterr(divide='ignore', invalid='ignore') # ignore divide by zero errors in the index calculation
-
-# Compute Normalized Difference Vegetation Index (NDVI) from the NIR(3) and RED (2) bands
-ndvi = (im_aligned[:,:,3] - im_aligned[:,:,2]) / (im_aligned[:,:,3] + im_aligned[:,:,2])
-
-# remove shadowed areas (mask pixels with NIR reflectance < 20%))
-if img_type == 'reflectance':
- ndvi = np.ma.masked_where(im_aligned[:,:,3] < 0.20, ndvi)
-elif img_type == 'radiance':
- lower_pct_radiance = np.percentile(im_aligned[:,:,3], 10.0)
- ndvi = np.ma.masked_where(im_aligned[:,:,3] < lower_pct_radiance, ndvi)
-
-# Compute and display a histogram
-ndvi_hist_min = np.min(ndvi)
-ndvi_hist_max = np.max(ndvi)
-fig, axis = plt.subplots(1, 1, figsize=(10,4))
-axis.hist(ndvi.ravel(), bins=512, range=(ndvi_hist_min, ndvi_hist_max))
-plt.title("NDVI Histogram")
-plt.show()
-
-min_display_ndvi = 0.45 # further mask soil by removing low-ndvi values
-#min_display_ndvi = np.percentile(ndvi.flatten(), 5.0) # modify with these percentilse to adjust contrast
-max_display_ndvi = np.percentile(ndvi.flatten(), 99.5) # for many images, 0.5 and 99.5 are good values
-masked_ndvi = np.ma.masked_where(ndvi < min_display_ndvi, ndvi)
-
-#reduce the figure size to account for colorbar
-figsize=np.asarray(figsize) - np.array([3,2])
-
-#plot NDVI over an RGB basemap, with a colorbar showing the NDVI scale
-fig, axis = plotutils.plot_overlay_withcolorbar(gamma_corr_rgb,
- masked_ndvi,
- figsize = figsize,
- title = 'NDVI filtered to only plants over RGB base layer',
- vmin = min_display_ndvi,
- vmax = max_display_ndvi)
-fig.savefig('ndvi_over_rgb.png')
-In the same manner, we can compute, filter, and display another index useful for the RedEdge camera, the Normalized Difference Red Edge (NDRE) index. We also filter out shadows and soil to ensure our display focuses only on the plant health.
- -# Compute Normalized Difference Red Edge Index from the NIR(3) and RedEdge(4) bands
-ndre = (im_aligned[:,:,3] - im_aligned[:,:,4]) / (im_aligned[:,:,3] + im_aligned[:,:,4])
-
-# Mask areas with shadows and low NDVI to remove soil
-masked_ndre = np.ma.masked_where(ndvi < min_display_ndvi, ndre)
-
-# Compute a histogram
-ndre_hist_min = np.min(masked_ndre)
-ndre_hist_max = np.max(masked_ndre)
-fig, axis = plt.subplots(1, 1, figsize=(10,4))
-axis.hist(masked_ndre.ravel(), bins=512, range=(ndre_hist_min, ndre_hist_max))
-plt.title("NDRE Histogram (filtered to only plants)")
-plt.show()
-
-min_display_ndre = np.percentile(masked_ndre, 5)
-max_display_ndre = np.percentile(masked_ndre, 99.5)
-
-fig, axis = plotutils.plot_overlay_withcolorbar(gamma_corr_rgb,
- masked_ndre,
- figsize=figsize,
- title='NDRE filtered to only plants over RGB base layer',
- vmin=min_display_ndre,vmax=max_display_ndre)
-fig.savefig('ndre_over_rgb.png')
-If our image is from an Altum and includes a thermal band, we can display the re-sampled and aligned thermal data over the RGB data to maintain the context of the thermal information.
-In the image below, it's very clear based on the average temperature where the soil is wet and dry, and even in the middle of the road we can find some wetter areas of soil.
- -if im_aligned.shape[2] >= 5:
-
- # by default we don't mask the thermal, since it's native resolution is much lower than the MS
- masked_thermal = im_aligned[:,:,5]
- # Alternatively we can mask the thermal only to plants here, which is useful for large contiguous areas
- # masked_thermal = np.ma.masked_where(ndvi < 0.45, im_aligned[:,:,5])
-
-
- # Compute a histogram
- fig, axis = plt.subplots(1, 1, figsize=(10,4))
- axis.hist(masked_thermal.ravel(), bins=512, range=(np.min(masked_thermal), np.max(masked_thermal)))
- plt.title("Thermal Histogram")
- plt.show()
-
- min_display_therm = np.percentile(masked_thermal, 1)
- max_display_therm = np.percentile(masked_thermal, 99)
-
- fig, axis = plotutils.plot_overlay_withcolorbar(gamma_corr_rgb,
- masked_thermal,
- figsize=figsize,
- title='Temperature over True Color',
- vmin=min_display_therm,vmax=max_display_therm,
- overlay_alpha=0.25,
- overlay_colormap='jet',
- overlay_steps=16,
- display_contours=True,
- contour_steps=16,
- contour_alpha=.4,
- contour_fmt="%.0fC")
- fig.savefig('thermal_over_rgb.png')
-Copyright (c) 2017-2019 MicaSense, Inc. For licensing information see the project git repository
- -In most use cases, each band of a multispectral capture must be aligned with the other bands in order to create meaningful data. In this tutorial, we show how to align the band to each other using open source OpenCV utilities.
-Image alignment allows the combination of images into true-color (RGB) and false color (such as CIR) composites, useful for scouting using single images as well as for display and management uses. In addition to composite images, alignment allows the calculation of pixel-accurate indices such as NDVI or NDRE at the single image level which can be very useful for applications like plant counting and coverage estimations, where mosaicing artifacts may otherwise skew analysis results.
-The image alignment method described below tends to work well on images with abundant image features, or areas of significant contrast. Cars, buildings, parking lots, and roads tend to provide the best results. This approach may not work well on images which contain few features or very repetitive features, such as full canopy row crops or fields of repetitive small crops such lettuce or strawberries. We will disscuss more about the advantages and disadvantages of these methods below.
-The functions behind this alignment process can work with most versions of RedEdge and Altum firmware. They will work best with versions above 3.2.0 which include the "RigRelatives" tags. These tags provide a starting point for the image transformation and can help to ensure convergence of the algorithm.
-As we have done in previous examples, we use the micasense.capture class to open, radiometrically correct, and visualize the 5 bands of a RedEdge capture.
-First, we'll load the autoreload extension. This lets us change underlying code (such as library functions) without having to reload the entire workbook and kernel. This is useful in this workbook because the cell that runs the alignment can take a long time to run, so with autoreload extension we can update the code after the alignment step for analysis and visualization without needing to re-compute the alignments each time.
%load_ext autoreload
-%autoreload 2
-import os, glob
-import micasense.capture as capture
-%matplotlib inline
-from pathlib import Path
-import matplotlib.pyplot as plt
-plt.rcParams["figure.facecolor"] = "w"
-
-panelNames = None
-
-imagePath = Path("./data/REDEDGE-MX")
-
-# these will return lists of image paths as strings
-imageNames = list(imagePath.glob('IMG_0020_*.tif'))
-imageNames = [x.as_posix() for x in imageNames]
-
-panelNames = list(imagePath.glob('IMG_0001_*.tif'))
-panelNames = [x.as_posix() for x in panelNames]
-
-# Allow this code to align both radiance and reflectance images; bu excluding
-# a definition for panelNames above, radiance images will be used
-# For panel images, efforts will be made to automatically extract the panel information
-# but if the panel/firmware is before Altum 1.3.5, RedEdge 5.1.7 the panel reflectance
-# will need to be set in the panel_reflectance_by_band variable.
-# Note: radiance images will not be used to properly create NDVI/NDRE images below.
-if panelNames is not None:
- panelCap = capture.Capture.from_filelist(panelNames)
-else:
- panelCap = None
-
-capture = capture.Capture.from_filelist(imageNames)
-
-if panelCap is not None:
- if panelCap.panel_albedo() is not None:
- panel_reflectance_by_band = panelCap.panel_albedo()
- else:
- panel_reflectance_by_band = [0.49, 0.49, 0.49, 0.49, 0.49] #RedEdge band_index order
- panel_irradiance = panelCap.panel_irradiance(panel_reflectance_by_band)
- img_type = "reflectance"
- capture.plot_undistorted_reflectance(panel_irradiance)
-else:
- if capture.dls_present():
- img_type='reflectance'
- capture.plot_undistorted_reflectance(capture.dls_irradiance())
- else:
- img_type = "radiance"
- capture.plot_undistorted_radiance()
-Alignment is a three step process:
-We provide utilities to find the alignement transformations within a single capture. Our experience shows that once a good alignmnet transformation is found, it tends to be stable over a flight and, in most cases, over many flights. The transformation may change if the camera undergoes a shock event (such as a hard landing or drop) or if the temperature changes substantially between flights. In these events a new transformation may need to be found.
-Further, since this approach finds a 2-dimensional (affine) transformation between images, it won't work when the parallax between bands results in a 3-dimensional depth field. This can happen if very close to the target or when targets are visible at significantly different ranges, such as a nearby tree or building against a background much farther way. In these cases it will be necessary to use photogrammetry techniques to find a 3-dimensional mapping between images.
-For best alignment results it's good to select a capture which has features which visible in all bands. Man-made objects such as cars, roads, and buildings tend to work very well, while captures of only repeating crop rows tend to work poorly. Remember, once a good transformation has been found for flight, it can be generally be applied across all of the images.
-It's also good to use an image for alignment which is taken near the same level above ground as the rest of the flights. Above approximately 35m AGL, the alignement will be consistent. However, if images taken closer to the ground are used, such as panel images, the same alignment transformation will not work for the flight data.
- -import cv2
-import numpy as np
-import matplotlib.pyplot as plt
-import micasense.imageutils as imageutils
-import micasense.plotutils as plotutils
-
-## Alignment settings
-match_index = 1 # Index of the band
-max_alignment_iterations = 10
-warp_mode = cv2.MOTION_HOMOGRAPHY # MOTION_HOMOGRAPHY or MOTION_AFFINE. For Altum images only use HOMOGRAPHY
-pyramid_levels = 0 # for images with RigRelatives, setting this to 0 or 1 may improve alignment
-
-print("Aligning images. Depending on settings this can take from a few seconds to many minutes")
-# Can potentially increase max_iterations for better results, but longer runtimes
-warp_matrices, alignment_pairs = imageutils.align_capture(capture,
- ref_index = match_index,
- max_iterations = max_alignment_iterations,
- warp_mode = warp_mode,
- pyramid_levels = pyramid_levels)
-
-print("Finished Aligning, warp matrices={}".format(warp_matrices))
-Aligning images. Depending on settings this can take from a few seconds to many minutes -Finished aligning band 1 -Finished aligning band 2 -Finished aligning band 4 -Finished aligning band 0 -Finished aligning band 3 -Finished Aligning, warp matrices=[array([[ 1.0109243e+00, 1.3997733e-03, 7.1472993e+00], - [ 2.3857178e-03, 1.0134553e+00, -1.1758372e+01], - [ 1.2079964e-06, 4.8187062e-06, 1.0000000e+00]], dtype=float32), array([[1.0000000e+00, 0.0000000e+00, 1.1368684e-13], - [0.0000000e+00, 1.0000000e+00, 0.0000000e+00], - [0.0000000e+00, 0.0000000e+00, 1.0000000e+00]], dtype=float32), array([[ 1.0093619e+00, 3.3198819e-03, -3.2410549e+01], - [ 8.4111997e-04, 1.0132477e+00, 1.3615667e+01], - [-4.9676191e-07, 7.0111369e-06, 1.0000000e+00]], dtype=float32), array([[ 1.0142691e+00, 6.4295256e-03, -2.3179104e+01], - [ 4.6660731e-05, 1.0160753e+00, -1.2280207e+01], - [ 2.7585847e-06, 7.7613631e-06, 1.0000000e+00]], dtype=float32), array([[ 1.0161121e+00, 2.1327983e-03, -2.0565905e+01], - [ 1.3823286e-03, 1.0168507e+00, -6.6128030e+00], - [ 2.2218899e-06, 3.5776138e-06, 1.0000000e+00]], dtype=float32)] --
After finding image alignments we may need to remove pixels around the edges which aren't present in every image in the capture. To do this we use the affine transforms found above and the image distortions from the image metadata. OpenCV provides a couple of handy helpers for this task in the cv2.undistortPoints() and cv2.transform() methods. These methods takes a set of pixel coordinates and apply our undistortion matrix and our affine transform, respectively. So, just as we did when registering the images, we first apply the undistortion process the coordinates of the image borders, then we apply the affine transformation to that result. The resulting pixel coordinates tell us where the image borders end up after this pair of transformations, and we can then crop the resultant image to these coordinates.
cropped_dimensions, edges = imageutils.find_crop_bounds(capture, warp_matrices, warp_mode=warp_mode, reference_band=match_index)
-print(cropped_dimensions)
-im_aligned = imageutils.aligned_capture(capture, warp_matrices, warp_mode, cropped_dimensions, match_index, img_type=img_type)
-(35.0, 19.0, 1222.0, 912.0) --
Once the transformation has been found, it can be verified by composting the aligned images to check alignment. The image 'stack' containing all bands can also be exported to a multi-band TIFF file for viewing in extrernal software such as QGIS. Usef ul componsites are a naturally colored RGB as well as color infrared, or CIR.
- -# figsize=(30,23) # use this size for full-image-resolution display
-figsize=(16,13) # use this size for export-sized display
-
-rgb_band_indices = [capture.band_names_lower().index('red'),
- capture.band_names_lower().index('green'),
- capture.band_names_lower().index('blue')]
-cir_band_indices = [capture.band_names_lower().index('nir'),
- capture.band_names_lower().index('red'),
- capture.band_names_lower().index('green')]
-
-# Create a normalized stack for viewing
-im_display = np.zeros((im_aligned.shape[0],im_aligned.shape[1],im_aligned.shape[2]), dtype=np.float32 )
-
-im_min = np.percentile(im_aligned[:,:,rgb_band_indices].flatten(), 0.5) # modify these percentiles to adjust contrast
-im_max = np.percentile(im_aligned[:,:,rgb_band_indices].flatten(), 99.5) # for many images, 0.5 and 99.5 are good values
-
-# for rgb true color, we use the same min and max scaling across the 3 bands to
-# maintain the "white balance" of the calibrated image
-for i in rgb_band_indices:
- im_display[:,:,i] = imageutils.normalize(im_aligned[:,:,i], im_min, im_max)
-
-rgb = im_display[:,:,rgb_band_indices]
-
-# for cir false color imagery, we normalize the NIR,R,G bands within themselves, which provides
-# the classical CIR rendering where plants are red and soil takes on a blue tint
-for i in cir_band_indices:
- im_display[:,:,i] = imageutils.normalize(im_aligned[:,:,i])
-
-cir = im_display[:,:,cir_band_indices]
-fig, axes = plt.subplots(1, 2, figsize=figsize)
-axes[0].set_title("Red-Green-Blue Composite")
-axes[0].imshow(rgb)
-axes[1].set_title("Color Infrared (CIR) Composite")
-axes[1].imshow(cir)
-plt.show()
-There are many techniques for image enhancement, but one which is commonly used to improve the visual sharpness of imagery is the unsharp mask. Here we apply an unsharp mask to the RGB image to improve the visualization, and then apply a gamma curve to make the darkest areas brighter.
- -# Create an enhanced version of the RGB render using an unsharp mask
-gaussian_rgb = cv2.GaussianBlur(rgb, (9,9), 10.0)
-gaussian_rgb[gaussian_rgb<0] = 0
-gaussian_rgb[gaussian_rgb>1] = 1
-unsharp_rgb = cv2.addWeighted(rgb, 1.5, gaussian_rgb, -0.5, 0)
-unsharp_rgb[unsharp_rgb<0] = 0
-unsharp_rgb[unsharp_rgb>1] = 1
-
-# Apply a gamma correction to make the render appear closer to what our eyes would see
-gamma = 1.4
-gamma_corr_rgb = unsharp_rgb**(1.0/gamma)
-fig = plt.figure(figsize=figsize)
-plt.imshow(gamma_corr_rgb, aspect='equal')
-plt.axis('off')
-plt.show()
-Composite images can be exported to JPEG or PNG format using the imageio package. These images may be useful for visualization or thumbnailing, and creating RGB thumbnails of a set of images can provide a convenient way to browse the imagery in a more visually appealing way that browsing the raw imagery.
import imageio
-imtype = 'png' # or 'jpg'
-imageio.imwrite('rgb.'+imtype, (255*gamma_corr_rgb).astype('uint8'))
-imageio.imwrite('cir.'+imtype, (255*cir).astype('uint8'))
-We can export the image easily stacks using the gdal library (http://www.glal.org). Once exported, these image stacks can be opened in software such as QGIS and raster operations such as NDVI or NDRE computation can be done in that software. At this time the stacks don't include any geographic information.
from osgeo import gdal, gdal_array
-rows, cols, bands = im_display.shape
-driver = gdal.GetDriverByName('GTiff')
-filename = "bgrne" #blue,green,red,nir,redEdge
-
-if im_aligned.shape[2] == 6:
- filename = filename + "t" #thermal
-outRaster = driver.Create(filename+"B.tiff", cols, rows, im_aligned.shape[2], gdal.GDT_UInt16)
-
-normalize = (img_type == 'radiance') # normalize radiance images to fit with in UInt16
-
-# Output a 'stack' in the same band order as RedEdge/Alutm
-# Blue,Green,Red,NIR,RedEdge[,Thermal]
-# reflectance stacks are output with 32768=100% reflectance to provide some overhead for specular reflections
-# radiance stacks are output with 65535=100% radiance to provide some overhead for specular reflections
-
-# NOTE: NIR and RedEdge are not in wavelength order!
-
-multispec_min = np.min(im_aligned[:,:,1:5])
-multispec_max = np.max(im_aligned[:,:,1:5])
-
-for i in range(0,5):
- outband = outRaster.GetRasterBand(i+1)
- if normalize:
- outdata = imageutils.normalize(im_aligned[:,:,i],multispec_min,multispec_max)
- else:
- outdata = im_aligned[:,:,i]
- outdata[outdata<0] = 0
- outdata[outdata>2] = 2
-
- outdata = outdata*32767
- outdata[outdata<0] = 0
- outdata[outdata>65535] = 65535
- outband.WriteArray(outdata)
- outband.FlushCache()
-
-if im_aligned.shape[2] == 6:
- outband = outRaster.GetRasterBand(6)
- outdata = im_aligned[:,:,5] * 100 # scale to centi-C to fit into uint16
- outdata[outdata<0] = 0
- outdata[outdata>65535] = 65535
- outband.WriteArray(outdata)
- outband.FlushCache()
-outRaster = None
-"Stacks" as described above are useful in a number of processing cases. For example, at the time of this writing, many photogrammetry suites could import and process stack files without significantly impacting the radiometric processing which has already been accomplished.
-Running photogrammetry on stack files instead of raw image files has both advantages and drawbacks. The primary advantage has been found to be an increase in processing speed and a reduction in program memory usage. As the photogrammetric workflow generally operates on luminance images and may not use color information, stacked images may require similar resources and be processed at a similar speed as single-band images. This is because one band of the stack can be used to generate the matching feature space while the others are ignored for matching purposes. This reduces the feature space 5-fold over matching using all images separately.
-One disadvantage is that stacking images outside of the photogrammetric workflow may result in poor image matching. The RedEdge is known to have stable lens characteristics over the course of normal operation, but variations in temperature or impacts to the camera through handling or rough landings may change the image alignment parameters. For this reason, we recommend finding a matching transformation for each flight (each take-off and landing). Alignment transformations from multiple images within a flight can be compared to find the best transformation to apply to the set of the flight. While not described or supported in this generic implementation, some matching algorithms can use a "seed" value as a starting point to speed up matching. For most cases, this seed could be the transformation found in a previous flight, or another source of a known good transformation.
- -For raw index computation on single images, the numpy package provides a simple way to do math and simple visualizatoin on images. Below, we compute and visualize an image histogram and then use that to pick a colormap range for visualizing the NDVI of an image.
-After computing the NDVI and prior to displaying it, we use a very rudimentary method for focusing on the plants and removing the soil and shadow information from our images and histograms. Below we remove non-plant pixels by setting to zero any pixels in the image where the NIR reflectance is less than 20%. This helps to ensure that the NDVI and NDRE histograms aren't skewed substantially by soil noise.
- -from micasense import plotutils
-import matplotlib.pyplot as plt
-
-nir_band = capture.band_names_lower().index('nir')
-red_band = capture.band_names_lower().index('red')
-
-np.seterr(divide='ignore', invalid='ignore') # ignore divide by zero errors in the index calculation
-
-# Compute Normalized Difference Vegetation Index (NDVI) from the NIR(3) and RED (2) bands
-ndvi = (im_aligned[:,:,nir_band] - im_aligned[:,:,red_band]) / (im_aligned[:,:,nir_band] + im_aligned[:,:,red_band])
-
-# remove shadowed areas (mask pixels with NIR reflectance < 20%))
-if img_type == 'reflectance':
- ndvi = np.ma.masked_where(im_aligned[:,:,nir_band] < 0.20, ndvi)
-elif img_type == 'radiance':
- lower_pct_radiance = np.percentile(im_aligned[:,:,3], 10.0)
- ndvi = np.ma.masked_where(im_aligned[:,:,nir_band] < lower_pct_radiance, ndvi)
-
-# Compute and display a histogram
-ndvi_hist_min = np.min(ndvi)
-ndvi_hist_max = np.max(ndvi)
-fig, axis = plt.subplots(1, 1, figsize=(10,4))
-axis.hist(ndvi.ravel(), bins=512, range=(ndvi_hist_min, ndvi_hist_max))
-plt.title("NDVI Histogram")
-plt.show()
-
-min_display_ndvi = 0.45 # further mask soil by removing low-ndvi values
-#min_display_ndvi = np.percentile(ndvi.flatten(), 5.0) # modify with these percentilse to adjust contrast
-max_display_ndvi = np.percentile(ndvi.flatten(), 99.5) # for many images, 0.5 and 99.5 are good values
-masked_ndvi = np.ma.masked_where(ndvi < min_display_ndvi, ndvi)
-
-#reduce the figure size to account for colorbar
-figsize=np.asarray(figsize) - np.array([3,2])
-
-#plot NDVI over an RGB basemap, with a colorbar showing the NDVI scale
-fig, axis = plotutils.plot_overlay_withcolorbar(gamma_corr_rgb,
- masked_ndvi,
- figsize = figsize,
- title = 'NDVI filtered to only plants over RGB base layer',
- vmin = min_display_ndvi,
- vmax = max_display_ndvi)
-fig.savefig('ndvi_over_rgb.png')
-In the same manner, we can compute, filter, and display another index useful for the RedEdge camera, the Normalized Difference Red Edge (NDRE) index. We also filter out shadows and soil to ensure our display focuses only on the plant health.
- -# Compute Normalized Difference Red Edge Index from the NIR(3) and RedEdge(4) bands
-rededge_band = capture.band_names_lower().index('red edge')
-ndre = (im_aligned[:,:,nir_band] - im_aligned[:,:,rededge_band]) / (im_aligned[:,:,nir_band] + im_aligned[:,:,rededge_band])
-
-# Mask areas with shadows and low NDVI to remove soil
-masked_ndre = np.ma.masked_where(ndvi < min_display_ndvi, ndre)
-
-# Compute a histogram
-ndre_hist_min = np.min(masked_ndre)
-ndre_hist_max = np.max(masked_ndre)
-fig, axis = plt.subplots(1, 1, figsize=(10,4))
-axis.hist(masked_ndre.ravel(), bins=512, range=(ndre_hist_min, ndre_hist_max))
-plt.title("NDRE Histogram (filtered to only plants)")
-plt.show()
-
-min_display_ndre = np.percentile(masked_ndre, 5)
-max_display_ndre = np.percentile(masked_ndre, 99.5)
-
-fig, axis = plotutils.plot_overlay_withcolorbar(gamma_corr_rgb,
- masked_ndre,
- figsize=figsize,
- title='NDRE filtered to only plants over RGB base layer',
- vmin=min_display_ndre,vmax=max_display_ndre)
-fig.savefig('ndre_over_rgb.png')
-If our image is from an Altum and includes a thermal band, we can display the re-sampled and aligned thermal data over the RGB data to maintain the context of the thermal information.
-In the image below, it's very clear based on the average temperature where the soil is wet and dry, and even in the middle of the road we can find some wetter areas of soil.
- -if len(capture.lw_indices()) > 0:
-
- # by default we don't mask the thermal, since it's native resolution is much lower than the MS
- masked_thermal = im_aligned[:,:,5]
- # Alternatively we can mask the thermal only to plants here, which is useful for large contiguous areas
- # masked_thermal = np.ma.masked_where(ndvi < 0.45, im_aligned[:,:,5])
-
-
- # Compute a histogram
- fig, axis = plt.subplots(1, 1, figsize=(10,4))
- axis.hist(masked_thermal.ravel(), bins=512, range=(np.min(masked_thermal), np.max(masked_thermal)))
- plt.title("Thermal Histogram")
- plt.show()
-
- min_display_therm = np.percentile(masked_thermal, 1)
- max_display_therm = np.percentile(masked_thermal, 99)
-
- fig, axis = plotutils.plot_overlay_withcolorbar(gamma_corr_rgb,
- masked_thermal,
- figsize=figsize,
- title='Temperature over True Color',
- vmin=min_display_therm,vmax=max_display_therm,
- overlay_alpha=0.25,
- overlay_colormap='jet',
- overlay_steps=16,
- display_contours=True,
- contour_steps=16,
- contour_alpha=.4,
- contour_fmt="%.0fC")
- fig.savefig('thermal_over_rgb.png')
-Finally, we show a classic agricultural remote sensing output in the tassled cap plot. This plot can be useful for visualizing row crops and plots the Red Reflectance channel on the X-axis against the NIR reflectance channel on the Y-axis. This plot also clearly shows the line of the soil in that space. The tassled cap view isn't very useful for this arid data set; however, we can see the "badge of trees" of high NIR reflectance and relatively low red reflectance. This provides an example of one of the uses of aligned images for single capture analysis.
- -x_band = red_band
-y_band = nir_band
-x_max = np.max(im_aligned[:,:,x_band])
-y_max = np.max(im_aligned[:,:,y_band])
-
-fig = plt.figure(figsize=(12,12))
-plt.hexbin(im_aligned[:,:,x_band],im_aligned[:,:,y_band],gridsize=640,bins='log',extent=(0,x_max,0,y_max))
-ax = fig.gca()
-ax.set_xlim([0,x_max])
-ax.set_ylim([0,y_max])
-plt.xlabel("{} Reflectance".format(capture.band_names()[x_band]))
-plt.ylabel("{} Reflectance".format(capture.band_names()[y_band]))
-plt.show()
-Last, we output the warp_matrices that we got for this image stack for usage elsewhere. Currently these can be used in the Batch Processing.ipynb notebook to save reflectance-compensated stacks of images to a directory.
print(warp_matrices)
-[array([[ 1.0109243e+00, 1.3997733e-03, 7.1472993e+00], - [ 2.3857178e-03, 1.0134553e+00, -1.1758372e+01], - [ 1.2079964e-06, 4.8187062e-06, 1.0000000e+00]], dtype=float32), array([[1.0000000e+00, 0.0000000e+00, 1.1368684e-13], - [0.0000000e+00, 1.0000000e+00, 0.0000000e+00], - [0.0000000e+00, 0.0000000e+00, 1.0000000e+00]], dtype=float32), array([[ 1.0093619e+00, 3.3198819e-03, -3.2410549e+01], - [ 8.4111997e-04, 1.0132477e+00, 1.3615667e+01], - [-4.9676191e-07, 7.0111369e-06, 1.0000000e+00]], dtype=float32), array([[ 1.0142691e+00, 6.4295256e-03, -2.3179104e+01], - [ 4.6660731e-05, 1.0160753e+00, -1.2280207e+01], - [ 2.7585847e-06, 7.7613631e-06, 1.0000000e+00]], dtype=float32), array([[ 1.0161121e+00, 2.1327983e-03, -2.0565905e+01], - [ 1.3823286e-03, 1.0168507e+00, -6.6128030e+00], - [ 2.2218899e-06, 3.5776138e-06, 1.0000000e+00]], dtype=float32)] --
Copyright (c) 2017-2018 MicaSense, Inc. For licensing information see the project git repository
- -In this example, we use the micasense.imageset class to load a set of directories of images into a list of micasense.capture objects, and we iterate over that list, saving out each image as an aligned stack of images as separate bands in a single tiff file each. Part of this process (via imageutils.write_exif_to_stack) injects that the GPS, capture datetime, camera model, etc into the processed images, allowing us to stitch those images using commercial software such as Pix4DMapper or Agisoft Metashape.
Note: for this example to work, the images must have a valid RigRelatives tag. This requires RedEdge (3/M/MX) version of at least 3.4.0, or any version of RedEdge-P/Altum-PT/Altum/RedEdge-MX Dual. If your images don't meet that spec, you can also follow this support article to add the RigRelatives tag to your imagery: https://support.micasense.com/hc/en-us/articles/360006368574-Modifying-older-collections-for-Pix4Dfields-support
- -%load_ext autoreload
-%autoreload 2
-from ipywidgets import FloatProgress, Layout
-from IPython.display import display
-import micasense.imageset as imageset
-import micasense.capture as capture
-import os, glob
-import multiprocessing
-from pathlib import Path
-
-# set to True if you have an Altum-PT
-# or RedEdge-P and wish to output pan-sharpened stacks
-panSharpen = True
-
-panelNames = None
-useDLS = True
-
-# set your image path here. See more here: https://docs.python.org/3/library/pathlib.html
-imagePath = Path("./data/REDEDGE-MX")
-
-# these will return lists of image paths as strings
-panelNames = list(imagePath.glob('IMG_0001_*.tif'))
-panelNames = [x.as_posix() for x in panelNames]
-
-panelCap = capture.Capture.from_filelist(panelNames)
-
-# destinations on your computer to put the stacks
-# and RGB thumbnails
-outputPath = imagePath / '..' / 'stacks'
-thumbnailPath = outputPath / 'thumbnails'
-
-cam_model = panelCap.camera_model
-cam_serial = panelCap.camera_serial
-
-# determine if this sensor has a panchromatic band
-if cam_model == 'RedEdge-P' or cam_model == 'Altum-PT':
- panchroCam = True
-else:
- panchroCam = False
- panSharpen = False
-
-# if this is a multicamera system like the RedEdge-MX Dual,
-# we can combine the two serial numbers to help identify
-# this camera system later.
-if len(panelCap.camera_serials) > 1:
- cam_serial = "_".join(panelCap.camera_serials)
- print("Serial number:",cam_serial)
-else:
- cam_serial = panelCap.camera_serial
- print("Serial number:",cam_serial)
-
-overwrite = False # can be set to set to False to continue interrupted processing
-generateThumbnails = True
-
-# Allow this code to align both radiance and reflectance images; but excluding
-# a definition for panelNames above, radiance images will be used
-# For panel images, efforts will be made to automatically extract the panel information
-# but if the panel/firmware is before Altum 1.3.5, RedEdge 5.1.7 the panel reflectance
-# will need to be set in the panel_reflectance_by_band variable.
-# Note: radiance images will not be used to properly create NDVI/NDRE images below.
-if panelNames is not None:
- panelCap = capture.Capture.from_filelist(panelNames)
-else:
- panelCap = None
-
-if panelCap is not None:
- if panelCap.panel_albedo() is not None and not any(v is None for v in panelCap.panel_albedo()):
- panel_reflectance_by_band = panelCap.panel_albedo()
- else:
- panel_reflectance_by_band = [0.49]*len(panelCap.eo_band_names()) #RedEdge band_index order
-
- panel_irradiance = panelCap.panel_irradiance(panel_reflectance_by_band)
- img_type = "reflectance"
-else:
- if useDLS:
- img_type='reflectance'
- else:
- img_type = "radiance"
-Serial number: RX02-2023065-SC --
## This progress widget is used for display of the long-running process
-f = FloatProgress(min=0, max=1, layout=Layout(width='100%'), description="Loading")
-display(f)
-def update_f(val):
- if (val - f.value) > 0.005 or val == 1: #reduces cpu usage from updating the progressbar by 10x
- f.value=val
-
-%time imgset = imageset.ImageSet.from_directory(imagePath, progress_callback=update_f)
-update_f(1.0)
-CPU times: user 31.5 ms, sys: 7.14 ms, total: 38.6 ms -Wall time: 278 ms --
We can map out the capture GPS locations to ensure we are processing the right data. A GeoJSON of the captures will later be saved to the outputPath.
- -import math
-import numpy as np
-from mapboxgl.viz import *
-from mapboxgl.utils import df_to_geojson, create_radius_stops, scale_between
-from mapboxgl.utils import create_color_stops
-import pandas as pd
-
-data, columns = imgset.as_nested_lists()
-df = pd.DataFrame.from_records(data, index='timestamp', columns=columns)
-
-#Insert your mapbox token here
-token = 'pk.eyJ1IjoibWljYXNlbnNlIiwiYSI6ImNqYWx5dWNteTJ3cWYzMnBicmZid3g2YzcifQ.Zrq9t7GYocBtBzYyT3P4sw'
-color_property = 'dls-yaw'
-num_color_classes = 8
-
-min_val = df[color_property].min()
-max_val = df[color_property].max()
-
-import jenkspy
-# breaks = jenkspy.jenks_breaks(df[color_property], nb_class=num_color_classes)
-
-# color_stops = create_color_stops(breaks,colors='YlOrRd')
-geojson_data = df_to_geojson(df,columns[3:],lat='latitude',lon='longitude')
-
-viz = CircleViz(geojson_data, access_token=token, color_property=color_property,
-# color_stops=color_stops,
- center=[df['longitude'].median(),df['latitude'].median()],
- zoom=16, height='600px',
- style='mapbox://styles/mapbox/satellite-streets-v9')
-viz.show()
-/Users/stephen/anaconda3/envs/micasense/lib/python3.7/site-packages/IPython/core/display.py:724: UserWarning: Consider using IPython.display.IFrame instead
- warnings.warn("Consider using IPython.display.IFrame instead")
-
-For newer data sets with RigRelatives tags (images captured with RedEdge (3/M/MX) version 3.4.0 or greater with a valid calibration load, see https://support.micasense.com/hc/en-us/articles/360005428953-Updating-RedEdge-for-Pix4Dfields), we can use the RigRelatives for a simple alignment. To use this simple alignment, simply set warp_matrices=None
For sets without those tags, or sets that require a RigRelatives optimization, we can go through the Alignment.ipynb notebook and get a set of warp_matrices that we can use here to align.
- -from numpy import array
-from numpy import float32
-from skimage.transform import ProjectiveTransform
-
-if panchroCam:
- warp_matrices_filename = cam_serial + "_warp_matrices_SIFT.npy"
-else:
- warp_matrices_filename = cam_serial + "_warp_matrices_opencv.npy"
-
-if Path('./' + warp_matrices_filename).is_file():
- print("Found existing warp matrices for camera", cam_serial)
- load_warp_matrices = np.load(warp_matrices_filename, allow_pickle=True)
- loaded_warp_matrices = []
- for matrix in load_warp_matrices:
- if panchroCam:
- transform = ProjectiveTransform(matrix=matrix.astype('float64'))
- loaded_warp_matrices.append(transform)
- else:
- loaded_warp_matrices.append(matrix.astype('float32'))
-
- if panchroCam:
- warp_matrices_SIFT = loaded_warp_matrices
- else:
- warp_matrices = loaded_warp_matrices
- print("Loaded warp matrices from",Path('./' + warp_matrices_filename).resolve())
-else:
- print("No warp matrices found at expected location:",warp_matrices_filename)
-
-Found existing warp matrices for camera RX02-2023065-SC -Loaded warp matrices from /Users/stephen/programs/gitlab/imageProcessingTutorial/RX02-2023065-SC_warp_matrices_opencv.npy --
import exiftool
-import datetime
-## This progress widget is used for display of the long-running process
-f2 = FloatProgress(min=0, max=1, layout=Layout(width='100%'), description="Saving")
-display(f2)
-def update_f2(val):
- f2.value=val
-
-if not os.path.exists(outputPath):
- os.makedirs(outputPath)
-if generateThumbnails and not os.path.exists(thumbnailPath):
- os.makedirs(thumbnailPath)
-
-# Save out geojson data so we can open the image capture locations in our GIS
-with open(os.path.join(outputPath,'imageSet.json'),'w') as f:
- f.write(str(geojson_data))
-
-try:
- irradiance = panel_irradiance+[0]
-except NameError:
- irradiance = None
-
-start = datetime.datetime.now()
-for i,capture in enumerate(imgset.captures):
- outputFilename = capture.uuid+'.tif'
- thumbnailFilename = capture.uuid+'.jpg'
- fullOutputPath = os.path.join(outputPath, outputFilename)
- fullThumbnailPath= os.path.join(thumbnailPath, thumbnailFilename)
- if (not os.path.exists(fullOutputPath)) or overwrite:
- if(len(capture.images) == len(imgset.captures[0].images)):
- if panchroCam:
- capture.radiometric_pan_sharpened_aligned_capture(warp_matrices=warp_matrices_SIFT,irradiance_list=irradiance)
- else:
- capture.create_aligned_capture(irradiance_list=irradiance, warp_matrices=warp_matrices)
- capture.save_capture_as_stack(fullOutputPath, pansharpen=panSharpen,sort_by_wavelength=False)
- if generateThumbnails:
- capture.save_capture_as_rgb(fullThumbnailPath)
- capture.clear_image_data()
- update_f2(float(i)/float(len(imgset.captures)))
-update_f2(1.0)
-end = datetime.datetime.now()
-
-print("Saving time: {}".format(end-start))
-print("Alignment+Saving rate: {:.2f} images per second".format(float(len(imgset.captures))/float((end-start).total_seconds())))
-Saving time: 0:00:04.008601 -Alignment+Saving rate: 0.50 images per second --
-In this example, we use the micasense.imageset class to load a set of directories of images into a list of micasense.capture objects, and we iterate over that list saving out each image as an aligned stack of images as separate bands in a single tiff file each. Next, we use the metadata from the original captures to write out a log file of the captures and their locations. Finally, we use exiftool from the command line to inject that metadata into the processed images, allowing us to stitch those images using commercial software such as Pix4D or Agisoft.
Note: for this example to work, the images must have a valid RigRelatives tag. This requires RedEdge version of at least 3.4.0 or any version of Altum. If your images don't meet that spec, you can also follow this support ticket to add the RigRelatives tag to them: https://support.micasense.com/hc/en-us/articles/360006368574-Modifying-older-collections-for-Pix4Dfields-support
- -%load_ext autoreload
-%autoreload 2
-The autoreload extension is already loaded. To reload it, use: - %reload_ext autoreload --
from ipywidgets import FloatProgress, Layout
-from IPython.display import display
-import micasense.imageset as imageset
-import micasense.capture as capture
-import os, glob
-import multiprocessing
-from pathlib import Path
-
-panelNames = None
-useDLS = True
-
-# set your image path here. See more here: https://docs.python.org/3/library/pathlib.html
-imagePath = Path("./data/REDEDGE-MX")
-
-# these will return lists of image paths as strings
-panelNames = list(imagePath.glob('IMG_0001_*.tif'))
-panelNames = [x.as_posix() for x in panelNames]
-
-panelCap = capture.Capture.from_filelist(panelNames)
-
-# destinations on your computer to put the stacks
-# and RGB thumbnails
-outputPath = imagePath / '..' / 'stacks'
-thumbnailPath = outputPath / 'thumbnails'
-
-overwrite = False # can be set to set to False to continue interrupted processing
-generateThumbnails = True
-
-# Allow this code to align both radiance and reflectance images; bu excluding
-# a definition for panelNames above, radiance images will be used
-# For panel images, efforts will be made to automatically extract the panel information
-# but if the panel/firmware is before Altum 1.3.5, RedEdge 5.1.7 the panel reflectance
-# will need to be set in the panel_reflectance_by_band variable.
-# Note: radiance images will not be used to properly create NDVI/NDRE images below.
-if panelNames is not None:
- panelCap = capture.Capture.from_filelist(panelNames)
-else:
- panelCap = None
-
-if panelCap is not None:
- if panelCap.panel_albedo() is not None and not any(v is None for v in panelCap.panel_albedo()):
- panel_reflectance_by_band = panelCap.panel_albedo()
- else:
- panel_reflectance_by_band = [0.49, 0.49, 0.49, 0.49, 0.49] #RedEdge band_index order
-
- panel_irradiance = panelCap.panel_irradiance(panel_reflectance_by_band)
- img_type = "reflectance"
-else:
- if useDLS:
- img_type='reflectance'
- else:
- img_type = "radiance"
-## This progress widget is used for display of the long-running process
-f = FloatProgress(min=0, max=1, layout=Layout(width='100%'), description="Loading")
-display(f)
-def update_f(val):
- if (val - f.value) > 0.005 or val == 1: #reduces cpu usage from updating the progressbar by 10x
- f.value=val
-
-%time imgset = imageset.ImageSet.from_directory(imagePath, progress_callback=update_f)
-update_f(1.0)
-FloatProgress(value=0.0, description='Loading', layout=Layout(width='100%'), max=1.0)-
CPU times: user 32.4 ms, sys: 10.3 ms, total: 42.8 ms -Wall time: 315 ms --
import math
-import numpy as np
-from mapboxgl.viz import *
-from mapboxgl.utils import df_to_geojson, create_radius_stops, scale_between
-from mapboxgl.utils import create_color_stops
-import pandas as pd
-
-data, columns = imgset.as_nested_lists()
-df = pd.DataFrame.from_records(data, index='timestamp', columns=columns)
-
-#Insert your mapbox token here
-token = 'pk.eyJ1IjoibWljYXNlbnNlIiwiYSI6ImNqYWx5dWNteTJ3cWYzMnBicmZid3g2YzcifQ.Zrq9t7GYocBtBzYyT3P4sw'
-color_property = 'dls-yaw'
-num_color_classes = 8
-
-min_val = df[color_property].min()
-max_val = df[color_property].max()
-
-import jenkspy
-# breaks = jenkspy.jenks_breaks(df[color_property], nb_class=num_color_classes)
-
-# color_stops = create_color_stops(breaks,colors='YlOrRd')
-geojson_data = df_to_geojson(df,columns[3:],lat='latitude',lon='longitude')
-
-viz = CircleViz(geojson_data, access_token=token, color_property=color_property,
-# color_stops=color_stops,
- center=[df['longitude'].median(),df['latitude'].median()],
- zoom=16, height='600px',
- style='mapbox://styles/mapbox/satellite-streets-v9')
-viz.show()
-/Users/stephen/anaconda3/envs/micasense/lib/python3.7/site-packages/IPython/core/display.py:724: UserWarning: Consider using IPython.display.IFrame instead
- warnings.warn("Consider using IPython.display.IFrame instead")
-
-For newer data sets with RigRelatives tags (images captured with RedEdge version 3.4.0 or greater with a valid calibration load, see https://support.micasense.com/hc/en-us/articles/360005428953-Updating-RedEdge-for-Pix4Dfields), we can use the RigRelatives for a simple alignment.
-For sets without those tags, or sets that require a RigRelatives optimization, we can go through the Alignment.ipynb notebook and get a set of warp_matrices that we can use here to align.
from numpy import array
-from numpy import float32
-
-# Set warp_matrices to none to align using RigRelatives
-# Or
-# Use the warp_matrices derived from the Alignment Tutorial for this RedEdge set without RigRelatives
-warp_matrices = [array([[ 1.0109243e+00, 1.3997733e-03, 7.1472993e+00],
- [ 2.3857178e-03, 1.0134553e+00, -1.1758372e+01],
- [ 1.2079964e-06, 4.8187062e-06, 1.0000000e+00]], dtype=float32), array([[1.0000000e+00, 0.0000000e+00, 1.1368684e-13],
- [0.0000000e+00, 1.0000000e+00, 0.0000000e+00],
- [0.0000000e+00, 0.0000000e+00, 1.0000000e+00]], dtype=float32), array([[ 1.0093619e+00, 3.3198819e-03, -3.2410549e+01],
- [ 8.4111997e-04, 1.0132477e+00, 1.3615667e+01],
- [-4.9676191e-07, 7.0111369e-06, 1.0000000e+00]], dtype=float32), array([[ 1.0142691e+00, 6.4295256e-03, -2.3179104e+01],
- [ 4.6660731e-05, 1.0160753e+00, -1.2280207e+01],
- [ 2.7585847e-06, 7.7613631e-06, 1.0000000e+00]], dtype=float32), array([[ 1.0161121e+00, 2.1327983e-03, -2.0565905e+01],
- [ 1.3823286e-03, 1.0168507e+00, -6.6128030e+00],
- [ 2.2218899e-06, 3.5776138e-06, 1.0000000e+00]], dtype=float32)]
-import exiftool
-import datetime
-## This progress widget is used for display of the long-running process
-f2 = FloatProgress(min=0, max=1, layout=Layout(width='100%'), description="Saving")
-display(f2)
-def update_f2(val):
- f2.value=val
-
-if not os.path.exists(outputPath):
- os.makedirs(outputPath)
-if generateThumbnails and not os.path.exists(thumbnailPath):
- os.makedirs(thumbnailPath)
-
-# Save out geojson data so we can open the image capture locations in our GIS
-with open(os.path.join(outputPath,'imageSet.json'),'w') as f:
- f.write(str(geojson_data))
-
-try:
- irradiance = panel_irradiance+[0]
-except NameError:
- irradiance = None
-
-start = datetime.datetime.now()
-for i,capture in enumerate(imgset.captures):
- outputFilename = capture.uuid+'.tif'
- thumbnailFilename = capture.uuid+'.jpg'
- fullOutputPath = os.path.join(outputPath, outputFilename)
- fullThumbnailPath= os.path.join(thumbnailPath, thumbnailFilename)
- if (not os.path.exists(fullOutputPath)) or overwrite:
- if(len(capture.images) == len(imgset.captures[0].images)):
- capture.create_aligned_capture(irradiance_list=irradiance, warp_matrices=warp_matrices)
- capture.save_capture_as_stack(fullOutputPath)
- if generateThumbnails:
- capture.save_capture_as_rgb(fullThumbnailPath)
- capture.clear_image_data()
- update_f2(float(i)/float(len(imgset.captures)))
-update_f2(1.0)
-end = datetime.datetime.now()
-
-print("Saving time: {}".format(end-start))
-print("Alignment+Saving rate: {:.2f} images per second".format(float(len(imgset.captures))/float((end-start).total_seconds())))
-FloatProgress(value=0.0, description='Saving', layout=Layout(width='100%'), max=1.0)-
Saving time: 0:00:05.053582 -Alignment+Saving rate: 0.40 images per second --
def decdeg2dms(dd):
- is_positive = dd >= 0
- dd = abs(dd)
- minutes,seconds = divmod(dd*3600,60)
- degrees,minutes = divmod(minutes,60)
- degrees = degrees if is_positive else -degrees
- return (degrees,minutes,seconds)
-
-header = "SourceFile,\
-GPSDateStamp,GPSTimeStamp,\
-GPSLatitude,GpsLatitudeRef,\
-GPSLongitude,GPSLongitudeRef,\
-GPSAltitude,GPSAltitudeRef,\
-FocalLength,\
-XResolution,YResolution,ResolutionUnits\n"
-
-lines = [header]
-for capture in imgset.captures:
- #get lat,lon,alt,time
- outputFilename = capture.uuid+'.tif'
- fullOutputPath = os.path.join(outputPath, outputFilename)
- lat,lon,alt = capture.location()
- #write to csv in format:
- # IMG_0199_1.tif,"33 deg 32' 9.73"" N","111 deg 51' 1.41"" W",526 m Above Sea Level
- latdeg, latmin, latsec = decdeg2dms(lat)
- londeg, lonmin, lonsec = decdeg2dms(lon)
- latdir = 'North'
- if latdeg < 0:
- latdeg = -latdeg
- latdir = 'South'
- londir = 'East'
- if londeg < 0:
- londeg = -londeg
- londir = 'West'
- resolution = capture.images[0].focal_plane_resolution_px_per_mm
-
- linestr = '"{}",'.format(fullOutputPath)
- linestr += capture.utc_time().strftime("%Y:%m:%d,%H:%M:%S,")
- linestr += '"{:d} deg {:d}\' {:.2f}"" {}",{},'.format(int(latdeg),int(latmin),latsec,latdir[0],latdir)
- linestr += '"{:d} deg {:d}\' {:.2f}"" {}",{},{:.1f} m Above Sea Level,Above Sea Level,'.format(int(londeg),int(lonmin),lonsec,londir[0],londir,alt)
- linestr += '{}'.format(capture.images[0].focal_length)
- linestr += '{},{},mm'.format(resolution,resolution)
- linestr += '\n' # when writing in text mode, the write command will convert to os.linesep
- lines.append(linestr)
-
-fullCsvPath = os.path.join(outputPath,'log.csv')
-with open(fullCsvPath, 'w') as csvfile: #create CSV
- csvfile.writelines(lines)
-import subprocess
-
-if os.environ.get('exiftoolpath') is not None:
- exiftool_cmd = os.path.normpath(os.environ.get('exiftoolpath'))
-else:
- exiftool_cmd = 'exiftool'
-
-cmd = '{} -csv="{}" -overwrite_original "{}"'.format(exiftool_cmd, fullCsvPath, outputPath)
-print(cmd)
-subprocess.check_call(cmd, shell=True)
-exiftool -csv="data/REDEDGE-MX/../stacks/log.csv" -overwrite_original "data/REDEDGE-MX/../stacks" --
1 directories scanned - 2 image files updated --
0-
-A micasense.capture.Capture object holds a set of images (one per band) which are captured at the same moment together in a MicaSense camera. Files which meet this criteria will tend to have the same filename except for the suffix, but that is not required to load a captures. Captures can be loaded by starting with one image by calling Capture.from_file('file_name.tif') and adding others using the append_file, or by providing a list of filenames or images. See capture.py for more creation methods.
import os, glob
-import micasense.capture as capture
-%matplotlib inline
-
-imagePath = os.path.join('.','data','REDEDGE-MX')
-imageNames = glob.glob(os.path.join(imagePath,'IMG_0001_*.tif'))
-
-capture = capture.Capture.from_filelist(imageNames)
-capture.plot_raw()
-capture.plot_vignette();
-capture.plot_undistorted_radiance();
-capture.plot_panels()
-Copyright (c) 2017-2019 MicaSense, Inc. For licensing information see the project git repository
- -A large group of images captured using a RedEdge over a central California orchard are available for download here.
-With this set extracted to a working folder, this notebook will load all of the images in the set and provide more usage examples of ImageSet data.
-%load_ext autoreload
-%autoreload 2
-from ipywidgets import FloatProgress, Layout
-from IPython.display import display
-import micasense.imageset as imageset
-import os
-
-## This progress widget is used for display of the long-running process
-f = FloatProgress(min=0, max=1, layout=Layout(width='100%'), description="Loading")
-display(f)
-def update_f(val):
- if (val - f.value) > 0.005 or val == 1: #reduces cpu usage from updating the progressbar by 10x
- f.value=val
-
-images_dir = os.path.expanduser(os.path.join('.','data','REDEDGE-MX'))
-%time imgset = imageset.ImageSet.from_directory(images_dir, progress_callback=update_f)
-FloatProgress(value=0.0, description='Loading', layout=Layout(width='100%'), max=1.0)-
CPU times: user 33.4 ms, sys: 9.42 ms, total: 42.8 ms -Wall time: 304 ms --
import pandas as pd
-data, columns = imgset.as_nested_lists()
-print("Columns: {}".format(columns))
-df = pd.DataFrame.from_records(data, index='timestamp', columns=columns)
-Columns: ['timestamp', 'latitude', 'longitude', 'altitude', 'capture_id', 'dls-yaw', 'dls-pitch', 'dls-roll', 'irr-475', 'irr-560', 'irr-668', 'irr-842', 'irr-717'] --
Below we use the mapboxgl extension to plot the measured DLS yaw (or heading) angle from each image's meatadata over the whole imageset. We draw circles at each image location, and then color the circle based on the yaw value.
import math
-import numpy as np
-from mapboxgl.viz import *
-from mapboxgl.utils import df_to_geojson, create_radius_stops, scale_between
-from mapboxgl.utils import create_color_stops
-
-#Insert your mapbox token here
-token = 'pk.eyJ1IjoibWljYXNlbnNlIiwiYSI6ImNqYWx5dWNteTJ3cWYzMnBicmZid3g2YzcifQ.Zrq9t7GYocBtBzYyT3P4sw'
-color_stops = create_color_stops(np.linspace(-math.pi,math.pi,num=8),colors='YlOrRd')
-data = df_to_geojson(df,columns[3:],lat='latitude',lon='longitude')
-viz = CircleViz(data, access_token=token, color_property='dls-yaw',
- color_stops=color_stops,
- center=[df['longitude'].median(),df['latitude'].median()],
- zoom=16, height='600px',
- style='mapbox://styles/mapbox/satellite-streets-v9')
-viz.show()
-/Users/stephen/anaconda3/envs/micasense/lib/python3.7/site-packages/IPython/core/display.py:724: UserWarning: Consider using IPython.display.IFrame instead
- warnings.warn("Consider using IPython.display.IFrame instead")
-
-import matplotlib.pyplot as plt
-%matplotlib inline
-ax=df.plot(y=columns[3:], subplots=True, figsize=(16, 10), style=['g','c','y','k','b','g','r','m','k']);
-for a in ax:
- a.legend(loc='right', bbox_to_anchor=(1.1, 0.5), ncol=1, fancybox=True, shadow=True);
-We can plot a histogram of the image altitudes for all images in the flight and use the basic statistics to guess at an altitude below which the panel images were likely taken. This can give us an altitude threshold we can apply to separate images taken in flight and on the ground.
- -# plot the histogram of the altitude data
-df.altitude.hist();
-# find the altitude above which the flight images occur
-cutoff_altitude = df.altitude.mean()-3.0*df.altitude.std()
-plt.axvline(x=cutoff_altitude,c='r');
-plt.xlabel('Capture altitude (m)');
-plt.ylabel('Number of occurances');
-Using the Panel class, we can automatically find and compute the radiance of panel regions within panel images. Each Capture object exposes a panel_radiance() method which can be used on panel images to get a list of the radiance of each band (in the native RedEdge band order).
import numpy as np
-flight = df.altitude>cutoff_altitude
-ground = ~flight
-ground_idx = np.arange(len(ground))[ground]
-flight_idx = np.arange(len(ground))[flight]
-ground_captures = np.array(imgset.captures)[ground_idx]
-flight_captures = np.array(imgset.captures)[flight_idx]
-
-panel_radiances = []
-dls_irradiances = []
-panel_timestamps = []
-for cap in ground_captures:
- if cap.panels_in_all_expected_images():
- panel_timestamps.append(cap.utc_time())
- panel_radiances.append(cap.panel_radiance())
- dls_irradiances.append(cap.dls_irradiance())
-
-dls_irradiances = np.asarray(dls_irradiances)
-panel_radiances = np.asarray(panel_radiances)
-
-###
-panel_reflectance_by_band = [0.49, 0.49, 0.49, 0.49, 0.49] #RedEdge band_index order
-panel_irradiance = ground_captures[0].panel_irradiance(panel_reflectance_by_band)
-plt.figure()
-plt.scatter(ground_captures[0].center_wavelengths(), panel_irradiance);
-plt.xscale('log')
-plt.ylabel("Irradiance (w/m^2/nm)");
-plt.xlabel("Wavelength (nm)");
-Once we have extracted the panel radiances from the ground images, we can convert them to irradiance as we did in tutorial #1 and plot the irradiances over time and extract radiance to reflectance conversion factors.
- -import matplotlib.pyplot as plt
-import math
-
-df[df.altitude>cutoff_altitude].plot(y=columns[8:13], figsize=(14,8))
-plt.ylabel("Irradiance (w/m^2/nm)");
-Copyright (c) 2017-2018 MicaSense, Inc. For licensing information see the project git repository
- -This notebook shows usage for the Image type. This type is useful for loading single image files, reading the embedded metadata, and performing radiometric calibrations.
- -from micasense.image import Image
-import os, glob
-%matplotlib inline
-
-imagePath = os.path.join('.','data','ALTUM')
-imageName = glob.glob(os.path.join(imagePath,'IMG_0021_1.tif'))[0]
-
-img = Image(imageName)
-img.plot_all();
-import numpy as np
-import matplotlib.pyplot as plt
-import cv2
-
-nbins = 1024
-vmin = 0
-vmax = 2**16
-bins = range(vmin,vmax, int(vmax/nbins))
-hist = cv2.calcHist([img.raw().ravel()],[0],None,[nbins],[vmin,vmax])
-plt.plot(bins,hist);
-plt.xlim(img.raw().min(),img.raw().max())
-plt.ylim(0,hist.max())
-plt.xlabel('Pixel Value')
-plt.ylabel('Frequency')
-plt.show()
-Copyright (c) 2017-2019 MicaSense, Inc. For licensing information see the project git repository
- -This series of tutorials will be a walk through on how to process RedEdge data from raw images through conversion to reflectance. In this first tutorial, we will cover the tools required to do this, get them installed, and verify that the installation works.
-Our tutorials are written using Python3. Python has great library support for image processing through libraries such as OpenCV, SciKit Image, and others. In this tutorial, we'll use python, OpenCV, numpy, and matplotlib, as well as the standalone exiftool and it's python wrapper to open and manipulate RedEdge images to transform raw digital number values into quantitative reflectance.
-This tutorial has been tested on Windows, MacOS, and Linux. It is likely to work on other platforms, especially unix-based platforms like OSX, but you will have to do the legwork to get the required software installed and working.
-The following softare and libraries are required for this tutorial:
-Below, we go through the options to download and install a full working python environment with these tools (and their dependencies). We're using the Anaconda or miniconda environments where possible to ease installation, but if you're already a python package management guru, you can use git to checkout this code repository and look at the micasense_conda_env.yml file for the dependencies you'll need in your virtual environment.
For linux (and Mac, to some extent) you can either install the libraries directly using pip or install miniconda or anaconda to create completely separate environments. We have had success installing miniconda locally -- it's a smaller install than anaconda and can be installed without using sudo and doesn't impact the system-installed python or python libraries. You will likely still need to use sudo to install
The following is what we had to do on a fresh Ubuntu 18.04 image to install the library. First we installed some system tools and libraries:
- -sudo apt install git
-sudo apt install libzbar0
-sudo apt install make
-
-
-Next we installed exiftool:
- -wget https://exiftool.org/Image-ExifTool-12.57.tar.gz
-tar -xvzf Image-ExifTool-12.57.tar.gz
-cd Image-ExifTool-12.57/
-perl Makefile.PL
-make test
-sudo make install
-
-
-Then we installed miniconda. Navigate to the miniconda download page and download the installer for your system and follow the installation instructions
-Once these tools are installed, you can check out this repository and create the micasense conda environment:
git clone https://github.com/micasense/imageprocessing.git
-cd imageprocessing
-conda env create -f micasense_conda_env.yml
-
-
-Finally, one way to verify our install by running the built in tests:
- -cd imageprocessing
-conda activate micasense
-pytest .
-
-
-Or, to start working with the notebooks (including running the test code below):
- -cd imageprocessing
-conda activate micasense
-jupyter notebook .
-
-
-
-When installing on Windows we rely on the Anaconda python environment to do most of the heavy lifting for us.
-Install Anaconda for your system by downloading the Python 3.7 version
-Download the exiftool windows package and unzip it to a permanent location such as c:\exiftool\. Now we need to tell the python code where to find exiftool (so we don't have to set it up in every script we write), and we do that by adding the path to exiftool as an environment variable.
exiftoolpath with a value of the full path to exiftool. For example, c:\exiftool\exiftool.exePath and click Edit Environment Variables for Your AccountNewexiftoolpathc:\exiftool\exiftool.exeOpen an Anaconda console from the start menu as an administrator by clicking Start->Anaconda, right-click Anaconda Console, choose Run as Administrator. Execute the following commands in the anaconda console:
cd to the directory you git cloned this repository toconda env create -f micasense_conda_env.ymlactivate micasense to activate the environment configuredactivate micasenseThe following steps to get going on MacOS worked for us.
-First we installed git by installing the xcode developer tools, or you can follow the instructions at the git site.
Next, we downloaded and installed exiftool using the MacOS installer.
-Third, we installed Homebrew and used it to install zbar:
/usr/bin/ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"
-brew install zbar
-
-
-Then we installed miniconda. If you're comfortable on the command line, navigate to the miniconda download page and download the installer for your system. Open an iTerm and follow the installation instructions.
-If instead you're more comfortable with graphical installers, the Anaconda version for Python 3.7 may be right for you.
-Once these tools are installed, you can check out this repository and create the micasense conda environment by opening an iTerm and running the following commands:
git clone https://github.com/micasense/imageprocessing.git
-cd imageprocessing
-conda env create -f micasense_conda_env.yml
-
-
-This will take a while (5-10 minutes isn't uncommon). Once it's done, one way to verify our install by running the built-in tests:
- -conda activate micasense
-pytest .
-
-
-Or, to start working with the notebooks (including running the test code below):
- -cd imageprocessing
-conda activate micasense
-jupyter notebook .
-
-
-
-Anaconda Prompt from the Start menu and type activate micasensecd to the imageprocessing checkout directoryjupyter notebook .The following python snippet can be run from a jupyter notebook, inside iPython, or by saving to a script and running from the command line. If you're on windows, make sure you have set the location of exiftool in the exiftoolpath environment variable. If this script succeeds, your system is ready to go! If not, check the installation documentation for the module import that is having issues.
import cv2 #openCV
-import exiftool
-import os, glob
-import numpy as np
-import pyzbar.pyzbar as pyzbar
-import matplotlib.pyplot as plt
-import mapboxgl
-
-print()
-print("Successfully imported all required libraries.")
-print()
-
-if os.name == 'nt':
- if os.environ.get('exiftoolpath') is None:
- print("Set the `exiftoolpath` environment variable as described above")
- else:
- if not os.path.isfile(os.environ.get('exiftoolpath')):
- print("The provided exiftoolpath isn't a file, check the settings")
-
-try:
- with exiftool.ExifTool(os.environ.get('exiftoolpath')) as exift:
- print('Successfully executed exiftool.')
-except Exception as e:
- print("Exiftool isn't working. Double check that you've followed the instructions above.")
- print("The execption text below may help to find the source of the problem:")
- print()
- print(e)
--Successfully imported all required libraries. - -Successfully executed exiftool. --
The above code checks for the proper libraries to be installed and verifies it can execute exiftool. This code opens an example image, reads the metadata, and then uses the pyzbar library to find a MicaSense panel in the image.
from micasense.image import Image
-imagePath = os.path.join('.','data','REDEDGE-MX')
-imageName = glob.glob(os.path.join(imagePath,'IMG_0001_1.tif'))[0]
-
-img = Image(imageName)
-img.plot_raw(figsize=(8.73,8.73));
-
-from micasense.panel import Panel
-panel = Panel(img)
-if not panel.panel_detected():
- raise IOError("Panel Not Detected! Check your installation of pyzbar")
-else:
- panel.plot(figsize=(8,8));
-
-print('Success! Now you are ready for Part 1 of the tutorial.')
-Success! Now you are ready for Part 1 of the tutorial. --
Copyright (c) 2017-2019 MicaSense, Inc. For licensing information see the project git repository
- -This tutorial assumes you have gone through the basic setup here and your system is set up and ready to go.
-In this tutorial, we will walk through how to convert RedEdge data from raw images to radiace and then to reflectance. We will cover the tools required to do this, and walk through some of the basic image processing and radiometric conversions.
-RedEdge 16-bit images can be read with pyplot directly into numpy arrays using the pyplot imread function or the matplotlib imread function, and then we can display the image inline using the imshow function of matplotlib.
import cv2
-import matplotlib.pyplot as plt
-import numpy as np
-import os,glob
-import math
-%matplotlib inline
-
-imagePath = os.path.join('.','data','REDEDGE-MX')
-imageName = os.path.join(imagePath,'IMG_0001_4.tif')
-
-# Read raw image DN values
-# reads 16 bit tif - converts 12 bit to 16 bit if necessary
-imageRaw=plt.imread(imageName)
-
-# Display the image
-fig, ax = plt.subplots(figsize=(8,6))
-ax.imshow(imageRaw, cmap='gray')
-plt.show()
-For many of the steps in the tutorial, we will use code from the MicaSense utilities module. The code is in the micasense directory and can be imported via normal python import commands using the syntax import micasense or import micasense.submodule as short_name for use in this and other scripts. While we will not cover all of the utility functions in this tutorial, they are available for reference and some will be used and discussed in future tutorials.
We will use start by using a plotting function in micasense.plotutils that adds a colorbar to the display, so that we can more easily see changes in the values in the images and also see the range of the image values after various conversions. This function also colorizes the grayscale images, so that changes can more easily be seen. Depending on your viewing style, you may prefer a different color map and you can also select that colormap here or browsing the colormaps on the matplotlib site.
import micasense.plotutils as plotutils
-
-# Optional: pick a color map that fits your viewing style
-# one of 'gray, viridis, plasma, inferno, magma, nipy_spectral'
-plotutils.colormap('viridis');
-
-fig = plotutils.plotwithcolorbar(imageRaw, title='Raw image values with colorbar')
-<Figure size 432x288 with 0 Axes>-
In order to perform various processing on the images, we need to read the metadata of each image. For this we use ExifTool. We can read standard image capture metadata such as location, UTC time, imager exposure and gain, but also RedEdge specific metadata which can make processing workflows easier.
-For example, each image contains a unique capture identifier. Capture identifiers are shared between all 5 images captured by RedEdge at the same moment, and can be used to unambiguously group images in post processing, regardless of how the images are named or stored on disk. Each image also contains a flight identifer which is the same for all images taken during a single power cycle of the camera. This can be used in post-processing workflows to group images and in many cases, more easily identify when the vehicle took off and landed.
- -import micasense.metadata as metadata
-exiftoolPath = None
-if os.name == 'nt':
- exiftoolPath = os.environ.get('exiftoolpath')
-# get image metadata
-meta = metadata.Metadata(imageName, exiftoolPath=exiftoolPath)
-cameraMake = meta.get_item('EXIF:Make')
-cameraModel = meta.get_item('EXIF:Model')
-firmwareVersion = meta.get_item('EXIF:Software')
-bandName = meta.get_item('XMP:BandName')
-print('{0} {1} firmware version: {2}'.format(cameraMake,
- cameraModel,
- firmwareVersion))
-print('Exposure Time: {0} seconds'.format(meta.get_item('EXIF:ExposureTime')))
-print('Imager Gain: {0}'.format(meta.get_item('EXIF:ISOSpeed')/100.0))
-print('Size: {0}x{1} pixels'.format(meta.get_item('EXIF:ImageWidth'),meta.get_item('EXIF:ImageHeight')))
-print('Band Name: {0}'.format(bandName))
-print('Center Wavelength: {0} nm'.format(meta.get_item('XMP:CentralWavelength')))
-print('Bandwidth: {0} nm'.format(meta.get_item('XMP:WavelengthFWHM')))
-print('Capture ID: {0}'.format(meta.get_item('XMP:CaptureId')))
-print('Flight ID: {0}'.format(meta.get_item('XMP:FlightId')))
-print('Focal Length: {0}'.format(meta.get_item('XMP:FocalLength')))
-MicaSense RedEdge-M firmware version: v7.5.0-beta6 -Exposure Time: 0.00036 seconds -Imager Gain: 1.0 -Size: 1280x960 pixels -Band Name: NIR -Center Wavelength: 842 nm -Bandwidth: 57 nm -Capture ID: Rb0pibHa08uHJwrTjf8Y -Flight ID: ePdxBBmSitgkTdpwZiM9 -Focal Length: None --
Ultimately most RedEdge users want to calibrate raw images from the camera into reflectance maps. This can be done using off-the-shelf software from third parties, but you are here because there is no fun in that! Along with this tutorial we have included some helper utilities that will handle much of this conversion for you, but here we will walk through a few of those functions to discuss what is happening inside.
-Any RedEdge workflow must include these common steps.
-All of these are handled by the micasense.utils.raw_image_to_radiance(metadata, raw_image) function. Let us take a look at that fuction in more detail.
First, we get the darkPixel values. These values come from optically-covered pixels on the imager which are exposed at the same time as the image pixels. They measure the small amount of random charge generation in each pixel, independent of incoming light, which is common to all semiconductor imaging devices.
-blackLevel = np.array(meta.get_item('Exif.BlackLevel'))
- darkLevel = blackLevel.mean()
-Now, we get the imager-specific calibrations.
-a1, a2, a3 = meta.get_item('XMP:RadiometricCalibration')
-We get the parameters of the optical chain (vignette) effects and create a vignette map. This map will be multiplied by the black-level corrected image values to reverse the darkening seen at the image corners. See the vignette_map function for the details of the vignette parameters and their use.
V, x, y = vignette_map(meta, xDim, yDim)
-Now we can calculate the imager-specfic radiometric correction function, which help to account for the radiometric inaccuracies of the CMOS imager pixels.
-# row gradient correction
- R = 1.0 / (1.0 + a2 * y / exposureTime - a3 * y)
-Finally, we apply these functions to the raw image to result in a corrected image
-# subtract the dark level and adjust for vignette and row gradient
- L = V * R * (imageRaw - darkLevel)
-Next, we get the exposure and gain settings (gain is represented in the photographic parameter ISO, with a base ISO of 100, so we divide the result to get a numeric gain).
-exposureTime = float(meta.get_item('EXIF:ExposureTime'))
- gain = float(meta.get_item('EXIF:ISOSpeed'))/100.0
-Now that we have a corrected image, we can apply a conversion from calibrated digital number values to radiance units (W/m^2/nm/sr). Note that in this conversion, we need to normalize by the image bitdepth (2^16 for 16 bit images, 2^12 for 12-bit images), because the calibration coefficients are scaled to work with normalized input values.
-# apply the radiometric calibration -
- # scale by the gain-exposure product and multiply with the radiometric calibration coefficient
- bitsPerPixel = meta.get_item('EXIF:BitsPerSample')
- dnMax = float(2**bitsPerPixel)
- radianceImage = L.astype(float)/(gain * exposureTime)*a1/dnMax
-For convenience, we have written the raw_image_to_radiance function to return the intermediate compensation images as well, so we can visualize them for the tutorial. These intermediate results are not required in most implementations and can be ommitted, if performance is a concern.
import micasense.utils as msutils
-radianceImage, L, V, R = msutils.raw_image_to_radiance(meta, imageRaw)
-plotutils.plotwithcolorbar(V,'Vignette Factor');
-plotutils.plotwithcolorbar(R,'Row Gradient Factor');
-plotutils.plotwithcolorbar(V*R,'Combined Corrections');
-plotutils.plotwithcolorbar(L,'Vignette and row gradient corrected raw values');
-plotutils.plotwithcolorbar(radianceImage,'All factors applied and scaled to radiance');
-Now that we have a flat and calibrated radiance image, we can convert into reflectance. To do this, we will use the radiance values of the panel image of known reflectance to determine a scale factor between radiance and reflectance.
-In this case, we have our MicaSense calibrated reflectance panel and its known reflectance of 49% in the band of interest. We will extract the area of the image containing the lambertian panel, determine it's radiance to reflectance scale factor, and then scale the whole image by that factor to get a reflectance image.
- -markedImg = radianceImage.copy()
-ulx = 610 # upper left column (x coordinate) of panel area
-uly = 610 # upper left row (y coordinate) of panel area
-lrx = 770 # lower right column (x coordinate) of panel area
-lry = 760 # lower right row (y coordinate) of panel area
-cv2.rectangle(markedImg,(ulx,uly),(lrx,lry),(0,255,0),3)
-
-# Our panel calibration by band (from MicaSense for our specific panel)
-panelCalibration = {
- "Blue": 0.49,
- "Green": 0.49,
- "Red": 0.49,
- "Red edge": 0.49,
- "NIR": 0.49
-}
-
-# Select panel region from radiance image
-panelRegion = radianceImage[uly:lry, ulx:lrx]
-plotutils.plotwithcolorbar(markedImg, 'Panel region in radiance image')
-meanRadiance = panelRegion.mean()
-print('Mean Radiance in panel region: {:1.3f} W/m^2/nm/sr'.format(meanRadiance))
-panelReflectance = panelCalibration[bandName]
-radianceToReflectance = panelReflectance / meanRadiance
-print('Radiance to reflectance conversion factor: {:1.3f}'.format(radianceToReflectance))
-
-reflectanceImage = radianceImage * radianceToReflectance
-plotutils.plotwithcolorbar(reflectanceImage, 'Converted Reflectane Image');
-Mean Radiance in panel region: 0.105 W/m^2/nm/sr -Radiance to reflectance conversion factor: 4.649 --
In some cases we might notice that some parts of a converted reflectance image show reflectances above 1.0, or 100%, and wonder how this is possible. In fact, reflectances higher than 100% are normal in specific cases of specular reflections. The panel area is a special material that reflects incident light equally well in all directions; however, some of the objects in a scene (especially man-made objects) instead reflect most incident light in one direction, more like a mirror. An example is the reflection of the sun off of the smooth surface of a car or the reflection off of a body of water.
-Now we will extract the same region and verify the reflectance in that region is what we expect. In the process, we will blur and visualize the extracted area to look for any trends. The area should have a very consistent reflectance. If a gradient or a high standard deviation (>3% absolute reflectance) is noticed across the panel area it is possible that the panel was captured under inconsistent lighting conditions (e.g. next to a wall or vehicle) or it was captured too close to the edge of the image where the optical calibration is the least accurate.
- -panelRegionRaw = imageRaw[uly:lry, ulx:lrx]
-panelRegionRefl = reflectanceImage[uly:lry, ulx:lrx]
-panelRegionReflBlur = cv2.GaussianBlur(panelRegionRefl,(55,55),5)
-plotutils.plotwithcolorbar(panelRegionReflBlur, 'Smoothed panel region in reflectance image')
-print('Min Reflectance in panel region: {:1.2f}'.format(panelRegionRefl.min()))
-print('Max Reflectance in panel region: {:1.2f}'.format(panelRegionRefl.max()))
-print('Mean Reflectance in panel region: {:1.2f}'.format(panelRegionRefl.mean()))
-print('Standard deviation in region: {:1.4f}'.format(panelRegionRefl.std()))
-Min Reflectance in panel region: 0.43 -Max Reflectance in panel region: 0.59 -Mean Reflectance in panel region: 0.49 -Standard deviation in region: 0.0152 --
In this case the panel is less uniform than we would like, but we can also notice that the full color scale is only 4% of absolute reflectance. Likewise, it's well below the standard deviation threshold we have set of 3% absolute reflectance. This panel has seen over two years of hard field use, so it may be time for it to retire.
-Reasons for a high standard deviation across a panel can include panel contamination or inconsistent lighting across the panel due to environmental conditions. Based on the context of the image, it is also clear that the user is taking the panel image facing the sun, which can cast reflected light from the operator's clothing on the panel and contaminate results. For this reason it is always best to capture panel images in an open area and with the operator's back to the sun.
-Finally, we need to remove lens distortion effects from images for some processing workflows, such as band-to-band image alignment. Generally for photogrammetry processes on raw (or radiance/reflectance) images, this step is not required, as the photogrammetry process will optimize a lens distortion model as part of it's bulk bundle adjustment. RedEdge has very low distortion lenses, so the changes to images in this step tend to be very small and noticeable only in pixels on the border of the image.
- -# correct for lens distortions to make straight lines straight
-undistortedReflectance = msutils.correct_lens_distortion(meta, reflectanceImage)
-plotutils.plotwithcolorbar(undistortedReflectance, 'Undistorted reflectance image');
-Now that we can convert from raw RedEge images to reflectance, we will use these methods to convert an image taken during the same campaign to a reflectance image, and look at a few interesting areas of the image to validate our conversion.
- -flightImageName = os.path.join(imagePath,'IMG_0001_4.tif')
-flightImageRaw=plt.imread(flightImageName)
-plotutils.plotwithcolorbar(flightImageRaw, 'Raw Image')
-
-flightRadianceImage, _, _, _ = msutils.raw_image_to_radiance(meta, flightImageRaw)
-flightReflectanceImage = flightRadianceImage * radianceToReflectance
-flightUndistortedReflectance = msutils.correct_lens_distortion(meta, flightReflectanceImage)
-plotutils.plotwithcolorbar(flightUndistortedReflectance, 'Reflectance converted and undistorted image');
-In this tutorial we have found that we can read MicaSense RedEdge images and their metadata, and use python and OpenCV to convert those images to radiance and then to reflectance using the standard scientific field method of imaging a lambertian reflector. We have corrected for both the electro-optical effects of the sensor and optical chain, as well as the incident light at the time of capture.
-In future tutorials, we will introduce the Downwelling Light Sensor (DLS) information into the calibration process in order to account for changing irradiance over time (e.g. such as clouds). However, since the panel method is straightforward and repeatable under constant illumination conditions, and is the standard scientific calibration method of surface reflectance, this process is useful and sufficient for many calibration needs.
-Looking for more? Try the second tutorial here.
- -Copyright (c) 2017-2019 MicaSense, Inc. For licensing information see the project git repository
- -This tutorial assumes you have gone through the basic setup and builds on the basic radiance, irradiance, and reflectance concepts and code covered in the first tutorial.
-In this tutorial, we will cover usage of the MicaSense python library to access images and groups of images. Most of the processing details are hidden away in the library, but the library code is open and available in the git repository.
-In the first tutorial, we introduced micasense.utils which provided some helper functions for single image manipulation, and micasense.plotutils which provided some plotting helpers.
For this second tutorial, we are going to introduce the usage of the included micasense libraries for opening, converting, and displaying images. This will allow us to discuss and visualize results at a high level, while the underlying source code is available for those interested in the implementation details. In some cases, the libraries themselves may be enough to implement a custom workflow without the need to re-implement or translate the code to another system or language.
-The library code provides some basic classes to manage image data. At the highest level is the ImageSet, which is able to load a list of files or recursively search a whole directory into data structures which are easy to access and manipulate. ImageSets are made up of Captures, which hold the set of (usually 5) images as they are simultaneously gathered by the RedEdge camera. Within Captures are Images, which hold a single image file and allow easy access to the image metadata. The Image class also provides the ability to extract metadata from individual images and to convert individual images in similar ways to those described in the first tutorial.
For the rest of this article, we will look at each of the objects available starting with the single Image object, and work our way up to the whole ImageSet. Each section in this article is standalone, and can be copied into another workbook or edited in place to explore more of the functions associated with that object.
An image is the lowest level object. It represents the data in a single tiff file as taken by the camera. Image objects expose a set of data retrieval methods which provide access to raw, radiance, and reflectance corrected images, and to undistort any of those images. Note that when retrieving image data from an Image object, the data is stored internally in the object, increasing the object's memory footprint. If operating on a large number of images, it may be necessary to release this data memory after each image is processed to limit the program memory footprint. This can be done by calling the Image.clear_image_data() method.
import os
-import micasense.image as image
-%matplotlib inline
-
-image_path = os.path.join('.','data','ALTUM','IMG_0000_1.tif')
-img = image.Image(image_path)
-img.plot_raw();
-Image Metadata¶Metadata for each image is available in the Image.meta parameter. This object is a micasense.Metadata object and can be accessed directly for image specific metadata extraction. Below, we print the same metadata values as we did in Tutorial #1, but using direct access to the Metadata object parameters.
A notebook for experimenting with the Image class can be found here.
print('{0} {1} firmware version: {2}'.format(img.meta.camera_make(),
- img.meta.camera_model(),
- img.meta.firmware_version()))
-print('Exposure Time: {0} seconds'.format(img.meta.exposure()))
-print('Imager Gain: {0}'.format(img.meta.gain()))
-print('Size: {0}x{1} pixels'.format(img.meta.image_size()[0],
- img.meta.image_size()[1]))
-print('Band Name: {0}'.format(img.meta.band_name()))
-print('Center Wavelength: {0} nm'.format(img.meta.center_wavelength()))
-print('Bandwidth: {0} nm'.format(img.meta.bandwidth()))
-print('Capture ID: {0}'.format(img.meta.capture_id()))
-print('Flight ID: {0}'.format(img.meta.flight_id()))
-MicaSense Altum firmware version: v2.5.0-beta6 -Exposure Time: 0.000423757 seconds -Imager Gain: 1.0 -Size: 2064x1544 pixels -Band Name: Blue -Center Wavelength: 475 nm -Bandwidth: 20 nm -Capture ID: CoDN4adNch7KFnclyi5t -Flight ID: jfCG6qSK0t7CDUvnJiOO --
The Capture class is a container for Images which allows access to metadata common to the group of images. The internal Image objects are accessible via the capture.images properties, and images in this list are kept sorted by the band property. Data which is different for each image can be accessed through composite methods, such as the capture.dls_irradiance() method, which returns a list of irradiances in band order.
import os, glob
-import micasense.capture as capture
-
-images_path = os.path.join('.','data','ALTUM')
-image_names = glob.glob(os.path.join(images_path,'IMG_0000_*.tif'))
-cap = capture.Capture.from_filelist(image_names)
-cap.plot_radiance();
-Capture metadata¶Metadata which is common to all captures can be accessed via methods on the Capture object. Metadata which varies between the images of the capture, such as DLS information, is available as lists accessed from the capture object.
Below we plot the raw and tilt compensated DLS irradiance by center wavelength and by band name.
- -import matplotlib.pyplot as plt
-
-print(cap.band_names())
-fig = plt.figure(figsize=(14,6))
-plt.subplot(1,2,1)
-plt.scatter(cap.center_wavelengths(), cap.dls_irradiance())
-plt.ylabel('Irradiance $(W/m^2/nm)$')
-plt.xlabel('Center Wavelength (nm)')
-plt.subplot(1,2,2)
-plt.scatter(cap.band_names(), [img.meta.exposure() for img in cap.images])
-plt.xlabel('Band Names')
-plt.ylim([0,2.5e-3])
-plt.ylabel('Exposure Time (s)')
-plt.show()
-['Blue', 'Green', 'Red', 'NIR', 'Red edge', 'LWIR'] --
A notebook for experimenting with the Capture class can be found here.
The Panel class is a helper class which can automatically extract panel information from MicaSense calibrated reflectance panels by finding the QR code within an image and using the QR Code location and orientation information to find the lambertian panel area. The class then allows extraction of statistics from the panel area such as mean raw values, mean radiance, standard deviation, and the number of saturated pixels in the panel region. The panel object can be included standalone, or used within the context of a Capture object.
import os, glob
-import micasense.image as image
-import micasense.panel as panel
-
-image_path = os.path.join('.','data','ALTUM','IMG_0000_1.tif')
-img = image.Image(image_path)
-
-pnl = panel.Panel(img)
-print("Panel found: {}".format(pnl.panel_detected()))
-print("Panel serial: {}".format(pnl.serial))
-print("QR Code Corners:\n{}".format(pnl.qr_corners()))
-mean, std, count, saturated_count = pnl.raw()
-print("Panel mean raw pixel value: {}".format(mean))
-print("Panel raw pixel standard deviation: {}".format(std))
-print("Panel region pixel count: {}".format(count))
-print("Panel region saturated pixel count: {}".format(count))
-
-pnl.plot();
-Panel found: True -Panel serial: RP06-2051037-OB -QR Code Corners: -None -Panel mean raw pixel value: 32852.294381134314 -Panel raw pixel standard deviation: 396.866996832778 -Panel region pixel count: 15181 -Panel region saturated pixel count: 15181 --
A notebook for experimenting with the Panel class can be found here
An ImageSet contains a group of Captures. The captures can be loaded from image object, from a list of files, or by recursively searching a directory for images.
Loading an ImageSet can be a time consuming process. It uses python multithreading under the hood to maximize cpu usage on multi-core machines.
from ipywidgets import FloatProgress
-from IPython.display import display
-f = FloatProgress(min=0, max=1)
-display(f)
-def update_f(val):
- f.value=val
-
-import micasense.imageset as imageset
-import os
-images_dir = os.path.join('.','data','ALTUM')
-
-imgset = imageset.ImageSet.from_directory(images_dir, progress_callback=update_f)
-
-for cap in imgset.captures:
- print ("Opened Capture {} with bands {}".format(cap.uuid,[str(band) for band in cap.band_names()]))
-FloatProgress(value=0.0, max=1.0)-
Opened Capture CoDN4adNch7KFnclyi5t with bands ['Blue', 'Green', 'Red', 'NIR', 'Red edge', 'LWIR'] -Opened Capture E9PO6BRS2ua7AnqabMMH with bands ['Blue', 'Green', 'Red', 'NIR', 'Red edge', 'LWIR'] --
A large group of images captured over a central California orchard are available for download here.
-With this set extracted to a working folder, the extended ImageSet example notebook provides more usages of ImageSet data.
-In this tutorial, we have introduced the MicaSense library and provided some examples of opening Images, Captures, and ImageSets, as well as detecting and extracting panel information from images.
-The next tutorial covers basic usage of DLS information, and is available here
- -Copyright (c) 2017-2019 MicaSense, Inc. For licensing information see the project git repository
- -This tutorial will walk through how to convert RedEdge data from raw images to radiace and then use the DLS information in the image metadata to convert that to reflectance.
- -This tutorial was originally written as a guide to basic processing of DLS1 data for RedEdge-3 and RedeEdge-M. Since October 2018, Altum and RedEdge-MX have been shipped standard with DLS2. DLS1 is a small, red square with a single diffuser on top. DLS2 is a larger, thinner black rectangle with two small diffusers on top and 8 small diffusers on different up-looking faces.
-For DLS2 data, we recommend using the Capture class to access the Capture.dls_irradiance() method. This will provide a compensated horizontal irradiance useful for radiometrically correcting imagery. We leave the below intact for legacy users and as a tutorial on remote sensing calibrations and terms.
For RedEdge and Altum with up-to-date firmware, the EXIF:HorizontalIrradiance tag and the appropriate scale factor will provide the necessary information for reflectance compensation. If necessary, this can be further field-calibrated by using a pane image to tie the DLS2 and camera calibrations together in-situ.
In Tutorial 1 we covered conversion of the image radiance to reflecance using the average radiance from an area of a specially made lambertian reflectance panel. As we did not get very far into remote sensing theory, we will cover that some here as we move on to the more complex problem of using the data from the RedEdge Downwelling Light Sensor (DLS). To get started, we will review some of the fundamental measurements in remote sensing.
-The DLS measures the energy coming from the sky, while the camera measures the energy coming from the fields below. The property that we are generally looking to estimate in agricultural remote sensing is the surface reflectance of the plants and ground below the sensor.
-The MicaSense RedEdge measures spectral radiance, a measure of the power incident on the CMOS image sensor at each pixel, and has units power per unit area per unit bandwidth per unit solid angle (W/m^2/nm/sr). We will call the radiance image L.
-The MicaSense Downwelling Light Sensor (DLS) measures incoming spectral irradiance on the top surface of the device for the same 5 spectral bands as the camera. Irradiance has the units of power per unit area per unit bandwidth, or Watts per meter squared per nanometer (W/m^2/nm). We will call this downwelling quantity I.
-If we look closely at the units, they differ by the last term - the unit over solid angle. This is because each pixel of the camera sensor only sees a specific solid angle, while the dls detector has a specific area.
-In order to use downwelling irradiance and upwelling radiance to determine reflectance, we simply need to divide pi times the irradiance by the radiance. Our units cancel out, and we are left with the unitless quantity of reflectance. However, there is one caveat. In most systems, the DLS is not held level to the ground at all times, so the amount of light that shines on it's surface changes as the aircraft moves about. In oder for that simple math to work, it's necessary that the irradiance and radiance be compensated for a very important property: the illumination angle.
-Both the radiance and irradiance are sensitive to changes in the direction of the lighting. Both quantites decrease with the cosine of the angle from the surface normal, eventually reaching zero when the light source is parallel to the surface. So when measuring irradiance, this effect needs to be corrected in order to estimate the actual irradiance.
-To complicate matters further, the light reaching the DLS and the ground is not coming from a single source, even under sunny conditions. The direct sunlight provides what we will call the direct source, resulting in a direct irradiance I_direct. The rest of the (usually blue) sky provides what we term the indirect, or diffuse, irradiance which we will call I_diffuse. Under even the sunniest conditions with the summer sun directly overhead, only about 85% of the light is direct, while the other 15% is diffuse, which we will express as the ratio 6:1. Under thin cloud conditions, this ratio drops to close 1:1, meaning half of the measured irradiance is coming from parts of the sky other than the direction of the sun. The ratio goes even lower under fully cloudy skies, becoming fully diffuse.
As we ultimately want to estimate the time-varying irradiance, we will need to estimate both the direct and diffuse irradiance terms as well as correct for the cosine effect of the DLS sensor itself.
-There's one final piece of theory required that pertains to irradiance sensing, and that is the Fresnel Effect. The Fresnel Effect describes a physical property of the diffuser that gathers light on top of the DLS. This effect results from the fact that much of the light striking the diffuser is reflected, and the amount reflected changes with the angle of the incident light. So, instead of our diffuser behaving as a perfect cosine detector, we need to also compensate for the material imperfections.
-An onboard orientation sensor provides a fairly rough version of the earth-fixed orientation of the DLS. The quality of this orientation varies based on the quality of the sensor installation and calibration, but in even the best cases it has a few degrees of error. If the DLS and the camera are mounted fixed to one another, the photogrammetric workflow can be used to augment the orientation of the irradiance sensor. Our examples here will use the sensor measurements directly, but significant improvements can be made using the very precise angles photorammetry can provide.
-As we've discussed in the first tutorial, the method of using a calibrated reflectance panel is the scientific standard for remote sensing measurements. Spectroradiometers which cost tens of thousands of dollars are calibrated in this method in the field. Panels images may also be used both before and after each flight to track slowly changing solar irradiance. In many cases of direct sunlight with no clouds, panel images can provide the best results when taken before and after each flight. As the diffuse lighting increases, the utility of the DLS increases as well.
-In cases of more varied irradiance, the DLS can be used for irradiance estimation, using only the internal pose estimation. The DLS irradiance can be be used in conjunction with panels to ensure that the DLS and panel measurements match at the time of panel capture. In order for this to be effective, it's important that the DLS and panel are exposed to the same irradiance. This is best accomplished by ensuring panel images are taken with the operator's back to the sun and the operator shadow next to the panel. The aircraft and operator should block as little of the sky as possible.
-Finally, in fixed mount systems, the photogrammetric angles can be used, along with the offset between the camera and DLS sensors, to most accurately estimate irradiance. This can again be coupled with one or more panel measurements for best radiometric results.
- -Below, we will use the micasense.Capture class to load a set of images and use tilt-compensated DLS irradiance data to roughly translate those images to reflectance. We will then compare these to the known reflectances of a lambertian panel, and develop a site-specific compensation factor to tie together the DLS irradiance and camera radiance measurements through the panel.
Below, we show how to access DLS data and compute the minimum required corrections to the DLS data for use. The code in this cell is present in the micasense.dls library, and is used by the Image and Capture class to expose corrected DLS data. If your intent is to simply use this code to correct your data, you can skip this section and access pose-corrected DLS data through the capture.dls_irradiance() method.
For the purpose of this tutorial, we assume a constant irradiance (over time) and a ratio of direct to diffuse irradiance of 6:1. These assumptions only hold approximately, and only for clear sky. Cloud cover can lead to dramatically different results.
-import numpy as np
-import micasense.dls as dls
-
-import os, glob
-import micasense.capture as capture
-
-images_path = os.path.join('.','data','REDEDGE-P')
-image_names = glob.glob(os.path.join(images_path,'IMG_0000_*.tif'))
-cap = capture.Capture.from_filelist(image_names)
-# Define DLS sensor orientation vector relative to dls pose frame
-dls_orientation_vector = np.array([0,0,-1])
-# compute sun orientation and sun-sensor angles
-(
- sun_vector_ned, # Solar vector in North-East-Down coordinates
- sensor_vector_ned, # DLS vector in North-East-Down coordinates
- sun_sensor_angle, # Angle between DLS vector and sun vector
- solar_elevation, # Elevation of the sun above the horizon
- solar_azimuth, # Azimuth (heading) of the sun
-) = dls.compute_sun_angle(cap.location(),
- cap.dls_pose(),
- cap.utc_time(),
- dls_orientation_vector)
-# Since the diffuser reflects more light at shallow angles than at steep angles,
-# we compute a correction for this
-fresnel_correction = dls.fresnel(sun_sensor_angle)
-
-# Now we can correct the raw DLS readings and compute the irradiance on level ground
-dls_irradiances = []
-center_wavelengths = []
-for img in cap.images:
- dir_dif_ratio = 6.0
- percent_diffuse = 1.0/dir_dif_ratio
- # measured Irradiance / fresnelCorrection
- sensor_irradiance = img.spectral_irradiance / fresnel_correction
- untilted_direct_irr = sensor_irradiance / (percent_diffuse + np.cos(sun_sensor_angle))
- # compute irradiance on the ground using the solar altitude angle
- dls_irr = untilted_direct_irr * (percent_diffuse + np.sin(solar_elevation))
- dls_irradiances.append(dls_irr)
- center_wavelengths.append(img.center_wavelength)
-
-import matplotlib.pyplot as plt
-plt.scatter(center_wavelengths,dls_irradiances)
-plt.xlabel('Wavelength (nm)')
-plt.ylabel('Irradiance ($W/m^2/nm$)')
-plt.show();
-
-cap.plot_undistorted_reflectance(dls_irradiances)
-Finally, if we want to use the panel as the ultimate reference tying the panel and the DLS readings together with a known lambertian reflector, we can calculate a DLS factor which we can then apply to DLS derived irradiances over the whole flight.
- -import math
-
-panel_reflectance_by_band = [0.49] * len(cap.eo_band_names()) # This panel has an average of 49% reflectance for each EO band
-
-panel_radiances = np.array(cap.panel_radiance())
-irr_from_panel = math.pi * panel_radiances / panel_reflectance_by_band
-dls_correction = irr_from_panel/dls_irradiances
-cap.plot_undistorted_reflectance(dls_irradiances*dls_correction)
-
-plt.scatter(cap.center_wavelengths(), cap.panel_reflectance())
-plt.title("Panel Reflectances")
-plt.xlabel("Wavelength (nm)")
-plt.ylabel("Reflectance")
-plt.show()
-While the need to compensate the DLS for movement increases the complexity of the problem, it also provides us more information than might be available for other diffuser types. This is because as the DLS moves about on the aircraft, it provides many measurements per second that should all have consistent effects assuming the light isn't changing over those very fast measurements. This movement can then be used to develop an irradiance series which uses all of the available information from the on board sensors. A future tutorial will provide an overview of some methods for improving irradiance sensing further.
- -Copyright (c) 2017-2019 MicaSense, Inc. For licensing information see the project git repository
- -This notebook shows usage for the Panel class. This type is useful for detecting MicaSense calibration panels and extracting information about the lambertian panel surface.
- -%load_ext autoreload
-%autoreload 2
-import os, glob
-from micasense.image import Image
-from micasense.panel import Panel
-%matplotlib inline
-
-imagePath = os.path.join('.','data','ALTUM')
-imageName = glob.glob(os.path.join(imagePath,'IMG_0000_4.tif'))[0]
-
-img = Image(imageName)
-panel = Panel(img)
-
-if not panel.panel_detected():
- raise IOError("Panel Not Detected!")
-
-print("Detected panel serial: {}".format(panel.serial))
-mean, std, num, sat_count = panel.raw()
-print("Extracted Panel Statistics:")
-print("Mean: {}".format(mean))
-print("Standard Deviation: {}".format(std))
-print("Panel Pixel Count: {}".format(num))
-print("Saturated Pixel Count: {}".format(sat_count))
-
-panel.plot();
-Detected panel serial: RP06-2051037-OB -Extracted Panel Statistics: -Mean: 31758.02988274707 -Standard Deviation: 386.3622284454486 -Panel Pixel Count: 14925 -Saturated Pixel Count: 0 --
# imagePath = os.path.join('.','data','ALTUM1SET','000')
-# imageName = glob.glob(os.path.join(imagePath,'IMG_0000_1.tif'))[0]
-
-imagePath = os.path.join('.','data','ALTUM')
-imageName = glob.glob(os.path.join(imagePath,'IMG_0000_4.tif'))[0]
-
-img = Image(imageName)
-if img.auto_calibration_image:
- print("Found automatic calibration image")
-panel = Panel(img)
-
-if not panel.panel_detected():
- raise IOError("Panel Not Detected!")
-
-print("Detected panel serial: {}".format(panel.serial))
-mean, std, num, sat_count = panel.raw()
-print("Extracted Panel Statistics:")
-print("Mean: {}".format(mean))
-print("Standard Deviation: {}".format(std))
-print("Panel Pixel Count: {}".format(num))
-print("Saturated Pixel Count: {}".format(sat_count))
-
-panel.plot();
-Found automatic calibration image -Detected panel serial: RP06-2051037-OB -Extracted Panel Statistics: -Mean: 31758.02988274707 -Standard Deviation: 386.3622284454486 -Panel Pixel Count: 14925 -Saturated Pixel Count: 0 --
Copyright (c) 2017-2018 MicaSense, Inc. For licensing information see the project git repository
- -