Skip to content

Instantly share code, notes, and snippets.

@brunifrancesco
Created April 16, 2019 11:02
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save brunifrancesco/bd5fe37238357e39d4d5bfa305c4c408 to your computer and use it in GitHub Desktop.
Save brunifrancesco/bd5fe37238357e39d4d5bfa305c4c408 to your computer and use it in GitHub Desktop.
Ingest and visualize landsat data
import rasterio
import geopyspark as gps
import numpy as np
import matplotlib
import matplotlib.cm as cm
from pyspark import SparkContext
conf = gps.geopyspark_conf(master="local[*]", appName="ingest-example", )
conf.set(key='spark.kryoserializer.buffer.max', value='256m')
conf.set(key='spark.driver.maxResultSize', value="8g")
conf.set(key='spark.driver.memory', value="8g")
pysc = SparkContext(conf=conf)
raster_layer_b3 = gps.geotiff.get(layer_type=gps.LayerType.SPATIAL, uri="file:///landsat_data/LC81070352015218LGN00_B3.TIF")
raster_layer_b5 = gps.geotiff.get(layer_type=gps.LayerType.SPATIAL, uri="file:///landsat_data/LC81070352015218LGN00_B5.TIF")
raster_layer_bqa = gps.geotiff.get(layer_type=gps.LayerType.SPATIAL, uri="file://landsat_data/LC81070352015218LGN00_BQA.TIF")
combined_layer = gps.combine_bands(
[raster_layer_b3, raster_layer_b5, raster_layer_bqa]
).convert_data_type(new_type=gps.CellType.INT32)
combined_layer.collect_metadata()
tiled_raster_layer = combined_layer.tile_to_layout(gps.GlobalLayout(), target_crs=3857)
for layer in tiled_raster_layer.pyramid().levels.values():
gps.write(uri="file:///tmp/bbb", layer_name="japan", tiled_raster_layer=layer)
cm = gps.ColorMap.nlcd_colormap()
layers = []
for zoom in range(10, 13):
layers.append(gps.query(uri="file:///tmp/bbb",
layer_name="japan",
layer_zoom=zoom))
def render_image(tiles):
print('here')
arr = tiles[0].cells[0]
arr1 = tiles[1].cells[0]
res = arr/arr1
norm = matplotlib.colors.Normalize(vmin=min(res), vmax=max(res), clip=True)
mapper = cm.ScalarMappable(norm=norm, cmap=cm.Greys_r)
return Image.fromarray(mapper.to_rgba(res), mode='RGBA')
landsat_pyramid = gps.Pyramid(layers)
tms = gps.TMS.build(source=landsat_pyramid, display=render_image)
tms.bind(host='0.0.0.0', requested_port=8889)
tms.url_pattern
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment