Mask a GeoTiff

In [1]:

    

import sys.process._
import geotrellis.proj4.CRS
import geotrellis.raster.io.geotiff.writer.GeoTiffWriter
import geotrellis.raster.io.geotiff.{SinglebandGeoTiff, _}
import geotrellis.raster.{CellType, DoubleArrayTile}
import geotrellis.spark.io.hadoop._
import geotrellis.vector.{Extent, ProjectedExtent}
import org.apache.spark.mllib.linalg.Vector
import org.apache.spark.rdd.RDD
import org.apache.spark.{SparkConf, SparkContext}

//Spire is a numeric library for Scala which is intended to be generic, fast, and precise.
import spire.syntax.cfor._


    

In [2]:

    

var geo_projected_extent = new ProjectedExtent(new Extent(0,0,0,0), CRS.fromName("EPSG:3857"))
var geo_num_cols_rows :(Int, Int) = (0, 0)
val geo_path = "hdfs:///user/hadoop/modis/Onset_Greenness_Increase/A2001001__Onset_Greenness_Increase.tif"

val geo_tiles_RDD = sc.hadoopGeoTiffRDD(geo_path).values

val geo_extents_withIndex = sc.hadoopMultibandGeoTiffRDD(geo_path).keys.zipWithIndex().map{case (e,v) => (v,e)}
geo_projected_extent = (geo_extents_withIndex.filter(m => m._1 == 0).values.collect())(0)

val geo_tiles_withIndex = geo_tiles_RDD.zipWithIndex().map{case (e,v) => (v,e)}
val geo_tile0 = (geo_tiles_withIndex.filter(m => m._1==0).values.collect())(0)
geo_num_cols_rows = (geo_tile0.cols, geo_tile0.rows)
val geo_cellT = geo_tile0.cellType


    

    






Waiting for a Spark session to start...






    






geo_projected_extent = ProjectedExtent(Extent(-171.1296209227216, 19.996701428244357, -42.56469205974807, 49.99999999547987),geotrellis.proj4.CRS$$anon$3@41d0d1b7)
geo_num_cols_rows = (17800,4154)
geo_path = hdfs:///user/hadoop/modis/Onset_Greenness_Increase/A2001001__Onset_Greenness_Increase.tif
geo_tiles_RDD = MapPartitionsRDD[2] at values at <console>:50
geo_extents_withIndex = MapPartitionsRDD[7] at map at <console>:52
geo_projected_extent = ProjectedExtent(Extent(-171.1296209227216, 19.996701428244357, -42.56469205974807, 49.99999999547987),geotrellis.proj4.CRS$$anon$3@41d0d1b...







    
Out[2]:





ProjectedExtent(Extent(-171.1296209227216, 19.996701428244357, -42.56469205974807, 49.99999999547987),geotrellis.proj4.CRS$$anon$3@41d0d1b7)

In [3]:

    

val classes = Array(32767)

val classesBro = sc.broadcast(classes)
val res_tile = geo_tiles_RDD.map(m => m.toArrayDouble()).map(m => m.map(v => if (classesBro.value.contains(v)) 1.0 else 0.0))


    

    






classes = Array(32767)
classesBro = Broadcast(4)
res_tile = MapPartitionsRDD[15] at map at <console>:52







    
Out[3]:





MapPartitionsRDD[15] at map at <console>:52

In [4]:

    

val clone_tile = DoubleArrayTile(res_tile.collect()(0), geo_num_cols_rows._1, geo_num_cols_rows._2)

val cloned = geotrellis.raster.DoubleArrayTile.empty(geo_num_cols_rows._1, geo_num_cols_rows._2)
cfor(0)(_ < geo_num_cols_rows._1, _ + 1) { col =>
    cfor(0)(_ < geo_num_cols_rows._2, _ + 1) { row =>
        val v = clone_tile.getDouble(col, row)
        cloned.setDouble(col, row, v)
    }
}

val geoTif = new SinglebandGeoTiff(cloned, geo_projected_extent.extent, geo_projected_extent.crs, Tags.empty, GeoTiffOptions(compression.DeflateCompression))

//Save GeoTiff to /tmp
val output = "hdfs:///user/pheno/modis/Onset_Greenness_Increase/A2001001__Onset_Greenness_Increase.tif"
val tmp_output = "/tmp/A2001001__Onset_Greenness_Increase.tif"

GeoTiffWriter.write(geoTif, tmp_output)

//Upload to HDFS
var cmd = "hadoop dfs -copyFromLocal -f " + tmp_output + " " + output
Process(cmd)!

cmd = "rm -fr " + tmp_output
Process(cmd)!


    

    




DEPRECATED: Use of this script to execute hdfs command is deprecated.
Instead use the hdfs command for it.







    






clone_tile = DoubleConstantNoDataArrayTile([D@2ea6dcc8,17800,4154)
cloned = DoubleConstantNoDataArrayTile([D@108690ef,17800,4154)
geoTif = SinglebandGeoTiff(DoubleConstantNoDataArrayTile([D@108690ef,17800,4154),Extent(-171.1296209227216, 19.996701428244357, -42.56469205974807, 49.99999999547987),geotrellis.proj4.CRS$$anon$3@41d0d1b7,Tags(Map(),List()),GeoTiffOptions(geotrellis.raster.io.geotiff.Striped@7a8ae72b,geotrellis.raster.io.geotiff.compression.DeflateCompression$@5230299e,1,None))
output = hdfs:///user/pheno/modis/Onset_Greenness_Increase/A2001001__Onset_Greenness_Increase.tif
tmp_output = /tmp/A2001001__Onset_Greenness_Increase.tif
cmd = r...







    






warning: there were two feature warnings; re-run with -feature for details







    
Out[4]:





rm -fr /tmp/A2001001__Onset_Greenness_Increase.tif

In [ ]: