pramsey · July 9, 2021 22:56
diff --git a/vrt2pgsql.py b/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'])
	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'])