This is an automated email from the ASF dual-hosted git repository.
fgreg pushed a commit to branch SDAP-149
in repository https://gitbox.apache.org/repos/asf/incubator-sdap-nexus.git
The following commit(s) were added to refs/heads/SDAP-149 by this push:
new f498e9c SDAP-149 Generate images directly from data
f498e9c is described below
commit f498e9cbd9eef702af138235656af6cb416b2420
Author: Kevin M. Gill <[email protected]>
AuthorDate: Mon Nov 5 09:55:14 2018 -0800
SDAP-149 Generate images directly from data
Initial work.
---
analysis/setup.py | 2 +-
.../algorithms/imaging/ColorTableHandler.py | 26 ++-
.../algorithms/imaging/LayerListHandler.py | 143 +++++++++++++-
.../webservice/algorithms/imaging/colortable.py | 7 +-
.../webservice/algorithms/imaging/colortables.py | 73 ++++++-
analysis/webservice/algorithms/imaging/layer.py | 135 +++++++++++++
.../webservice/algorithms/imaging/mapprocessing.py | 4 +
.../webservice/algorithms/imaging/mrf_reader.py | 85 ++++----
.../webservice/algorithms/imaging/reprojection.py | 220 +++++++++++++++------
data-access/requirements.txt | 2 +-
data-access/setup.py | 4 +-
docker/nexus-webapp/Dockerfile | 1 +
docker/spark-mesos-base/Dockerfile | 1 +
13 files changed, 582 insertions(+), 121 deletions(-)
diff --git a/analysis/setup.py b/analysis/setup.py
index 42903cb..1603073 100644
--- a/analysis/setup.py
+++ b/analysis/setup.py
@@ -16,7 +16,7 @@
import setuptools
-__version__ = '1.5'
+__version__ = '1.6'
setuptools.setup(
name="nexusanalysis",
diff --git a/analysis/webservice/algorithms/imaging/ColorTableHandler.py
b/analysis/webservice/algorithms/imaging/ColorTableHandler.py
index 9c798ac..0e06c0a 100644
--- a/analysis/webservice/algorithms/imaging/ColorTableHandler.py
+++ b/analysis/webservice/algorithms/imaging/ColorTableHandler.py
@@ -75,9 +75,23 @@ class ColorTableHandler(BaseHandler):
values.append((value, value_high))
return values
+ def _create_colortable(self, color_table, min_value, max_value, num_items,
units):
+ colors = ColorTableHandler._create_hex_color_list(color_table,
num_items)
+ labels = ColorTableHandler._create_labels(min_value, max_value, units,
num_items)
+ values = ColorTableHandler._create_values(min_value, max_value,
num_items)
+ return colors, labels, values
+
+ def _create_from_sdap(self, ds, color_table, force_min, force_max,
num_items, units):
+ stats = self._tile_service.get_dataset_overall_stats(ds)
+ min_value = force_min if force_min is not None else stats["minValue"]
+ max_value = force_max if force_max is not None else stats["maxValue"]
+ return self._create_colortable(color_table, min_value, max_value,
num_items, units)
+
def calc(self, computeOptions, **args):
ds = computeOptions.get_argument("ds")
+ generic = computeOptions.get_boolean_arg("generic", False)
+
color_table_identifier = computeOptions.get_argument("ct", "rainbow")
color_table = colortables.get_color_table(color_table_identifier)
@@ -88,16 +102,12 @@ class ColorTableHandler(BaseHandler):
num_items = np.array((2048, num_items)).min()
num_items = np.array((255, num_items)).max()
- stats = self._tile_service.get_dataset_overall_stats(ds)
-
units = computeOptions.get_argument("units", "°c")
- min_value = force_min if force_min is not None else stats["minValue"]
- max_value = force_max if force_max is not None else stats["maxValue"]
-
- colors = ColorTableHandler._create_hex_color_list(color_table,
num_items)
- labels = ColorTableHandler._create_labels(min_value, max_value, units,
num_items)
- values = ColorTableHandler._create_values(min_value, max_value,
num_items)
+ if generic is False:
+ colors, labels, values = self._create_from_sdap(ds, color_table,
force_min, force_max, num_items, units)
+ else:
+ colors, labels, values = self._create_colortable(color_table,
force_min, force_max, num_items, units)
return ColorTableResponse({
"units": units,
diff --git a/analysis/webservice/algorithms/imaging/LayerListHandler.py
b/analysis/webservice/algorithms/imaging/LayerListHandler.py
index f5dd078..fcb2f0d 100644
--- a/analysis/webservice/algorithms/imaging/LayerListHandler.py
+++ b/analysis/webservice/algorithms/imaging/LayerListHandler.py
@@ -17,6 +17,11 @@ import json
import ConfigParser
import pkg_resources
import uuid
+import sys
+import uuid
+import json
+import hashlib
+import layer
from webservice.NexusHandler import NexusHandler as BaseHandler
from webservice.NexusHandler import nexus_handler
@@ -47,6 +52,109 @@ class LayerListHandler(BaseHandler):
BaseHandler.__init__(self)
self.imagery_config = ConfigParser.RawConfigParser()
self.imagery_config.readfp(pkg_resources.resource_stream(__name__,
"config.ini"), filename='config.ini')
+ self.static_layers_base = self.imagery_config.get("tileserver",
"mrf.base")
+ self.parse_static_layers()
+
+ def build_colorbar_url(self, colorbarSpec):
+ if colorbarSpec is None:
+ return None
+ url =
"{sdap_base}/imaging/colortable?ds={ds}&ct={ct}&min={min}&max={max}&units={units}&generic=true".format(
+ sdap_base=self.imagery_config.get("imaging", "imaging.endpoint"),
+ ds=colorbarSpec["productLabel"],
+ ct=colorbarSpec["colorbarName"],
+ min=colorbarSpec["minValue"],
+ max=colorbarSpec["maxValue"],
+ units=colorbarSpec["units"]
+ )
+ return url
+
+ def build_layer_for_projection(self, layer_def, product_suffix,
endpoint_proj, layerProjection, tileMatrixSet):
+ layer = {}
+ layer["UUID"] = str(uuid.uuid4())
+ layer["DataID"] = layer_def["productLabel"]
+ layer["ProductLabel"] = "%s%s" % (layer_def["productLabel"],
product_suffix)
+ layer["CrmShortName"] = layer_def["crmShortName"] if
layer_def["crmShortName"] is not None else layer_def[
+ "productLabel"]
+ layer["enabled"] = layer_def["enabled"]
+ layer["mission"] = layer_def["mission"]
+ layer["instrument"] = layer_def["instrument"]
+ layer["parameter"] = layer_def["parameter"]
+ layer["ProductType"] = layer_def["productType"]
+ layer["ServiceProtocol"] = layer_def["serviceProtocol"]
+ layer["EndPoint"] = layer_def["endpoint"].format(proj=endpoint_proj)
+ layer["WMSEndPoint"] = ""
+ layer["wmsLayer"] = "0"
+ layer["LayerTitle"] = layer_def["layerTitle"]
+ layer["LayerService"] = layer_def["layerService"]
+ layer["LayerProjection"] = layerProjection
+ layer["availableForAnalysis"] = layer_def["availableForAnalysis"]
+ layer["analysisDatasetName"] = layer_def["analysisDatasetName"]
+ layer["ThumbnailImage"] = layer_def["thumbnailImage"]
+ layer["bounding"] = {
+ "westbc": layer_def["bounding"]["west"],
+ "eastbc": layer_def["bounding"]["east"],
+ "northbc": layer_def["bounding"]["north"],
+ "southbc": layer_def["bounding"]["south"]
+ }
+ layer["LayerSubtitle"] = layer_def["layerSubtitle"]
+ layer["legend"] = ""
+ layer["wmts"] = {
+ "tileMatrixSet": tileMatrixSet,
+ "tileLayerName": layer_def["wmtsLayerName"],
+ "tileFormat": layer_def["wmtsFormat"]
+ }
+
+ layer["nativeResolution"] = "9km"
+ if layer_def["availability"] is not None:
+ layer["availability"] = {
+ "start": layer_def["availability"]["start"],
+ "end": layer_def["availability"]["end"]
+ }
+ else:
+ layer["availability"] = None
+ layer["keywords"] = layer_def["keywords"]
+ layer["utilityLayer"] = layer_def["utilityLayer"]
+ layer["depths"] = layer_def["depthSpec"]
+ layer["type"] = layer_def["type"]
+ layer["colorbar"] = self.build_colorbar_url(layer_def["colorbar"])
+ return layer
+
+ def build_layer(self, layer_def, basemap):
+ layers = []
+ if layer_def["hasGlobal"] is True and basemap ==
layer_def["baseLayer"] and layer_def["overridesSdap"] is False:
+ layers.append(self.build_layer_for_projection(layer_def,
+ "",
+ "geo",
+
layer_def["layerProjectionGlobal"],
+
layer_def["wmtsGlobalMatrixSet"]))
+ if layer_def["hasNorth"] is True and basemap == layer_def["baseLayer"]
and layer_def["overridesSdap"] is False:
+ layers.append(self.build_layer_for_projection(layer_def,
+ "-arctic",
+ "arctic",
+
layer_def["layerProjectionNorth"],
+
layer_def["wmtsNorthMatrixSet"]))
+ if layer_def["hasSouth"] is True and basemap == layer_def["baseLayer"]
and layer_def["overridesSdap"] is False:
+ layers.append(self.build_layer_for_projection(layer_def,
+ "-antarctic",
+ "antarctic",
+
layer_def["layerProjectionSouth"],
+
layer_def["wmtsSouthMatrixSet"]))
+ return layers
+
+ def build_layer_config(self, layer_list, basemap):
+ compiled_layers = []
+
+ for entry in layer_list:
+ compiled_layers += self.build_layer(entry, basemap)
+
+ return compiled_layers
+
+ def parse_static_layers(self):
+ sys.path.insert(0, self.static_layers_base)
+ import layers
+
+
+
def _build_layer_spec(self, ds):
layer = {
@@ -93,13 +201,40 @@ class LayerListHandler(BaseHandler):
}
return layer
+
+ def apply_layer_overrides(self, sdap_layer):
+ layer_overrides = layer.getLayerByLabel(sdap_layer["ProductLabel"])
+ if layer_overrides["crmShortName"] is not None:
+ sdap_layer["CrmShortName"] = layer_overrides["crmShortName"]
+ if layer_overrides["layerTitle"] is not None:
+ sdap_layer["LayerTitle"] = layer_overrides["layerTitle"]
+ if layer_overrides["layerSubtitle"] is not None:
+ sdap_layer["LayerSubtitle"] = layer_overrides["layerSubtitle"]
+ if layer_overrides["keywords"] is not None:
+ sdap_layer["keywords"] = layer_overrides["keywords"]
+ if layer_overrides["instrument"] is not None:
+ sdap_layer["instrument"] = layer_overrides["instrument"]
+ if layer_overrides["mission"] is not None:
+ sdap_layer["mission"] = layer_overrides["mission"]
+ if layer_overrides["parameter"] is not None:
+ sdap_layer["parameter"] = layer_overrides["parameter"]
+
+
+
def calc(self, computeOptions, **args):
- ds_list = self._tile_service.get_dataseries_list()
+ basemaps = computeOptions.get_boolean_arg("basemaps", False)
layers = []
- for ds in ds_list:
- layer = self._build_layer_spec(ds)
- layers.append(layer)
+
+ if basemaps is False:
+ ds_list = self._tile_service.get_dataseries_list()
+ for ds in ds_list:
+ sdap_layer = self._build_layer_spec(ds)
+ self.apply_layer_overrides(sdap_layer)
+ layers.append(sdap_layer)
+
+ compiled_layers = self.build_layer_config(layer.getLayers(), basemaps)
+ layers += compiled_layers
return LayerListResponse({
"Layers": {
diff --git a/analysis/webservice/algorithms/imaging/colortable.py
b/analysis/webservice/algorithms/imaging/colortable.py
index 7592b4a..3356477 100644
--- a/analysis/webservice/algorithms/imaging/colortable.py
+++ b/analysis/webservice/algorithms/imaging/colortable.py
@@ -15,13 +15,18 @@
import numpy as np
import math
-
+import types
class ColorTable:
def __init__(self, identifier, name, spec):
self.identifier = identifier
self.name = name
+
+ # If the list of colors was passed in as a list of hex strings
("ABCDEFFF"), convert them to RGB values
+ if isinstance(spec[0], types.StringType):
+ spec = [[int(p[i:i+2], 16) for i in range(0, 6, 2)] for p in spec]
+
self.spec = np.array(spec)
def get_color(self, fraction):
diff --git a/analysis/webservice/algorithms/imaging/colortables.py
b/analysis/webservice/algorithms/imaging/colortables.py
index b1a6392..38afe13 100644
--- a/analysis/webservice/algorithms/imaging/colortables.py
+++ b/analysis/webservice/algorithms/imaging/colortables.py
@@ -1282,7 +1282,78 @@
COLORTABLES.append(colortable.ColorTable("NEO_modis_aer_od",
]
))
-
+COLORTABLES.append(colortable.ColorTable("IMERG_Rain_Rate",
+ "IMERG_Rain_Rate",
+ [
+ "bfbfbfff",
+ "007c53ff",
+ "00814dff",
+ "008548ff",
+ "008a41ff",
+ "008e3aff",
+ "009333ff",
+ "00972bff",
+ "009c22ff",
+ "00a019ff",
+ "00a50fff",
+ "00aa05ff",
+ "05ae00ff",
+ "11b300ff",
+ "1db700ff",
+ "29bc00ff",
+ "37c000ff",
+ "44c500ff",
+ "53c900ff",
+ "62ce00ff",
+ "71d200ff",
+ "81d700ff",
+ "92db00ff",
+ "a3e000ff",
+ "b5e400ff",
+ "c8e900ff",
+ "dbee00ff",
+ "eef200ff",
+ "f7eb00ff",
+ "fbdf00ff",
+ "ffd201ff",
+ "ffc305ff",
+ "ffb50aff",
+ "ffa70eff",
+ "ff9913ff",
+ "ff8d17ff",
+ "ff811cff",
+ "ff7520ff",
+ "ff6a25ff",
+ "ff6029ff",
+ "ff562eff",
+ "ff4c33ff",
+ "ff3628ff",
+ "ff1e1eff",
+ "ff1414ff",
+ "ff0a0aff",
+ "ff0000ff",
+ "f40000ff",
+ "ea0000ff",
+ "e00000ff",
+ "d60000ff",
+ "cc0000ff",
+ "c10000ff",
+ "b70000ff",
+ "ad0000ff",
+ "a30000ff",
+ "990000ff",
+ "8e0000ff",
+ "840000ff",
+ "7a0000ff",
+ "700000ff",
+ "660000ff",
+ "5b0000ff",
+ "510000ff",
+ "470000ff",
+ "3d0000ff",
+ "330000ff"
+ ]
+ ))
def get_color_table(identifier):
assert identifier is not None and (isinstance(identifier,
types.StringType) or isinstance(identifier, types.UnicodeType))
diff --git a/analysis/webservice/algorithms/imaging/layer.py
b/analysis/webservice/algorithms/imaging/layer.py
new file mode 100644
index 0000000..5fcb5be
--- /dev/null
+++ b/analysis/webservice/algorithms/imaging/layer.py
@@ -0,0 +1,135 @@
+import time
+from datetime import datetime
+
+__LAYERS__ = []
+
+
+def colorbar(productLabel,
+ minValue,
+ maxValue,
+ colorbarName,
+ units
+ ):
+ return {
+ "productLabel": productLabel,
+ "minValue": minValue,
+ "maxValue": maxValue,
+ "colorbarName": colorbarName,
+ "units": units
+ }
+
+def bounding(north=90.0,
+ south=-90.0,
+ east=180.0,
+ west=-180.0):
+ return {
+ "north": north,
+ "south": south,
+ "east": east,
+ "west": west
+ }
+
+def depthspec(title="Depth", units="m", labels=[], values=[]):
+
+ return {
+ "title": title,
+ "labels": labels,
+ "values": values,
+ "units": units,
+ }
+
+def addlayer(
+ productLabel,
+ crmShortName=None,
+ enabled=True,
+ mission="",
+ instrument="Unspecified",
+ parameter="Unspecified",
+ productType="Mosaic",
+ serviceProtocol="GIBS",
+ endpoint="",
+ layerTitle="",
+ layerService="wmtstiled",
+ hasGlobal=True,
+ hasNorth=True,
+ hasSouth=True,
+ thumbnailImage=None,
+ layerSubtitle="",
+ availableForAnalysis=False,
+ analysisDatasetName="",
+ layerProjectionGlobal="EPSG:4326",
+ layerProjectionNorth="EPSG:3413",
+ layerProjectionSouth="EPSG:3031",
+ wmtsGlobalMatrixSet="EPSG4326",
+ wmtsNorthMatrixSet="EPSG3413",
+ wmtsSouthMatrixSet="EPSG3031",
+ wmtsLayerName="",
+ wmtsFormat="image/png",
+ colorbar=None,
+ availability=None,
+ keywords=[],
+ utilityLayer=False,
+ baseLayer=False,
+ depthSpec=None,
+ bounding=bounding(),
+ type="observation",
+ overridesSdap=False):
+
+ global __LAYERS__
+
+ layer = {
+ "productLabel": productLabel,
+ "crmShortName": crmShortName,
+ "enabled": enabled,
+ "mission": mission,
+ "instrument": instrument,
+ "parameter": parameter,
+ "productType": productType,
+ "serviceProtocol": serviceProtocol,
+ "endpoint": endpoint,
+ "layerTitle": layerTitle,
+ "layerService": layerService,
+ "hasGlobal": hasGlobal,
+ "hasNorth": hasNorth,
+ "hasSouth": hasSouth,
+ "thumbnailImage": thumbnailImage,
+ "bounding": bounding,
+ "layerSubtitle": layerSubtitle,
+ "availableForAnalysis": availableForAnalysis,
+ "analysisDatasetName": analysisDatasetName,
+ "layerProjectionGlobal": layerProjectionGlobal,
+ "layerProjectionNorth": layerProjectionNorth,
+ "layerProjectionSouth": layerProjectionSouth,
+ "wmtsGlobalMatrixSet": wmtsGlobalMatrixSet,
+ "wmtsNorthMatrixSet": wmtsNorthMatrixSet,
+ "wmtsSouthMatrixSet": wmtsSouthMatrixSet,
+ "wmtsLayerName": wmtsLayerName,
+ "wmtsFormat": wmtsFormat,
+ "availability": availability,
+ "keywords": keywords,
+ "utilityLayer": utilityLayer,
+ "baseLayer": baseLayer,
+ "colorbar": colorbar,
+ "depthSpec": depthSpec,
+ "type":type,
+ "overridesSdap": overridesSdap
+ }
+
+ __LAYERS__.append(layer)
+
+
+
+def timestamp(year, month, day):
+ dt = datetime(year, month, day)
+ ts = time.mktime(dt.timetuple())
+ return int(ts * 1000)
+
+
+def getLayers():
+ return __LAYERS__
+
+def getLayerByLabel(label):
+ for l in __LAYERS__:
+ if l["productLabel"] == label:
+ return l
+ return None
\ No newline at end of file
diff --git a/analysis/webservice/algorithms/imaging/mapprocessing.py
b/analysis/webservice/algorithms/imaging/mapprocessing.py
index 3b97dc9..9b1721d 100644
--- a/analysis/webservice/algorithms/imaging/mapprocessing.py
+++ b/analysis/webservice/algorithms/imaging/mapprocessing.py
@@ -29,6 +29,7 @@ import multiprocessing
import colortables
import colorization
+import reprojection
NO_DATA_IMAGE = None
@@ -308,6 +309,9 @@ def process_tiles_to_map(nexus_tiles, stats, reqd_tllr,
width=None, height=None,
data = trim_map_to_requested_tllr(data, reqd_tllr, tiles_tllr)
data = expand_map_to_requested_tllr(data, reqd_tllr, tiles_tllr, x_res,
y_res)
+
+ #data = reprojection.reproject_rgba_map(data,
reprojection.TargetProjection.NORTH)
+
if width is not None and height is not None:
data = imresize(data, (height, width), interp=interpolation)
diff --git a/analysis/webservice/algorithms/imaging/mrf_reader.py
b/analysis/webservice/algorithms/imaging/mrf_reader.py
index cf0267c..addb5bd 100644
--- a/analysis/webservice/algorithms/imaging/mrf_reader.py
+++ b/analysis/webservice/algorithms/imaging/mrf_reader.py
@@ -17,8 +17,13 @@ import os
import math
import xml.etree.ElementTree as ET
import struct
+import numpy as np
+
class MrfMeta:
+ """
+ Implements a basic MRF metadata reader
+ """
def __init__(self, mrf_meta_path):
self.__mrf_meta_path = mrf_meta_path
@@ -62,13 +67,17 @@ class MrfMeta:
self.projection = proj_el.text
rsets = root.find("./Rsets")
- self.rsets_model = rsets.attrib["model"]
- self.rsets_scale = rsets.attrib["scale"]
-
-
+ if rsets is not None:
+ if "model" in rsets.attrib:
+ self.rsets_model = rsets.attrib["model"]
+ if "scale" in rsets.attrib:
+ self.rsets_scale = rsets.attrib["scale"]
class MrfReader:
+ """
+ Implements a basic MRF tile data reader
+ """
def __init__(self, mrf_path):
self.__mrf_path = mrf_path
@@ -87,68 +96,48 @@ class MrfReader:
len_base = w * h
low = 1
- idx_size = int(os.path.getsize(self.__idx_path))
- len_tiles = idx_size / 16
- print "Number of tiles:", len_tiles
self.__levels = []
self.__rows = []
self.__cols = []
- self.__levels.append(0)
- self.__rows.append(h)
- self.__cols.append(w)
-
while len_base > low:
self.__levels.append(len_base)
- w = int(math.ceil(float(w) / 2.0))
- h = int(math.ceil(float(h) / 2.0))
self.__rows.append(h)
self.__cols.append(w)
+
+ w = int(math.ceil(float(w) / 2.0))
+ h = int(math.ceil(float(h) / 2.0))
len_base = w * h
self.__levels.append(1)
+ self.__rows.append(h)
+ self.__cols.append(w)
+ def __read_offset_and_size(self, tile):
+ with open(self.__idx_path, "rb") as idx_file:
+ idx_file.seek(16 * tile, os.SEEK_SET)
+ offset, size = struct.unpack(">2Q", idx_file.read(16))
+ return offset, size
+ def __read_tile_data(self, offset, size):
+ with open(self.__data_path) as mrf_file:
+ mrf_file.seek(offset, os.SEEK_SET)
+ return mrf_file.read(size)
def get_tile_bytes(self, tilematrix, tilerow, tilecol):
- tilematrix = tilematrix + 2
- level = self.__levels[len(self.__levels) - tilematrix]
- row = self.__rows[len(self.__levels) - tilematrix]
- col = self.__cols[len(self.__levels) - tilematrix]
- #print tilerow, tilecol, level, row, col
- #assert tilerow > row - 1, "Tile row exceeds the maximum"
- #assert tilecol > col - 1, "Tole col exceeds the maximum"
-
- level_start = 0
-
- for idx in range(0, len(self.__levels)):
- val = self.__levels[idx]
- if idx > 0 and self.__levels[idx - 1] == level:
- break
- level_start += val
-
- tile = (tilerow * col) + tilecol + level_start
-
- offset = 0
- size = 0
- buffer = None
-
- with open(self.__idx_path) as idx_file:
- idx_file.seek(16 * tile)
- offset = struct.unpack('L', idx_file.read(8))[0]
- size = struct.unpack('L', idx_file.read(8))[0]
-
- with open(self.__mrf_path) as mrf_file:
- print offset, size
- mrf_file.seek(offset)
- buffer = mrf_file.read(size)
+ tilematrix += 1
- return buffer
+ cols = self.__cols[-(tilematrix + 1)]
+ rows = self.__rows[-(tilematrix + 1)]
+ assert tilecol < cols, "Tile column exceeds maximum columns for this
layer/tile matrix"
+ assert tilerow < rows, "Tile row exceeds maximum rows for this
layer/tile matrix"
-if __name__ == "__main__":
+ level_start = np.array(self.__levels[: - (tilematrix + 1) ]).sum()
+ tile = (tilerow * cols) + tilecol + level_start
- mrf_path =
"/Users/kgill/data/ASCATB-L2-Coastal/ASCATB-L2-Coastal/MRF-GEO/2018/ASCATBL2Coastal_2018254_.mrf"
+ offset, size = self.__read_offset_and_size(tile)
+ tile_data = self.__read_tile_data(offset, size)
+ return tile_data
- reader = MrfReader(mrf_path)
\ No newline at end of file
diff --git a/analysis/webservice/algorithms/imaging/reprojection.py
b/analysis/webservice/algorithms/imaging/reprojection.py
index 7c22f2b..9c4abe2 100644
--- a/analysis/webservice/algorithms/imaging/reprojection.py
+++ b/analysis/webservice/algorithms/imaging/reprojection.py
@@ -18,76 +18,186 @@ import numpy as np, sys
from osgeo import gdal, osr
from osgeo.gdalconst import *
import tempfile
+from PIL import Image
gdal.AllRegister()
-def reproject_2d_matrix(from_data):
+class TargetProjection:
+ NORTH = "north"
+ SOUTH = "south"
- #sr.SetWellKnownGeogCS("EPSG:4326")
- output_raster = "foo.tif"
+latitude_origin_north = 52.6
+latitude_origin_south = -52.6
+dst_proj4_north = '+proj=stere +lat_0=90 +lat_ts=52.6 +lon_0=-45 +k=1 +x_0=0
+y_0=0 +datum=WGS84 +units=m +no_defs'
+dst_proj4_south = '+proj=stere +lat_0=-90 +lat_ts=-52.6 +lon_0=0 +k=1 +x_0=0
+y_0=0 +datum=WGS84 +units=m +no_defs'
+
+
+def create_global_dataset(data):
driver = gdal.GetDriverByName("MEM")
- dataset = driver.Create(output_raster, from_data.shape[1],
from_data.shape[0], 1, gdal.GDT_Float32)
- gt = (-180.0, (360.0 / float(from_data.shape[1])), 0, -90.0, 0, (180.0 /
float(from_data.shape[0])))
- dataset.SetGeoTransform(gt)
+ src_ds = driver.Create("", data.shape[1], data.shape[0], 1, gdal.GDT_Byte)
+ gt = (-180.0,
+ (360.0 / float(data.shape[1])),
+ 0,
+ 90.0,
+ 0,
+ -(180.0 / float(data.shape[0])))
+
+ sr = osr.SpatialReference()
+ sr.ImportFromEPSG(4326)
+
+ src_ds.SetProjection(sr.ExportToWkt())
+ src_ds.SetGeoTransform(gt)
+ src_ds.GetRasterBand(1).WriteArray(data)
- dataset.GetRasterBand(1).WriteArray(np.flip(from_data, axis=0))
+ return src_ds
+def create_north_dataset(data, latitude_cutoff=latitude_origin_north):
+
+ lat_cutoff_pixel = int(round((90 - latitude_cutoff) / 180.0 *
float(data.shape[0])))
+ data_cutoff = data[:lat_cutoff_pixel]
+ driver = gdal.GetDriverByName("MEM")
+ src_ds = driver.Create("", data_cutoff.shape[1], data_cutoff.shape[0], 1,
gdal.GDT_Byte)
+ gt = (-180.0,
+ (360.0 / float(data_cutoff.shape[1])),
+ 0,
+ 90.0,
+ 0,
+ -((90 - latitude_cutoff) / float(data_cutoff.shape[0])))
+ print gt
sr = osr.SpatialReference()
sr.ImportFromEPSG(4326)
- dataset.SetProjection(sr.ExportToWkt())
-
- outSpatialRef = osr.SpatialReference()
- outSpatialRef.ImportFromEPSG(3031)
-
- from_sr = sr
- to_sr = outSpatialRef
-
- pixel_spacing = 105000.
- tx = osr.CoordinateTransformation(from_sr, to_sr)
- geo_t = dataset.GetGeoTransform()
- print geo_t
- x_size = dataset.RasterXSize # Raster xsize
- y_size = dataset.RasterYSize # Raster ysize
- # Work out the boundaries of the new dataset in the target projection
- (ulx, uly, ulz) = tx.TransformPoint(geo_t[0], geo_t[3])
- (lrx, lry, lrz) = tx.TransformPoint(geo_t[0] + geo_t[1] * x_size, geo_t[3]
+ geo_t[5] * y_size)
-
- #dest_cols = int((lrx - ulx) / pixel_spacing)
- #dest_rows = int((uly - lry) / pixel_spacing)
- #print dest_cols, dest_rows
- #print dest_cols, dest_rows
- dest = driver.Create('', from_data.shape[1], from_data.shape[0], 1,
gdal.GDT_Float32)
- # Calculate the new geotransform
- new_geo = (ulx, pixel_spacing, geo_t[2], \
- uly, geo_t[4], -pixel_spacing)
- # Set the geotransform
- dest.SetGeoTransform(new_geo)
- dest.SetProjection(to_sr.ExportToWkt())
-
-
- res = gdal.ReprojectImage(dataset, dest,
- from_sr.ExportToWkt(), """PROJCS["WGS 84 / NSIDC
Sea Ice Polar Stereographic North",GEOGCS["WGS
84",DATUM["WGS_1984",SPHEROID["WGS
84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",70],PARAMETER[
[...]
- gdal.GRA_Bilinear)
-
- a = dest.GetRasterBand(1).ReadAsArray()
- return a
-
-def reproject_rgba_map(from_data):
- from_data[:, :, 0] = reproject_2d_matrix(from_data[:, :, 0])
- from_data[:, :, 1] = reproject_2d_matrix(from_data[:, :, 1])
- from_data[:, :, 2] = reproject_2d_matrix(from_data[:, :, 2])
- from_data[:, :, 3] = reproject_2d_matrix(from_data[:, :, 3])
+ src_ds.SetProjection(sr.ExportToWkt())
+ src_ds.SetGeoTransform(gt)
-if __name__ == "__main__":
+ src_ds.GetRasterBand(1).WriteArray(data_cutoff)
+
+ return src_ds
+
+
+def create_south_dataset(data, latitude_cutoff=latitude_origin_south):
+ lat_cutoff_pixel = int(round((90 - latitude_cutoff) / 180.0 *
float(data.shape[0])))
+
+ data_cutoff = data[lat_cutoff_pixel:]
+ driver = gdal.GetDriverByName("GTiff")
+ src_ds = driver.Create("test_south.tif", data_cutoff.shape[1],
data_cutoff.shape[0], 1, gdal.GDT_Byte)
+ gt = (-180.0,
+ (360.0 / float(data_cutoff.shape[1])),
+ 0,
+ latitude_cutoff,
+ 0,
+ -((latitude_cutoff - -90) / float(data_cutoff.shape[0])))
+
+ print gt
+ sr = osr.SpatialReference()
+ sr.ImportFromEPSG(4326)
+
+ src_ds.SetProjection(sr.ExportToWkt())
+ src_ds.SetGeoTransform(gt)
+
+ src_ds.GetRasterBand(1).WriteArray(data_cutoff)
+
+ return src_ds
+
+
+
+def reproject_2d_matrix(from_data, target_proj):
+ srs = osr.SpatialReference()
- data = np.zeros((512, 1024, 4))
- data[250:260, 0:1024, :] = 255
- data[0:512, 512:522, :] = 255
+ if target_proj == TargetProjection.NORTH:
+ src_ds = create_north_dataset(from_data, latitude_origin_north)
+ srs.ImportFromProj4(dst_proj4_north)
+ dst_wkt = srs.ExportToWkt()
+ elif target_proj == TargetProjection.SOUTH:
+ src_ds = create_south_dataset(from_data, latitude_origin_south)
+ srs.ImportFromProj4(dst_proj4_south)
+ dst_wkt = srs.ExportToWkt()
+ else:
+ raise Exception("Unsupported target projection")
- reproject_rgba_map(data)
+ error_threshold = 0.125 # error threshold --> use same value as in
gdalwarp
+ resampling = gdal.GRA_NearestNeighbour
+
+ tmp_ds = gdal.AutoCreateWarpedVRT(src_ds,
+ None, # src_wkt : left to default value
--> will use the one from source
+ dst_wkt,
+ resampling,
+ error_threshold)
+
+
+ dst_xsize = tmp_ds.RasterXSize
+ dst_ysize = tmp_ds.RasterYSize
+ dst_gt = tmp_ds.GetGeoTransform()
+
+ tmp_ds = None
+
+ dst_ds = gdal.GetDriverByName('MEM').Create('', dst_xsize, dst_ysize,
src_ds.RasterCount)
+ dst_ds.SetProjection(dst_wkt)
+ dst_ds.SetGeoTransform(dst_gt)
+ dst_ds.GetRasterBand(1).SetNoDataValue(0)
+
+ res = gdal.ReprojectImage(src_ds,
+ dst_ds,
+ None, # src_wkt : left to default value -->
will use the one from source
+ None, # dst_wkt : left to default value -->
will use the one from destination
+ resampling,
+ 0, # WarpMemoryLimit : left to default value
+ error_threshold,
+ None, # Progress callback : could be left to
None or unspecified for silent progress
+ None,
+ ("SOURCE_EXTRA=125",)) # Progress callback
user data)
+
+ assert res == 0, 'Error in ReprojectImage'
+
+ return dst_ds.GetRasterBand(1).ReadAsArray()
+
+def reproject_rgba_map(from_data, target_proj):
+ data_0 = reproject_2d_matrix(from_data[:, :, 0], target_proj)
+ data_1 = reproject_2d_matrix(from_data[:, :, 1], target_proj)
+ data_2 = reproject_2d_matrix(from_data[:, :, 2], target_proj)
+ data_3 = reproject_2d_matrix(from_data[:, :, 3], target_proj)
+
+ dest_data = np.zeros((data_0.shape[0], data_0.shape[1], 4))
+ dest_data[:, :, 0] = data_0
+ dest_data[:, :, 1] = data_1
+ dest_data[:, :, 2] = data_2
+ dest_data[:, :, 3] = data_3
+
+
+ return dest_data
+
+
+
+
+def create_test_data(h=512, w=1024):
+ data = np.zeros((h, w, 4))
+
+ h_mid = h / 2
+ w_mid = w / 2
+ data[h_mid-10:h_mid+10, 0:w, 0:2] = 255
+ data[0:h, w_mid-5:w_mid+5, 0:2] = 255
+ data[:, :, 3] = 255
+ return data
+
+
+def read_test_data():
+ path = "test_pre.png"
+ image = Image.open(path)
+ im = np.copy(np.asarray(image, dtype=np.uint8))
+ return im
+
+
+if __name__ == "__main__":
+ data = read_test_data()
+ out_data = reproject_rgba_map(data, TargetProjection.SOUTH)
+ im = Image.fromarray(np.asarray(out_data, dtype=np.uint8))
+ im.save("test_warp.png")
+ #im = Image.fromarray(np.asarray(out_data, dtype=np.uint8))
+ #im.save("test_warp.tif")
+
+ #reproject_rgba_map(data)
diff --git a/data-access/requirements.txt b/data-access/requirements.txt
index 1901db2..651f776 100644
--- a/data-access/requirements.txt
+++ b/data-access/requirements.txt
@@ -14,7 +14,7 @@ futures==3.1.1
ipython==5.3.0
ipython-genutils==0.2.0
jmespath==0.9.3
-nexusproto==1.0.0-SNAPSHOT
+nexusproto==1.0.0
numpy==1.11.1
pathlib2==2.2.1
pexpect==4.2.1
diff --git a/data-access/setup.py b/data-access/setup.py
index 532dede..031ccb0 100644
--- a/data-access/setup.py
+++ b/data-access/setup.py
@@ -16,7 +16,7 @@
import setuptools
from Cython.Build import cythonize
-__version__ = '0.32'
+__version__ = '0.33'
setuptools.setup(
name="nexus-data-access",
@@ -36,7 +36,7 @@ setuptools.setup(
'cassandra-driver==3.5.0',
'pysolr==3.7.0',
'requests',
- 'nexusproto==1.0.0-SNAPSHOT',
+ 'nexusproto==1.0.0',
'shapely'
],
diff --git a/docker/nexus-webapp/Dockerfile b/docker/nexus-webapp/Dockerfile
index 3b2d5f9..94ea293 100644
--- a/docker/nexus-webapp/Dockerfile
+++ b/docker/nexus-webapp/Dockerfile
@@ -82,6 +82,7 @@ ARG
APACHE_NEXUSPROTO=https://github.com/apache/incubator-sdap-nexusproto.git
ARG APACHE_NEXUSPROTO_BRANCH=master
ARG APACHE_NEXUS=https://github.com/apache/incubator-sdap-nexus.git
ARG APACHE_NEXUS_BRANCH=master
+ARG REBUILD_CODE=1
RUN /tmp/install_nexusproto.sh $APACHE_NEXUSPROTO $APACHE_NEXUSPROTO_BRANCH &&
\
/tmp/install_nexus.sh $APACHE_NEXUS $APACHE_NEXUS_BRANCH $NEXUS_SRC
diff --git a/docker/spark-mesos-base/Dockerfile
b/docker/spark-mesos-base/Dockerfile
index 784f56c..4c27bd6 100644
--- a/docker/spark-mesos-base/Dockerfile
+++ b/docker/spark-mesos-base/Dockerfile
@@ -125,6 +125,7 @@ RUN yum install -y mesa-libGL.x86_64
# Install nexusproto
ARG APACHE_NEXUSPROTO=https://github.com/apache/incubator-sdap-nexusproto.git
ARG APACHE_NEXUSPROTO_BRANCH=master
+ARG REBUILD_CODE=1
COPY install_nexusproto.sh ./install_nexusproto.sh
RUN ./install_nexusproto.sh $APACHE_NEXUSPROTO $APACHE_NEXUSPROTO_BRANCH
