[
https://issues.apache.org/jira/browse/SDAP-149?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=16675495#comment-16675495
]
ASF GitHub Bot commented on SDAP-149:
-------------------------------------
fgreg closed pull request #49: SDAP-149 Generate images directly from data
URL: https://github.com/apache/incubator-sdap-nexus/pull/49
This is a PR merged from a forked repository.
As GitHub hides the original diff on merge, it is displayed below for
the sake of provenance:
As this is a foreign pull request (from a fork), the diff is supplied
below (as it won't show otherwise due to GitHub magic):
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 @@ def _create_values(min_value, max_value, num_values=255):
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 @@ def calc(self, computeOptions, **args):
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 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 @@ def __init__(self):
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 @@ def _build_layer_spec(self, ds):
}
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("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 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 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 @@ def __init__(self, mrf_meta_path):
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 @@ def __init__(self, mrf_path):
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 @@
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["central_meridian",-45],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],AUTHORITY["EPSG","3413"],AXIS["X",UNKNOWN],AXIS["Y",UNKNOWN]]""",
- 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 @@
'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
----------------------------------------------------------------
This is an automated message from the Apache Git Service.
To respond to the message, please log on GitHub and use the
URL above to go to the specific comment.
For queries about this service, please contact Infrastructure at:
[email protected]
> Generate images directly from data
> ----------------------------------
>
> Key: SDAP-149
> URL: https://issues.apache.org/jira/browse/SDAP-149
> Project: Apache Science Data Analytics Platform
> Issue Type: New Feature
> Components: nexus
> Reporter: Frank Greguska
> Assignee: Kevin M. Gill
> Priority: Major
>
> Develop the capability to generate images directly from Nexus data.
> Create services:
> /map - Generate colorized map given a dataset, time, and bounding box
> /wmts - Emulate a simple WMTS service
> - EPSG4326 (equirectangular) support initially
> - In Memory & S3 caching support
> - Support for static, pregenerated layers (tiles, MRF)
> /colortables - Lists supported color tables and their specifications as
> CMC-compatible JSON
> /layers - Creates a CMC-compatible layer configuration file
> /thumbnail - Creates a small thumbnail for a layer
--
This message was sent by Atlassian JIRA
(v7.6.3#76005)
