OpenTopography SRTM VRT Parser

vrt2pgsql.py

import sys
import xml.etree.ElementTree as ET
import argparse

######################################################################

parser = argparse.ArgumentParser()
parser.add_argument('-f', '--file',
                    required=True,
                    dest='file',
                    help='VRT file to read',
                    default='SRTM_GL1_srtm.vrt'
                    )
parser.add_argument('-u', '--url',
                    dest='url',
                    help='Base URL of VRT file',
                    default='https://opentopography.s3.sdsc.edu/raster/SRTM_GL1'
                    )

args = parser.parse_args()

######################################################################

def img_to_geo(img_x, img_y, gt):
    # offset_x, scale_x, skew_x
    geo_x = gt[0] + img_x * gt[1] + img_y * gt[2]
    # offset_y, scale_y, skew_y
    geo_y = gt[3] + img_y * gt[5] + img_x * gt[4] 
    return geo_x, geo_y


root_node = ET.parse(args.file).getroot()

raster_xsize = root_node.attrib['rasterXSize']
raster_ysize = root_node.attrib['rasterYSize']

srs = root_node.find('SRS').text

# GeoTransform is a 6-member affine to convert from image coordinates
# to geographic coordinates
geotransform = root_node.find('GeoTransform').text
gt = [float(s.strip()) for s in geotransform.split(',')]

#print("CREATE TABLE filebounds (fileurl text PRIMARY KEY, geom geometry(polygon, 4326)")

first = True
for tag in root_node.findall('VRTRasterBand/ComplexSource'):

    # <SourceFilename relativeToVRT="1">SRTM_GL1_srtm/S06E117.tif</SourceFilename>
    # <SourceBand>1</SourceBand>
    # <SourceProperties RasterXSize="3601" RasterYSize="3601" DataType="Int16" BlockXSize="512" BlockYSize="512" />
    # <SrcRect xOff="0" yOff="0" xSize="3601" ySize="3601" />
    # <DstRect xOff="1069200" yOff="234000" xSize="3601" ySize="3601" />
    # <NODATA>-32768</NODATA>
    
    source_filename = tag.find('SourceFilename')
    source_band = tag.find('SourceBand')
    source_properties = tag.find('SourceProperties')
    src_rect = tag.find('SrcRect')
    nodata = int(tag.find('NODATA').text)

    # Add URL prefix as necessary
    fileurl = source_filename.text
    if source_filename.attrib['relativeToVRT'] == '1':
        fileurl = args.url + '/' + source_filename.text

    dst_rect = tag.find('DstRect')
    ix = int(dst_rect.attrib['xOff'])
    iy = int(dst_rect.attrib['yOff'])
    (gx, gy) = img_to_geo(ix, iy, gt)

    # Fill in metadata about this file for generating SQL and shell commands
    src = {
        'width': int(dst_rect.attrib['xSize']),
        'height': int(dst_rect.attrib['ySize']),
        'ulx': gx,
        'uly': gy,
        'scalex': gt[1],
        'scaley': gt[5],
        'skewx': gt[2],
        'skewy': gt[4],
        'nodata': nodata,
        'fileurl': fileurl,
        'tbl': 'dem',
        'database_url': 'dbname=opentopo'
    }

    src['lrx'] = src['ulx'] + src['width'] * src['scalex']
    src['lry'] = src['uly'] + src['height'] * src['scaley']

    # Old entries for other loading tricks

    # Get the bounds of each file and insert them into a query table
    # src['mktile'] = """
    #     INSERT INTO filebounds (fileurl, geom)
    #     VALUES ('%(fileurl)s', ST_SetSRID(ST_MakeEnvelope(%(ulx)s, %(lry)s, %(lrx)s, %(uly)s), 4326));
    #     """ % src

    # src['mkraster'] = """
    #     ST_MakeEmptyRaster(width => %(width)s, height => %(height)s, 
    #         upperleftx => %(ulx)s, upperlefty => %(uly)s,
    #         scalex => %(scalex)s, scaley => %(scaley)s,
    #         skewx => %(skewx)s, skewy => %(skewy)s,
    #         srid => 4326)""" % src

    # src['mkband'] = """
    #     ST_AddBand(%(mkraster)s, 1, '%(fileurl)s'::text, ARRAY[1::int4], %(nodata)s)""" % src

    # Create table on first run, append after that
    if first:
        src['tblflag'] = '-c'
        first = False
    else:
        src['tblflag'] = '-a'

    src['sh'] = "raster2pgsql %(tblflag)s -F -Y -k -N %(nodata)s -s 4326 -t 512x512 -I -R /vsicurl/%(fileurl)s %(tbl)s | psql -d %(database_url)s" % src

    print(src['sh'])