0

I have an application that uses raster elevation data extensively in its computations. I would like to have a tool that allows users to define the minimum and maximum latitude-longitude coordinates (as a bounding box) and download publicly available elevation data from a web server as a single TIFF file or some other raster format that supports elevations. I don't know whether this would be a Web Map Service or a Web Coverage Service. I have searched online for something like this and have not been able to find anything. Of course, I know about the National Map for US data. I am also familiar with Google Earth and have used NASA World Wind in the past. I am also familiar with ESRI's World Elevation Services.

Does anybody have details of a web server that does something like this? Free would be nice, but maybe there's something available for a fee based on how much the server is used. If not, any suggestions of how I could set up a server that does something like this myself if I download all the elevation data I need?

  • GeoServer installs by default one demo layer that is DEM. DescribeCoverage will be found at localhost:8080/geoserver/…. Acquire your DEM files, make an image mosaic, publish and you've done. – user30184 May 8 '17 at 21:20
  • Maybe mapzen.com/documentation/terrain-tiles see elevatr R pkg for an interface, though that needs a wrapper to do extent – mdsumner May 8 '17 at 23:44
  • @mdsumner: I decided that mapzen.com option might be the best way for us to go. It looks like it's unlimited freely available terrain tiles for now according to their web site link. I'll try to report how this works for us after we're done with implementation. – Burrow May 9 '17 at 21:29
  • It's cool to finally have this available. I haven't explored it deeply, a big missing part of elevatr is to be able specify an extent/grain (i.e. give it a raster object to populate), you have to figure out that interactively from zoom and the default dimensions returned. – mdsumner May 9 '17 at 23:21
0

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})

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.