Download example landsat data

dl-landsat-example.py

import os
from subprocess import call

images = ["LC80160292015156LGN00",
"LC80160292015124LGN00",
"LC80160292015092LGN00",
"LC80160292015060LGN00",
"LC80160292015028LGN00",
"LC80160292014361LGN00",
"LC80160292014329LGN00",
"LC80160292014297LGN00",
"LC80160292014265LGN00",
"LC80160292014233LGN00",
"LC80160292014201LGN00",
"LC80160292014169LGN00",
"LC80160292014137LGN00",
"LC80160292014089LGN00",
"LC80160292014057LGN00",
"LC80160292014009LGN00",
"LC80160292013342LGN00",
"LC80160292013310LGN00",
"LC80160292013278LGN00",
"LC80160292013246LGN00",
"LC80160292013214LGN00",
"LC80160292013182LGN00",
"LC80160292013150LGN00",
"LC80160292013118LGN02",
"LC80160292015140LGN00",
"LC80160292015108LGN00",
"LC80160292015076LGN00",
"LC80160292015044LGN00",
"LC80160292015012LGN00",
"LC80160292014345LGN00",
"LC80160292014313LGN00",
"LC80160292014281LGN00",
"LC80160292014249LGN00",
"LC80160292014217LGN00",
"LC80160292014185LGN00",
"LC80160292014153LGN00",
"LC80160292014121LGN00",
"LC80160292014073LGN00",
"LC80160292014041LGN00",
"LC80160292013358LGN00",
"LC80160292013326LGN00",
"LC80160292013294LGN00",
"LC80160292013262LGN00",
"LC80160292013230LGN00",
"LC80160292013198LGN00",
"LC80160292013166LGN04",
"LC80160292013134LGN03",
"LC80160292013102LGN01"
]

for image in images:
  cmd = "landsat download " + image
  call(cmd, shell=True)

find-common-extent.scala

import geotrellis.raster._
import geotrellis.raster.io.geotiff._
import geotrellis.raster.io.geotiff.reader._
import geotrellis.raster.io.json._
import geotrellis.vector._
import spray.json._
import spray.json.DefaultJsonProtocol._

def fname(i: Int) = {
  val f = f"/home/rob/landsat/downloads/image$i%02d/B2.TIF"
  if(new java.io.File(f).exists) Some(f) else None
}

def combine(re1: RasterExtent, re2: RasterExtent): RasterExtent = {
  val e1 = re1.extent
  val e2 = re2.extent
  e1 & e2 match {
    case Some(e) =>
      val x = re1.createAligned(e)
      if(!e2.contains(x.extent)) re2.createAligned(e)
      else x
    case None =>
      sys.error("Some things done match.")
  }
}

val RasterExtent(Extent(xmin, ymin, xmax, ymax), _, _, cols, rows) = 
  (0 until 4)
    .map(fname)
    .flatten
    .map { path =>
      println(s"Reading $path")
      val tags = TiffTagsReader.read(path)
      RasterExtent(tags.extent, tags.cols, tags.rows)
     }
    .reduce(combine)

println(f"gdalwarp -te $xmin%.3f $ymin%.3f $xmax%.3f $ymax%.3f -ts $cols $rows")

unzip.py

import os, shutil, rasterio
from subprocess import call

l = enumerate(os.listdir('.'))
#l = [(0, 'LC80160292013294LGN00.tar.bz')]
#l = [(0, 'LC80160292013310LGN00.tar.bz')]

# Bands:
#  blue - 2
#  green - 3
#  red - 4
# near ir - 5
# cloud - 9
def keepers(n):
  bs = ['%s_B2.TIF', '%s_B3.TIF', '%s_B4.TIF', '%s_B5.TIF', '%s_B9.TIF']
  return map(lambda x: x % n, bs)

def process_unzipped(n):
  unzipped = os.listdir('.')
  kprs = keepers(n)
  for a in unzipped:
    if not a in kprs:
      os.remove(a)
    else:
      cmd = "gdal_translate -co compress=deflate -co predictor=2 %s %s" % (a, a[-6:])
      call(cmd, shell=True)
  call("rm LC*", shell=True)
    


def process_zip(n, i):
  p = "image%02d" % i
  if not os.path.exists(p):
    os.mkdir(p)
  new_path = os.path.join(p, x)
  shutil.move(x, new_path)
  os.chdir(p)
  cmd = "tar -xvjf %s" % n
  call(cmd, shell=True)
  shutil.move(x, '..')
  process_unzipped(n.replace(".tar.bz",""))

r = os.path.abspath('.')
for i, x in l:
    process_zip(x, i)
    os.chdir(r)

warp.py

import os
from subprocess import call

warpcmd = "gdalwarp -te 263700.000 4826370.000 492630.000 5053800.000 -ts 7631 7581 -co compress=deflate -co predictor=2 -r bilinear %s %s"

i = 0
ptoi = { }
for root, folders, files in os.walk('.'):
    for f in filter(lambda x: x.startswith('B') and x.endswith('.TIF'), files):
        p = os.path.abspath(os.path.join(root, f))
        b = os.path.basename(os.path.dirname(p))
        if not b in ptoi:
            ptoi[b] = "image%02d" % i
            i += 1
        d = "/home/rob/landsat/downloads/processed/%s" % ptoi[b]
        if not os.path.exists(d):
            os.mkdir(d)
        p2 = os.path.join(d, f)
        cmd = warpcmd % (p, p2)
        print "Working on %d %s/%s" % (i, b, f)
        call(cmd, shell=True)
        
        #print cmd

# 7891, 8021
# 7831, 8021