Geotrellis入门
2018-12-04 本文已影响0人
炼狱的吹笛人
一. 看起来很牛逼的几个链接
二. Geotrellis github
三. 使用
1. 建立工程
按道理来讲，用sbt建立一个scala工程，并在build里引入依赖，就可以自动下载geotrellis的相关包，但我的sbt不知道怎么回事，依赖包下载不来，最后还是通过maven来解决的。
geotrellis的maven列表
先用raster试试水，pom中加入
<!-- https://mvnrepository.com/artifact/org.locationtech.geotrellis/geotrellis-spark -->
<dependency>
<groupId>org.locationtech.geotrellis</groupId>
<artifactId>geotrellis-spark_2.11</artifactId>
<version>2.1.0</version>
</dependency>
<!-- https://mvnrepository.com/artifact/org.locationtech.geotrellis/geotrellis-raster -->
<dependency>
<groupId>org.locationtech.geotrellis</groupId>
<artifactId>geotrellis-raster_2.11</artifactId>
<version>2.1.0</version>
</dependency>
然后便是照抄官方文档里的hello raster示例。
package demo
import geotrellis.raster._
import geotrellis.raster.mapalgebra.focal.Square
import geotrellis.spark._
object Main {
def helloSentence = "Hello GeoTrellis"
def helloRaster(): Unit = {
val nd = NODATA //-2147483648
val input = Array[Int](
nd, 7, 1, 1, 3, 5, 9, 8, 2,
9, 1, 1, 2, 2, 2, 4, 3, 5,
3, 8, 1, 3, 3, 3, 1, 2, 2,
2, 4, 7, 1, nd, 1, 8, 4, 3)
//将数组转化为4*9矩阵
val iat = IntArrayTile(input, 9, 4)
//用一个n*n的窗口对矩阵做卷积，设中心值为平均值
//Square(i) => n = 2 * i + 1
val focalNeighborhood = Square(1)
println(focalNeighborhood)
val meanTile = iat.focalMean(focalNeighborhood)
for (i <- 0 to 3) {
for (j <- 0 to 8) {
print(meanTile.getDouble(j, i) + " ")
}
println()
}
}
def main(args: Array[String]): Unit = {
helloRaster()
}
}
输出结果如下
0 0 0
0 0 0
0 0 0
5.666666666666667 3.8 2.1666666666666665 1.6666666666666667 2.5 4.166666666666667 5.166666666666667 5.166666666666667 4.5
5.6 3.875 2.7777777777777777 1.8888888888888888 2.6666666666666665 3.5555555555555554 4.111111111111111 4.0 3.6666666666666665
4.5 4.0 3.111111111111111 2.5 2.125 3.0 3.111111111111111 3.5555555555555554 3.1666666666666665
4.25 4.166666666666667 4.0 3.0 2.2 3.2 3.1666666666666665 3.3333333333333335 2.75
2.读本地Geotiff文件
一共四行代码
import geotrellis.raster.io.geotiff.reader.GeoTiffReader
import geotrellis.raster.io.geotiff._
val path: String = "filepath/filename.tif"
val geoTiff: SinglebandGeoTiff = GeoTiffReader.readSingleband(path)
这种方法用于读取但波段tif影像，要读取多波段影像，可以用
val geoTiff: MultibandGeoTiff = GeoTiffReader.readMultiband(path)
如果强行用GeoTiffReader.readSingleband(path)方法去读取一个多波段影像，则最后的返回结果是一个单波段影像，且其中的数据为原始影像中的第一个波段。
此外，也可以用影像路径为参数直接构造一个Geotiff对象
import geotrellis.raster.io.geotiff.SinglebandGeoTiff
val geoTiff: SinglebandGeoTiff = SinglebandGeoTiff(path)
3.ETL工具
3.1 调用工具，生成金字塔图层
这个工具貌似是用来进行数据转换的，试着按照各路文档做了一个将landsat8的tif影像打成金字塔的过程。需要在maven中引入
<!-- https://mvnrepository.com/artifact/org.locationtech.geotrellis/geotrellis-spark-etl -->
<dependency>
<groupId>org.locationtech.geotrellis</groupId>
<artifactId>geotrellis-spark-etl_2.11</artifactId>
<version>2.1.0</version>
</dependency>
函数主体很短
def main(args: Array[String]): Unit = {
var args = Array[String](
"--input",
"D:\\Javaworkspace\\geotrellis\\study\\testdata\\input.json",
"--output",
"D:\\Javaworkspace\\geotrellis\\study\\testdata\\output.json",
"--backend-profiles",
"D:\\Javaworkspace\\geotrellis\\study\\testdata\\backend-profiles.json"
);
Logger.getLogger("org").setLevel(Level.ERROR)
implicit val sc = SparkUtils.createSparkContext("ETL", new SparkConf(true).setMaster("local[*]"))
try {
Etl.ingest[ProjectedExtent, SpatialKey, Tile](args)
} finally {
sc.stop
}
}
重点是args里面指定的三个json文件的配置，首先是input.json
[{
"format":"geotiff",
"name":"etlTest1",
"cache":"NONE",
"backend":{
"type":"hadoop",
"path":"D:\\Javaworkspace\\geotrellis\\study\\testdata\\hzpan.TIF"
}
}]
这里面制定了输入数据的格式（geotiff），要输出的图层名（etlTest1），读取的方式（type）和数据路径（path）等参数。对于output.json
{
"backend": {
"type": "file",
"path": "D:\\Javaworkspace\\geotrellis\\study\\testdata\\output1"
},
"reprojectMethod": "buffered",
"pyramid": true,
"tileSize": 256,
"keyIndexMethod": {
"type": "zorder"
},
"resampleMethod": "nearest-neighbor",
"layoutScheme": "zoomed",
"crs": "EPSG:4326"
}
我的windows本地跑spark和hadoop有个问题一直没解决，导致无法用hadoop的输出接口，因此将type属性指定为file。其他属性都可根据实际情况配置。这里官方文档有一个坑，当pyramid设置为true的时候，resampleMethod就不能设置为cubic-spline，否则会报错。最后是backend-profiles.json
{
"backend-profiles":[]
}
貌似是当你使用了accumulo或cassandra来进行输入输出的时候需要配置，一般情况下设个空数组就行。运行成功后在输出目录下新建了两个文件夹
- attributes
- etlTest1
第一个文件夹里存放的似乎是切片后每一个层级的元数据信息，第二个文件夹的命名来自于input.json中的name，其中存放了切片后的图层文件。
然而，上面的代码只能用于处理单波段影响，即便输入的是多波段影像，该工具也只会对第一个波段进行转换。如果想要处理多波段影像，至少应该修改两个地方：
（1）首先是input.json的format：
"format":"multiband-geotiff"
（2）在调用ETL的核心代码中，需要修改切片类型（Tile -> MultibandTile）：
Etl.ingest[ProjectedExtent, SpatialKey, MultibandTile](args)
3.2 渲染切片
上一步中生成的切片文件不能用普通查看图片的方法打开，需要借助Geotrellis提供的类，发布成地图服务或者渲染成图片格式文件（如jpg、png等）。这一步中，尝试将其渲染为图片。
第一步是读取图层。Geotrellis支持多种读取方式，如Accumulo、h3等，这里使用最基本的文件读取。
val zoomId = 9 //要读取的图层层级
implicit val sc = SparkUtils.createSparkContext("ReadLayer", new SparkConf(true).setMaster("local[*]"))
val path = "D:\\Javaworkspace\\geotrellis\\study\\testdata\\output1" //图层文件根目录
val store = FileAttributeStore(path)
val reader = FileLayerReader(path)
val layerId = LayerId("etlTest1", zoomId) //设置图层名称
val layers: TileLayerRDD[SpatialKey] = reader.read[SpatialKey, Tile, TileLayerMetadata[SpatialKey]](layerId)
正确读进来，layers将会是个TileLayerRDD对象。，这里有个坑，在import的时候，要把geotrellis.spark.io下的所有类文件都引进来，不要只引用到的AttributeStore, LayerReader，否则会报找不到隐式参数的错。
//错误示范，会提示 no implicits found for parameter
import geotrellis.spark.io.{AttributeStore, LayerReader}
//正确方法
import geotrellis.spark.io._
获取到RDD后，可以利用spark的foreach算子将其输出到本地的文件目录中
//定义色带，非必须
val colorMap1 = ColorMap(Map(
0 -> RGB(0,0,0),
1 -> RGB(255,255,255)
))
val colorRamp = ColorRamp(RGB(0,0,0), RGB(255,255,255))
.stops(100)
.setAlphaGradient(0xFF, 0xAA)
val outputPath = "D:\\Javaworkspace\\geotrellis\\study\\testdata\\output1\\render\\" + zoomId //图片输出路径
val zoomDir: File = new File(outputPath)
if (!zoomDir.exists()) {
zoomDir.mkdirs()
}
layers.foreach(layer => {
val key = layer._1
val tile = layer._2
val layerPath = outputPath + "\\" + key.row + "_" + key.col + ".jpg"
tile.renderJpg(colorRamp).write(layerPath) //调用渲染方法，colorRamp为非必须参数
})
sc.stop()
程序正确运行后，在相应的输出目录可以看到如下文件
图层渲染结果
将上面的代码稍加修改，就可以将切片导出成tiff影像。首先，需要构造一个GeoTiff对象，它有这样一个构造函数
def apply(tile: Tile, extent: Extent, crs: CRS): SinglebandGeoTiff =
SinglebandGeoTiff(tile, extent, crs)
tile类型的对象可以像上面一样通过RDD的foreach算子取得，CRS参数表示坐标系，可以通过EPSG取得，也可以自己构造，同时在geotrellis里也有预定义几个：
/**
* 4326
*/
object LatLng extends CRS with ObjectNameToString {
lazy val proj4jCrs = factory.createFromName("EPSG:4326")
def epsgCode: Option[Int] = Some(4326)
}
/**
* 3857
*/
object WebMercator extends CRS with ObjectNameToString {
lazy val proj4jCrs = factory.createFromName("EPSG:3857")
def epsgCode: Option[Int] = Some(3857)
}
/**
* Sinusoidal projection, commonly used with MODIS-based data products.
*/
object Sinusoidal extends CRS with ObjectNameToString {
lazy val proj4jCrs: CoordinateReferenceSystem = factory.createFromParameters(null,
"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")
val epsgCode: Option[Int] = None
}
最后的Extent对象，也可以通过TileLayerRDD取得，但首先需要取得该切片对应的图层定义：
val layoutDefinition: LayoutDefinition = layers.metadata.layout
layoutDefinition中有一个mapTransform方法，传入一个SpatialKey，就能得到其对应的格网范围Extend，因此，将字节码图层文件转换为带有空间信息的tiff影像的代码如下：
import geotrellis.raster.io.geotiff.GeoTiff
val layers: TileLayerRDD[SpatialKey] = /** 读取图层代码同上，此处略去 **/
val layoutDefinition: LayoutDefinition = layers.metadata.layout
layers.foreach(layer => {
val key = layer._1
val tile = layer._2
val tiffPath = outputPath + "\\" + key.row + "_" + key.col + ".tiff
GeoTiff(tile, layoutDefinition.mapTransform(key), LatLng).write(tiffPath)
})
采用的是4326即WGS84坐标，输出结果与上面的渲染成jpg类似，但具有空间信息，能被EnVI等遥感影像处理软件正确加载。
4. GeoJson的序列化与反序列化
import geotrellis.vector.{Polygon}
import geotrellis.vector.io._
import geotrellis.vector.io.json._
//serializing
val json = Polygon((10.0,10.0),(10.0,20.0),(30.0,30.0),(10.0,10.0)).toGeoJson
println(json)
//deserializing
val fc: String = """{
| "type": "FeatureCollection",
| "features": [
| {
| "type": "Feature",
| "geometry": { "type": "Point", "coordinates": [1.0, 2.0] },
| "properties": { "someProp": 14 },
| "id": "target_12a53e"
| }, {
| "type": "Feature",
| "geometry": { "type": "Point", "coordinates": [2.0, 7.0] },
| "properties": { "someProp": 5 },
| "id": "target_32a63e"
| }
| ]
|}""".stripMargin
val geos = fc.parseGeoJson[JsonFeatureCollectionMap]
val points1 = geos.getAllPoints()
for (k <- points1.keySet){
println(k + ": " + points1.get(k).get)
}
输出结果如下
{"type":"Polygon","coordinates":[[[10.0,10.0],[10.0,20.0],[30.0,30.0],[10.0,10.0]]]}
target_12a53e: POINT (1 2)
target_32a63e: POINT (2 7)
