Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Cloud Removal code

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")
}
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)
}
}
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))
}
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.