Cloud Removal code

cloud-removal.md

Cloud-removal

A recurring problem with raster processing is the need to remove clouds. With this in mind, code for generating cloud-free satellite imagery from a set of images of a geographical location is included within geotrellis.raster

cloudRemovalMultiband is the important function here. It takes in an array of MultibandTile, which are the GeoTIFF tiles we need to operate on and an optional threshold parameter (which specifies the pixel intensity value below which the resultant cloudless-pixels' intensities would lie). The function returns a processed MultibandTile that can be rendered as a PNG.

Here's an example of its use:

import geotrellis.raster._
import spire.syntax.cfor._

def main(args: Array[String]) : Unit = {
  val dirRed = new File(args(0))
  val dirGreen = new File(args(1))
  val dirBlue = new File(args(2))

  val fileListRed = dirRed.listFiles.filter(_.isFile).toList.toArray
  val fileListGreen = dirGreen.listFiles.filter(_.isFile).toList.toArray
  val fileListBlue = dirBlue.listFiles.filter(_.isFile).toList.toArray

  val numImages = fileListRed.length

  // Should have an equal number of R, G, B tiles
  assert(numImages == fileListBlue.length && numImages == fileListGreen.length)

  val multibands = Array.ofDim[MultibandTile](numImages)

  cfor(0)(_ < numImages, _ + 1) { i =>
    val red = SinglebandGeoTiff(fileListRed(i).toString).tile
    val green = SinglebandGeoTiff(fileListGreen(i).toString).tile
    val blue = SinglebandGeoTiff(fileListBlue(i).toString).tile

    multibands(i) = ArrayMultibandTile(Array(red, green, blue))
  }

  val cloudless = cloudRemovalMultiband(multibands)
  cloudless.renderPng().write("/tmp/cloudless.png")
}

CloudRemoval.scala

package geotrellis.raster.imagery

import geotrellis.raster._
import geotrellis.raster.render._
import geotrellis.raster.io.geotiff.SinglebandGeoTiff

import java.io.File
import spire.syntax.cfor._


/**
  * An object which houses various functions related to cloud removal.
  */
object CloudRemoval {

  /**
    * Attempt to remove clouds by averaging the values that are below
    * a given threshold at a particular point over a series of images.
    */
  def cloudRemovalSingleband(images: Array[Tile], threshold: Int) : Tile = {
    val headImage = images(0)
    val result = ArrayTile.empty(headImage.cellType, headImage.cols, headImage.rows)

    cfor(0)(_ < result.rows, _ + 1) { row =>
      cfor(0)(_ < result.cols, _ + 1) { col =>
        var sum = 0
        var count = 0
        cfor(0)(_ < images.length, _ + 1) { i =>
          val v = images(i).get(col, row)
          if(isData(v) && v < threshold) {
            sum += v
            count += 1
          }
        }
        result.set(col, row, sum / count)
      }
    }

    result
  }

  /**
    * Attempt to remove clouds by averaging the values that are below
    * a given threshold at a particular point over a series of images.
    * This is done on a band-by-band basis.
    */
  def cloudRemovalMultiband(images: Array[MultibandTile], threshold: Int): MultibandTile = {
    val numBands = images(0).bandCount
    val numImages = images.length

    val cloudlessTiles = new Array[Tile](numBands)

    cfor(0)(i => i < numBands, i => i + 1) { i =>
      val singleTiles = new Array[Tile](numImages)

      cfor(0)(j => j < numImages, j => j + 1) { j =>
        singleTiles(j) = images(j).band(i)
      }
      cloudlessTiles(i) = cloudRemovalSingleband(singleTiles, threshold)
    }

    ArrayMultibandTile(cloudlessTiles)
  }

  /**
    * Attempt to remove clouds by averaging the values that are below
    * a given threshold at a particular point over a series of images.
    * This is done on a band-by-band basis.
    */
  def cloudRemovalMultiband(images: Array[MultibandTile]): MultibandTile = {
    cloudRemovalMultiband(images, 10000)
  }

}

CloudRemovalSpec.scala

package geotrellis.raster.imagery

import geotrellis.raster.{ArrayMultibandTile, MultibandTile}
import geotrellis.raster.io.geotiff.{GeoTiffTestUtils, SinglebandGeoTiff}
import geotrellis.raster.testkit.RasterMatchers
import org.scalatest.FunSpec
import spire.syntax.cfor._

class CloudRemovalSpec extends FunSpec
    with RasterMatchers
    with GeoTiffTestUtils {

  describe("Checking cloud removal") {
    it("Pixel value should be less than original cloudy image") {

      val numImages = 3
      val multibands = Array.ofDim[MultibandTile](numImages)

      cfor(0)(_ < numImages, _ + 1) { i =>
        val red = SinglebandGeoTiff(geoTiffPath("cloud_images/red/" + (i+1) + ".TIF")).tile
        val green = SinglebandGeoTiff(geoTiffPath("cloud_images/green/" + (i+1) + ".TIF")).tile
        val blue = SinglebandGeoTiff(geoTiffPath("cloud_images/blue/" + (i+1) + ".TIF")).tile

        multibands(i) = ArrayMultibandTile(Array(red, green, blue))
      }

      // A cloudy pixel
      //print(multibands(1).band(0).get(400, 100), multibands(1).band(1).get(400, 100), multibands(1).band(2).get(400, 100))

      val cloudless = CloudRemoval.cloudRemovalMultiband(multibands)

      // Pixel value after cloud-removal
      assert(cloudless.band(0).get(400, 100) <= multibands(1).band(0).get(400, 100) &&
            cloudless.band(1).get(400, 100) <= multibands(1).band(1).get(400, 100) &&
            cloudless.band(2).get(400, 100) <= multibands(1).band(2).get(400, 100))

    }

    it("Pixel value should be less than threshold") {

      val numImages = 3
      val multibands = Array.ofDim[MultibandTile](numImages)

      cfor(0)(_ < numImages, _ + 1) { i =>
        val red = SinglebandGeoTiff(geoTiffPath("cloud_images/red/" + (i+1) + ".TIF")).tile
        val green = SinglebandGeoTiff(geoTiffPath("cloud_images/green/" + (i+1) + ".TIF")).tile
        val blue = SinglebandGeoTiff(geoTiffPath("cloud_images/blue/" + (i+1) + ".TIF")).tile

        multibands(i) = ArrayMultibandTile(Array(red, green, blue))
      }

      val threshold = 15000
      val cloudless = CloudRemoval.cloudRemovalMultiband(multibands, threshold)

      // Pixel value after cloud-removal
      assert(cloudless.band(0).get(400, 100) <= threshold &&
        cloudless.band(1).get(400, 100) <= threshold &&
        cloudless.band(2).get(400, 100) <= threshold)
    }

    it("Overloaded functions should give the same result for a specific threshold") {

      val numImages = 3
      val multibands = Array.ofDim[MultibandTile](numImages)

      cfor(0)(_ < numImages, _ + 1) { i =>
        val red = SinglebandGeoTiff(geoTiffPath("cloud_images/red/" + (i+1) + ".TIF")).tile
        val green = SinglebandGeoTiff(geoTiffPath("cloud_images/green/" + (i+1) + ".TIF")).tile
        val blue = SinglebandGeoTiff(geoTiffPath("cloud_images/blue/" + (i+1) + ".TIF")).tile

        multibands(i) = ArrayMultibandTile(Array(red, green, blue))
      }

      val threshold = 10000
      val cloudless1 = CloudRemoval.cloudRemovalMultiband(multibands)
      val cloudless2 = CloudRemoval.cloudRemovalMultiband(multibands, threshold)

      // Pixel value after cloud-removal
      assert(cloudless1.band(0).get(400, 100) == cloudless2.band(0).get(400, 100) &&
        cloudless1.band(1).get(400, 100) == cloudless2.band(1).get(400, 100) &&
        cloudless1.band(2).get(400, 100) == cloudless2.band(2).get(400, 100))
    }
  }
}