I think GeoServer would have worked, but I decided to go with the Mapzen Terrain Tiles since Amazon maintains the servers and it's free of charge. I originally implemented code to download the GeoTIFF tiles, but I can only get elevations to the nearest meter with this format since the elevations are stored as 16-bit integers. This did not provide the accuracy we needed. So I decided to go with the terrarium files, which are PNG files that can be converted from RGB values to elevation values at each pixel. These files work well and provide the accuracy we need. We download the tiles given a bounding box and convert each tile to a GeoTIFF file using GDAL. We then merge and crop the tiles to the bounding box using the gdal_merge.py script. Most of the code that does this can be found by searching around on the Mapzen site, but I will post the important parts of the Python code one can use to do the conversion here (I converted most of this to C++ for my application, but we don't use GDAL for everything so some of the C++ calls might not make sense. That's why I'm posting the Python):
def array2raster(raster_fname, raster_origin, pixel_width, pixel_height, array):
import gdal
import osr
cols = array.shape[1]
rows = array.shape[0]
origin_x = raster_origin[0]
origin_y = raster_origin[1]
driver = gdal.GetDriverByName('GTiff')
out_raster = driver.Create(raster_fname, cols, rows, 1, gdal.GDT_Float32)
out_raster.SetGeoTransform((origin_x, pixel_width, 0, origin_y, 0, pixel_height))
out_band = out_raster.GetRasterBand(1)
out_band.WriteArray(array)
out_band.SetNoDataValue(-9999.0)
out_raster_srs = osr.SpatialReference()
out_raster_srs.ImportFromEPSG(4326)
out_raster.SetProjection(out_raster_srs.ExportToWkt())
out_band.FlushCache()
def num2deg(x_tile, y_tile, zoom):
import math
n = 2.0 ** zoom
lon_deg = x_tile / n * 360.0 - 180.0
lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * y_tile / n)))
lat_deg = math.degrees(lat_rad)
return lon_deg, lat_deg
def deg2num(lon_deg, lat_deg, zoom):
import math
lat_rad = math.radians(lat_deg)
n = 2.0 ** zoom
x_tile = int((lon_deg + 180.0) / 360.0 * n)
y_tile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
return x_tile, y_tile
def test_mapzen(request):
if request.is_ajax() and request.method == 'GET':
from requests import get
from PIL import Image
import numpy as np
_image_width = 256
_image_height = 256
api_key = 'mapzen-avAZED5'
url_template = 'https://tile.mapzen.com/mapzen/terrain/v1/terrarium/{z}/{x}/{y}.png?api_key={api_key}'
params = request.GET
extent = loads(params.get('extent'))
zoom = int(params.get('zoom'))
north_edge = extent[3]
west_edge = extent[0]
south_edge = extent[1]
east_edge = extent[2]
left_tile, top_tile = deg2num(west_edge, north_edge, zoom)
right_tile, bottom_tile = deg2num(east_edge, south_edge, zoom)
pixel_width = None
pixel_height = None
# Get user workspace file paths
user_workspace = EpanetApp.get_user_workspace(request.user)
test_workspace_path = os.path.join(user_workspace.path, 'test')
if not os.path.exists(test_workspace_path):
os.makedirs(test_workspace_path)
tiffs_array = []
for x_tile in range(left_tile, right_tile + 1):
for y_tile in range(top_tile, bottom_tile + 1):
url = url_template.format(z=zoom, x=x_tile, y=y_tile, api_key=api_key)
png_name = '{z}_{x}_{y}.{ext}'.format(z=zoom, x=x_tile, y=y_tile, ext='png')
r = get(url)
png_path = os.path.join(test_workspace_path, png_name)
with open(png_path, 'w+') as f:
for chunk in r.iter_content(1024):
f.write(chunk)
raster_array = []
image = Image.open(png_path)
for y_pix in range(0, _image_height):
row_array = []
for x_pix in range(0, _image_width):
r, g, b = image.getpixel((x_pix, y_pix))
elevation = (r * 256.0 + g + b / 256.0) - 32768.0
row_array.append(elevation)
raster_array.append(row_array)
tif_name = '{z}_{x}_{y}.{ext}'.format(z=zoom, x=x_tile, y=y_tile, ext='tif')
tif_path = os.path.join(test_workspace_path, tif_name)
tif_origin = num2deg(x_tile, y_tile, zoom)
if not pixel_width:
next_x_tif_origin = num2deg(x_tile + 1, y_tile, zoom)
pixel_width = (next_x_tif_origin[0] - tif_origin[0]) / _image_width
if not pixel_height:
next_y_tif_origin = num2deg(x_tile, y_tile + 1, zoom)
pixel_height = (next_y_tif_origin[1] - tif_origin[1]) / _image_height
numpy_raster_array = np.array(raster_array)
array2raster(tif_path, tif_origin, pixel_width, pixel_height, numpy_raster_array)
tiffs_array.append(tif_path)
merged_tif_path = os.path.join(user_workspace.path, 'merged.tif')
process = Popen(['gdal_merge.py', '-o', merged_tif_path, '-ul_lr', str(west_edge), str(north_edge),
str(east_edge), str(south_edge), '-n', '-9999.0'] + tiffs_array)
process.communicate()
rmtree(test_workspace_path)
return JsonResponse({'success': True})
return JsonResponse({'success': False})
