Skip to content

Commit

Permalink
Additional fixes re: ONEARTH-516
Browse files Browse the repository at this point in the history
  • Loading branch information
@jdrodjpl
jdrodjpl committed Aug 31, 2018
1 parent 9e6b7ab commit 5e12ef2
Show file tree
Hide file tree
Showing 3 changed files with 105 additions and 25 deletions.
1 change: 1 addition & 0 deletions doc/config_reproject.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,3 +68,4 @@ The reprojected layer config tool scrapes data from the GetCapabilities file of

`<IncludeLayer>` (optional, can be multiple) -- Include one of these elements for each layer that you want to exclusively configure. **If this element is found, only the layers specified in each element will be configured, and all others will be ignored.**

`<TargetEpsgCode>` (generally required) -- Sets the output projection for the reprojected layers. Currently, only compatible with EPSG:4326 and EPSG:3857. If this parameter is omitted, it can be used to generate mapfiles that use a tile source but do not reproject the tiles.
2 changes: 2 additions & 0 deletions src/layer_config/layers/layer_configuration_file_reproject.xml.sample
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,4 +2,6 @@
<SrcWMTSGetCapabilitiesURI>https://gibs.earthdata.nasa.gov/wmts/epsg4326/best/1.0.0/WMTSCapabilities.xml</SrcWMTSGetCapabilitiesURI>
<SrcLocationRewrite internal="/outside" external="https://gibs.earthdata.nasa.gov" />
<EnvironmentConfig>/etc/onearth/config/conf/environment_webmercator.xml</EnvironmentConfig>
<IncludeLayer>AMSR2_Snow_Water_Equivalent</IncludeLayer>
<TargetEpsgCode>3857</TargetEpsgCode>
</ReprojectLayerConfig>
127 changes: 102 additions & 25 deletions src/scripts/oe_configure_reproject_layer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,14 +46,15 @@
import re
import shutil
from optparse import OptionParser
from osgeo import osr
from osgeo import osr, ogr
import png
from io import BytesIO
from time import asctime
import cgi
import hashlib
from decimal import Decimal
import copy
import json

reload(sys)
sys.setdefaultencoding('utf8')
Expand All @@ -78,13 +79,13 @@
STATUS ON
METADATA
"wms_title" "{layer_title}"
"wms_srs" "EPSG:3857"
"wms_extent" "-180 -85.0511 180 85.0511"
"wms_srs" "EPSG:{target_epsg}"
"wms_extent" "{target_bbox}"
{dimension_info}
END
DATA '{data_xml}'
PROJECTION
"init=epsg:4326"
"init=epsg:{src_epsg}"
END
{validation_info}
END
Expand Down Expand Up @@ -114,32 +115,92 @@ def bulk_replace(source_str, replace_list):
return out_str


def make_gdal_tms_xml(layer, bands):
def get_epsg_code_for_proj_string(proj_string):
# For some reason OSR can't parse this version of the 3857 CRS.
if proj_string == 'urn:ogc:def:crs:EPSG:6.18:3:3857':
return '3857'
proj = osr.SpatialReference()
proj.SetFromUserInput(proj_string)
return proj.GetAuthorityCode(None)


def get_bbox_for_proj_string(proj_string, use_oe_tms=False, get_in_map_units=False):
epsg_code = get_epsg_code_for_proj_string(proj_string)

if epsg_code == '4326':
if use_oe_tms:
return [-180, -198, 396.0, 90]
else:
return [-180, -90, 180, 90]
bbox = get_proj_bbox(epsg_code)
if not get_in_map_units:
return bbox
src_proj = osr.SpatialReference()
src_proj.SetFromUserInput('EPSG:4326')
target_proj = osr.SpatialReference()
target_proj.ImportFromEPSG(int(epsg_code))
transform = osr.CoordinateTransformation(src_proj, target_proj)

point_coords = [(bbox[1], bbox[0]), (bbox[3], bbox[2])]
new_bbox = []
for coords in point_coords:
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(coords[0], coords[1])
point.Transform(transform)
new_bbox.append(point.GetX())
new_bbox.append(point.GetY())
return [new_bbox[0], new_bbox[3], new_bbox[1], new_bbox[2]]


def get_proj_bbox(epsg_code):
bbox = [-180, -90, 180, 90]
uri = 'http://epsg.io/?q={0}&format=json'.format(epsg_code)
try:
r = requests.get(uri)
if r.status_code == 200:
res = json.loads(r.content)
bbox = res['results'][0]['bbox']
except:
pass
return bbox


def make_gdal_tms_xml(layer, bands, src_epsg):
tms = layer.find('{*}TileMatrixSetLink').findtext('{*}TileMatrixSet')

bbox = map(str, get_bbox_for_proj_string(
'EPSG:' + src_epsg, use_oe_tms=True, get_in_map_units=src_epsg != '4326'))

resource_url = layer.find('{*}ResourceURL')
template_string = resource_url.get('template')

out_root = etree.Element('GDAL_WMS')

service_element = etree.SubElement(out_root, 'Service')
service_element.set('name', 'TMS')
etree.SubElement(service_element, 'ServerUrl').text = template_string.replace(
'{TileMatrixSet}', tms).replace('{Time}', '%time%').replace('{TileMatrix}', '${z}').replace('{TileRow}', '${y}').replace('{TileCol}', '${x}')
etree.SubElement(service_element, 'ServerUrl').text = bulk_replace(template_string, [
('{TileMatrixSet}', tms), ('{Time}', '%time%'), ('{TileMatrix}',
'${z}'), ('{TileRow}', '${y}'), ('{TileCol}', '${x}')])

data_window_element = etree.SubElement(out_root, 'DataWindow')
etree.SubElement(data_window_element, 'UpperLeftX').text = '-180.0'
etree.SubElement(data_window_element, 'UpperLeftY').text = '90'
etree.SubElement(data_window_element, 'LowerRightX').text = '396.0'
etree.SubElement(data_window_element, 'LowerRightY').text = '-198'
etree.SubElement(data_window_element, 'TileLevel').text = TILE_LEVELS[tms]
etree.SubElement(data_window_element, 'TileCountX').text = '2'
etree.SubElement(data_window_element, 'UpperLeftX').text = bbox[0]
etree.SubElement(data_window_element, 'UpperLeftY').text = bbox[3]
etree.SubElement(data_window_element, 'LowerRightX').text = bbox[2]
etree.SubElement(data_window_element, 'LowerRightY').text = bbox[1]
if src_epsg == '3857':
tile_levels = tms.split('GoogleMapsCompatible_Level')[1]
else:
tile_levels = TILE_LEVELS[tms]
etree.SubElement(data_window_element, 'TileLevel').text = tile_levels
etree.SubElement(data_window_element,
'TileCountX').text = '2' if src_epsg == '4326' else '1'
etree.SubElement(data_window_element, 'TileCountY').text = '1'
etree.SubElement(data_window_element, 'YOrigin').text = 'top'

etree.SubElement(out_root, 'Projection').text = 'EPSG:4326'
etree.SubElement(out_root, 'BlockSizeX').text = '512'
etree.SubElement(out_root, 'BlockSizeY').text = '512'
etree.SubElement(out_root, 'Projection').text = 'EPSG:' + src_epsg
tile_size = '512' if src_epsg == '4326' else '256'
etree.SubElement(out_root, 'BlockSizeX').text = tile_size
etree.SubElement(out_root, 'BlockSizeY').text = tile_size
etree.SubElement(out_root, 'BandsCount').text = str(bands)

etree.SubElement(out_root, 'Cache')
Expand Down Expand Up @@ -202,6 +263,9 @@ def build_reproject_configs(layer_config_path, tilematrixsets_config_path, wmts=
name.text for name in config_xml.findall('ExcludeLayer')]
layer_include_list = [
name.text for name in config_xml.findall('IncludeLayer')]
target_epsg = config_xml.findtext('TargetEpsgCode')
if target_epsg and target_epsg.startswith('EPSG:'):
target_epsg = target_epsg.split(':')[1]

# Parse the tilematrixsets definition file
try:
Expand Down Expand Up @@ -494,7 +558,7 @@ def build_reproject_configs(layer_config_path, tilematrixsets_config_path, wmts=
dest_tilematrixset_name = dest_tilematrixset.findtext(
ows + 'Identifier')
dest_scale_denominator = float(sort_tilematrixset(dest_tilematrixset)[
0].findtext('{*}ScaleDenominator'))
0].findtext('{*}ScaleDenominator'))
dest_width = dest_height = int(
round(2 * math.pi * EARTH_RADIUS / (dest_scale_denominator * 0.28E-3)))
dest_levels = len(dest_tilematrixset.findall('{*}TileMatrix'))
Expand All @@ -505,7 +569,7 @@ def build_reproject_configs(layer_config_path, tilematrixsets_config_path, wmts=
dest_top_left_corner = [Decimal(value) for value in sort_tilematrixset(
dest_tilematrixset)[0].findtext('{*}TopLeftCorner').split(' ')]
dest_bbox = '{0},{1},{2},{3}'.format(dest_top_left_corner[0], dest_top_left_corner[
0], -dest_top_left_corner[0], -dest_top_left_corner[0])
0], -dest_top_left_corner[0], -dest_top_left_corner[0])
dest_resource_url_elem = layer.find('{*}ResourceURL')
dest_url = dest_resource_url_elem.get('template')
for location in src_locations:
Expand Down Expand Up @@ -861,9 +925,9 @@ def build_reproject_configs(layer_config_path, tilematrixsets_config_path, wmts=
patterns = ''
for tilematrix in sorted(out_tilematrixsets[0].findall('{*}TileMatrix'), key=lambda matrix: float(matrix.findtext('{*}ScaleDenominator'))):
resx = (dest_top_left_corner[
1] - dest_top_left_corner[0]) / Decimal(tilematrix.findtext('MatrixWidth'))
1] - dest_top_left_corner[0]) / Decimal(tilematrix.findtext('MatrixWidth'))
resy = (dest_top_left_corner[
1] - dest_top_left_corner[0]) / Decimal(tilematrix.findtext('MatrixHeight'))
1] - dest_top_left_corner[0]) / Decimal(tilematrix.findtext('MatrixHeight'))
local_xmax = dest_top_left_corner[0] + resx
local_xmax_str = str(local_xmax) if local_xmax else '0'
local_ymax = dest_top_left_corner[1] - resy
Expand All @@ -890,14 +954,27 @@ def build_reproject_configs(layer_config_path, tilematrixsets_config_path, wmts=
if not static:
default_datetime = dest_dim_elem.findtext('{*}Default')
dimension_info = bulk_replace(DIMENSION_TEMPLATE, [('{periods}', dest_dim_elem.findtext("{*}Value")),
('{default}', default_datetime)])
validation_info = VALIDATION_TEMPLATE.replace('{default}', default_datetime)
('{default}', default_datetime)])
validation_info = VALIDATION_TEMPLATE.replace(
'{default}', default_datetime)

# Mapserver automatically converts to RGBA and works better if we specify that for png layers
# Mapserver automatically converts to RGBA and works better if we
# specify that for png layers
mapserver_bands = 4 if 'image/png' in src_format else 3

src_crs = layer_tilematrixsets[0].findtext('{*}SupportedCRS')
src_epsg = get_epsg_code_for_proj_string(src_crs)

if not target_epsg:
target_epsg = src_epsg
target_bbox = map(
str, get_bbox_for_proj_string('EPSG:' + target_epsg, get_in_map_units=src_epsg != '4326'))
target_bbox = [target_bbox[1], target_bbox[
0], target_bbox[3], target_bbox[2]]

mapfile_snippet = bulk_replace(
MAPFILE_TEMPLATE, [('{layer_name}', identifier), ('{data_xml}', make_gdal_tms_xml(src_layer, mapserver_bands)), ('{layer_title}', cgi.escape(src_title)),
('{dimension_info}', dimension_info), ('{validation_info}', validation_info)])
MAPFILE_TEMPLATE, [('{layer_name}', identifier), ('{data_xml}', make_gdal_tms_xml(src_layer, mapserver_bands, src_epsg)), ('{layer_title}', cgi.escape(src_title)),
('{dimension_info}', dimension_info), ('{validation_info}', validation_info), ('{src_epsg}', src_epsg), ('{target_epsg}', target_epsg), ('{target_bbox}', ', '.join(target_bbox))])

mapfile_name = os.path.join(
mapfile_staging_location, identifier + '.map')
Expand Down

0 comments on commit 5e12ef2

Please sign in to comment.