Skip to content

Instantly share code, notes, and snippets.

@stepankuzmin
Created October 2, 2013 11:44
Show Gist options
  • Save stepankuzmin/6792465 to your computer and use it in GitHub Desktop.
Save stepankuzmin/6792465 to your computer and use it in GitHub Desktop.
See http://erouault.blogspot.ru/2011/12/seamless-access-to-remote-global-multi.html Don't know when this happened, but the machine 'igskmncngs506' is now 'gmted'.
#!/usr/bin/env python
###############################################################################
# $Id$
#
# Purpose: Builds global VRT for GMTED2010 data.
# Author: Even Rouault, <even dot rouault at mines dash paris dot org>
#
###############################################################################
# Copyright (c) 2011, Even Rouault, <even dot rouault at mines dash paris dot org>
#
# Permission is hereby granted, free of charge, to any person obtaining a
# copy of this software and associated documentation files (the "Software"),
# to deal in the Software without restriction, including without limitation
# the rights to use, copy, modify, merge, publish, distribute, sublicense,
# and/or sell copies of the Software, and to permit persons to whom the
# Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included
# in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
# DEALINGS IN THE SOFTWARE.
###############################################################################
from osgeo import gdal
import os
def build_src_filename(res, typ, lat_hemisphere, abs_lat, lon_hemisphere, abs_lon):
return '/vsicurl/http://gmted.cr.usgs.gov/gmted/Global_tiles_GMTED/%sdarcsec/%s/%s%03d/%02d%s%03d%s_20101117_gmted_%s%s.tif' % (res, typ, lon_hemisphere, abs_lon, abs_lat, lat_hemisphere, abs_lon, lon_hemisphere, typ, res)
def build_target_filename(res, typ, lat_hemisphere, abs_lat, lon_hemisphere, abs_lon):
return 'gmted/%02d%s%03d%s_20101117_gmted_%s%s.vrt' % (abs_lat, lat_hemisphere, abs_lon, lon_hemisphere, typ, res)
def add_overviews(filename, ovr_name_list):
f = open(filename, 'rt')
lines = f.readlines()
f.close()
for line in lines:
if line.find('<Overview') != -1:
return
f = open(filename, 'wt')
for line in lines:
if line.find('</VRTRasterBand>') != -1:
for ovr_name in ovr_name_list:
if ovr_name.find('/vsicurl/') == -1:
relative = 1
else:
relative = 0
f.write(' <Overview>\n')
f.write(' <SourceFilename relativeToVRT="%d">%s</SourceFilename>\n' % (relative, ovr_name))
f.write(' <SourceBand>1</SourceBand>\n')
f.write(' </Overview>\n')
f.write('%s' % line)
f.close()
def build_vrt(vrt_name, tile_list):
f = open('vrt_list.txt', 'wt')
for filename in tile_list:
f.write('%s\n' % filename)
f.close()
os.system('gdalbuildvrt -resolution highest -input_file_list vrt_list.txt %s' % vrt_name)
os.unlink('vrt_list.txt')
try:
os.stat('gmted')
except:
os.mkdir('gmted')
gdal.SetConfigOption('CPL_VSIL_CURL_ALLOWED_EXTENSIONS', ".tif")
vrt_list075 = []
vrt_list150 = []
vrt_list300 = []
lon_step = 30
lat_step = 20
lat = 90 - lat_step
while lat >= -90:
abs_lat = abs(lat)
if lat < 0:
lat_hemisphere = 'S'
else:
lat_hemisphere = 'N'
lon = -180
while lon < 180:
abs_lon = abs(lon)
if lon < 0:
lon_hemisphere = 'W'
else:
lon_hemisphere = 'E'
# Special case for southern most tiles : there's only data at resolution 300 and for mea type
if lat == -90:
typ = 'mea'
res = '300'
else:
# We can use the following values for the type :
# bln: Breakline Emphasis
# dsc: Systematic Subsample
# mea: Mean
# med: Median
# min: Minimum
# max: Maximum
# std: Standard Deviation
typ = 'bln'
res = '075'
src_name = build_src_filename(res, typ, lat_hemisphere, abs_lat, lon_hemisphere, abs_lon)
target_name = build_target_filename(res, typ, lat_hemisphere, abs_lat, lon_hemisphere, abs_lon)
try:
os.stat(target_name)
vrt_list075.append(target_name)
except:
print('Generating %s' % target_name)
src_ds = gdal.Open(src_name)
if src_ds is not None:
vrt_ds = gdal.GetDriverByName('VRT').CreateCopy(target_name, src_ds)
vrt_ds = None
vrt_list075.append(target_name)
if target_name in vrt_list075 and res == '075':
ovr_name_list = []
for ovr_res in ['150', '300']:
ovr_name = build_src_filename(ovr_res, typ, lat_hemisphere, abs_lat, lon_hemisphere, abs_lon)
ovr_name_list.append(ovr_name)
add_overviews(target_name, ovr_name_list)
res = '150'
src_name = build_src_filename(res, typ, lat_hemisphere, abs_lat, lon_hemisphere, abs_lon)
target_name = build_target_filename(res, typ, lat_hemisphere, abs_lat, lon_hemisphere, abs_lon)
try:
os.stat(target_name)
vrt_list150.append(target_name)
except:
print('Generating %s' % target_name)
src_ds = gdal.Open(src_name)
if src_ds is not None:
vrt_ds = gdal.GetDriverByName('VRT').CreateCopy(target_name, src_ds)
vrt_ds = None
vrt_list150.append(target_name)
res = '300'
src_name = build_src_filename(res, typ, lat_hemisphere, abs_lat, lon_hemisphere, abs_lon)
target_name = build_target_filename(res, typ, lat_hemisphere, abs_lat, lon_hemisphere, abs_lon)
try:
os.stat(target_name)
vrt_list300.append(target_name)
except:
print('Generating %s' % target_name)
src_ds = gdal.Open(src_name)
if src_ds is not None:
vrt_ds = gdal.GetDriverByName('VRT').CreateCopy(target_name, src_ds)
vrt_ds = None
vrt_list300.append(target_name)
elif target_name in vrt_list075 and res == '300':
vrt_list150.append(target_name)
vrt_list300.append(target_name)
lon += lon_step
lat -= lat_step
build_vrt('gmted/all075.vrt', vrt_list075)
build_vrt('gmted/all150.vrt', vrt_list150)
build_vrt('gmted/all300.vrt', vrt_list300)
add_overviews('gmted/all075.vrt', ['all150.vrt', 'all300.vrt'])
add_overviews('gmted/all150.vrt', ['all300.vrt'])
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment