Needs to be able to autmoatically handle DEM boundaries

[tjlang]

For example, you get a half view if you place a radar near your DEM borders (e.g., 40 N).

[nguy]

@tjlang I didn't consider that. I have an idea to deal with this and will try to get this implemented soon.

[nguy]

So the wradlib port is taking advantage of the gdal package. There appears to be a way to merge files (at least tiff images) into a mosaic. I will try to start playing around with this.

[kmuehlbauer]

Just as a side note: GDAL and its several (standalone-) utility functions can make use of many raster file formats, depending on compiled libraries.

Using the GTOPO30 data for instance, you could cut germany out of the W020N90 tile by:

gdalwarp -te 5.5 47. 15.5 55. -of EHdr w020n90.dem germany.dem

with -te beeing the bounding box and -of the destination raster format. It creates a germany.dem and the according germany.hdr.

If you need two tiles (if the radar is located near a tile border) you just have to give this second tile together with the first tile as input parameter. And, if no input file format is specified, gdal tries to guess the format.

[nguy]

@kmuehlbauer Thanks! I hadn't looked into the functions yet and that's great that many formats can be managed. Since it's a limited data set in only using GTOPO30 at the moment, I hope to be able to automate the search.

[kmuehlbauer]

@nguy If you allow the use of subprocess here is an approach using gdal vrt:

import subprocess
import numpy as np
import matplotlib as mpl
import glob
from osgeo import gdal

# precalculated radar bounding box
corners = np.array([[lonmin, latmin], [lonmin, latmax], [lonmax,latmax], [lonmax, latmin]])

# get all .dem files except antarcps.dem
path_to_gtopo = '/work/radproc/data/radar/gtopo30/data/'
flist = sorted(glob.glob(path_to_gtopo + '[e,w]*.dem'))
found = []

# iterate over dem files
for file in flist:
    ds = gdal.Open(filename)
    size = (ds.RasterXSize, ds.RasterYSize)
    geot = ds.GetGeoTransform()
    # get upper left corner and scale
    ulx = geot[0]
    uly = geot[3]
    scale = (geot[1], geot[5])
    # calculate lower right corner 
    lrx = ulx + size[0] * scale[0]
    lry = uly + size[1] * scale[1]

    # create bounding box and matplotlib Path
    bbox = [[ulx, lry], [ulx, uly], [lrx, uly], [lrx, lry]]
    bbpath = mpl.path.Path(bbox)

    # Path contains radar corner points ?
    pinside = bbpath.contains_points(corners)

    # if any point is inside append file to found
    if np.any(pinside):
        found.append(file)
            # if all points are inside, we can break 
            if np.all(pinside):
                break
print(found)

# setup gdalbuildvrt command to create a interim vrt file clipped with the radar bounding box 
command=["gdalbuildvrt",  '-overwrite', '-te', str(lonmin-scale[0]),  str(latmin-scale[1]),
             str(lonmax+scale[0]), str(latmax+scale[1]), path_to_gtopo+ "interim.vrt"]
# add needed files to vrt command
for file in found:
    command.append(file)
# run gdalbuildvirt
subprocess.call(command)
# use the created file just as you would use any other gdal raster file
rasterfile = path_to_gtopo+'interim.vrt'

I have just tested it, it performs fine. And it will find every tile which has at least one radar bounding box corner point inside. If there is a radar located near a gtopo30 grid point, it will find all four needed dem-files. The advantage using vrt is no intermediate data, and no fuzz with gdal-python internals.

[kmuehlbauer]

Downside is the need of gdal as a python dependency. But gdal utils need to be installed anyway.

To overcome this you can use the .hdr files to extract the needed information to create the bounding box etc (just as you do now). This actually works, but you have to change some signs in the calculation of the dem bounding box corners in respect to the gdal approach above.

And then you would have to call e.g gdalwarp using subprocess:

gdalwarp -of Ehdr interim.vrt interim.dem

Then you can read interim.dem just as you are reading the gtopo30 files.

[kmuehlbauer]

see pull request #7

[qxang]

If the dem is very big,It will produce problems.

[nguy]

@shanyuk Can you elaborate on what kind of problems?