diff --git a/.gitignore b/.gitignore index 1ff6a80..49530c9 100644 --- a/.gitignore +++ b/.gitignore @@ -18,6 +18,7 @@ gmted/ etopo1/ ned/ ned_topobathy/ +dmr1/ # tile generation directories tiles/ diff --git a/config.example.yaml b/config.example.yaml index 0066cea..de0bb3c 100644 --- a/config.example.yaml +++ b/config.example.yaml @@ -29,8 +29,11 @@ sources: tries: 100 - type: greatlakes - type: srtm - url: http://e4ftl01.cr.usgs.gov/SRTM/SRTMGL1.003/2000.02.11/ - mask-url: http://e4ftl01.cr.usgs.gov/SRTM/SRTMSWBD.003/2000.02.11/ + url: http://e4ftl01.cr.usgs.gov/MEASURES/SRTMGL1.003/2000.02.11 + mask-url: http://e4ftl01.cr.usgs.gov/MEASURES/SRTMSWBD.003/2000.02.11 + username: 'earthdata username' + password: 'earthdata password' + auth-url: https://urs.earthdata.nasa.gov - type: ned13 ftp_server: rockyftp.cr.usgs.gov base_path: vdelivery/Datasets/Staged/NED/13/IMG @@ -40,5 +43,8 @@ sources: - type: ned_topobathy ftp_server: rockyftp.cr.usgs.gov base_path: vdelivery/Datasets/Staged/NED/19/IMG + - type: dmr1 + uri: http://gis.arso.gov.si/lidar/dmr1 + fishnet_url: https://raw.githubusercontent.com/DavixDevelop/TerraLidar/master/mapzen_data/whitelist.yaml logging: config: logging.example.config diff --git a/config.slovenia.yaml b/config.slovenia.yaml new file mode 100644 index 0000000..fd3c6d1 --- /dev/null +++ b/config.slovenia.yaml @@ -0,0 +1,49 @@ +--- +regions: + # san-francisco-bay_california: + # bbox: + # top: 38.719 + # left: -123.640 + # bottom: 36.791 + # right: -121.025 + # zoom_range: [0, 16] + slovenia: + bbox: + top: 46.88333333 + left: 13.39027778 + bottom: 45.40750000 + right: 16.62694444 + zoom_range: [10, 16] +outputs: + - type: skadi + - type: tiff + - type: terrarium + - type: normal +sources: + - type: etopo1 + url: https://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/bedrock/grid_registered/georeferenced_tiff/ETOPO1_Bed_g_geotiff.zip + - type: srtm + url: http://e4ftl01.cr.usgs.gov/MEASURES/SRTMGL1.003/2000.02.11 + mask-url: http://e4ftl01.cr.usgs.gov/MEASURES/SRTMSWBD.003/2000.02.11 + username: 'earthdata username' + password: 'earthdata password' + auth-url: https://urs.earthdata.nasa.gov + - type: dmr1 + uri: http://gis.arso.gov.si/lidar/dmr1 + fishnet_url: https://raw.githubusercontent.com/DavixDevelop/TerraLidar/master/mapzen_data/whitelist.yaml + bbox: + top: 46.88333333 + left: 13.39027778 + bottom: 45.40750000 + right: 16.62694444 +logging: + config: logging.example.config +cluster: + queue: + type: fake +store: + type: file + base_dir: 'render' +source_store: + type: file + base_dir: 'source' \ No newline at end of file diff --git a/config.slovenia_test.yaml b/config.slovenia_test.yaml new file mode 100644 index 0000000..8206336 --- /dev/null +++ b/config.slovenia_test.yaml @@ -0,0 +1,49 @@ +--- +regions: + # san-francisco-bay_california: + # bbox: + # top: 38.719 + # left: -123.640 + # bottom: 36.791 + # right: -121.025 + # zoom_range: [0, 16] + slovenia-test: + bbox: + top: 46.239346 + left: 15.237007 + bottom: 46.205735 + right: 15.288334 + zoom_range: [10, 16] +outputs: + - type: skadi + - type: tiff + - type: terrarium + - type: normal +sources: + - type: etopo1 + url: https://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/bedrock/grid_registered/georeferenced_tiff/ETOPO1_Bed_g_geotiff.zip + - type: srtm + url: http://e4ftl01.cr.usgs.gov/MEASURES/SRTMGL1.003/2000.02.11 + mask-url: http://e4ftl01.cr.usgs.gov/MEASURES/SRTMSWBD.003/2000.02.11 + username: 'earthdata username' + password: 'earthdata password' + auth-url: https://urs.earthdata.nasa.gov + - type: dmr1 + uri: http://gis.arso.gov.si/lidar/dmr1 + fishnet_url: https://raw.githubusercontent.com/DavixDevelop/TerraLidar/master/mapzen_data/whitelist.yaml + bbox: + top: 46.239346 + left: 15.237007 + bottom: 46.205735 + right: 15.288334 +logging: + config: logging.example.config +cluster: + queue: + type: fake +store: + type: file + base_dir: 'render' +source_store: + type: file + base_dir: 'source' \ No newline at end of file diff --git a/docs/attribution.md b/docs/attribution.md index 39d7c8c..60afeec 100644 --- a/docs/attribution.md +++ b/docs/attribution.md @@ -200,3 +200,18 @@ Attribution statement: Attribution statement: > © Commonwealth of Australia (Geoscience Australia) 2017. + +### Digital Elevation Model (DEM) of Slovenia derived from LiDAR 1 Meter Grid (DMR1) + +[Released](http://www.evode.gov.si/index.php?id=69) by Slovenia's [GURS](https://www.e-prostor.gov.si/), licensed under [Creative Commons 4.0](https://creativecommons.org/licenses/by/4.0/deed.sl) + +> Use Limitations: +> General conditions for the use of geodetic data: https://www.e-prostor.gov.si/dostop-do-podatkov/dostop-do-podatkov/#tab1-1029 +> +> Public access Restrictions: +> No restrictions + +Attribution statement: + +> Geodetska uprava Republike Slovenije, LiDAR, DMR1 2015 +> © 2017 MOP - Geodetska uprava Republike Slovenije - Vse pravice pridržane. diff --git a/docs/data-sources.md b/docs/data-sources.md index bfecb46..2212743 100644 --- a/docs/data-sources.md +++ b/docs/data-sources.md @@ -21,6 +21,7 @@ The underlying data sources are a mix of: - [Kartverket](http://data.kartverket.no/download/content/digital-terrengmodell-10-m-utm-33)'s Digital Terrain Model, 10 meters over Norway - [LINZ](https://data.linz.govt.nz/layer/1768-nz-8m-digital-elevation-model-2012/), 8 meters over New Zealand - [SRTM](https://lta.cr.usgs.gov/SRTM) globally except high latitudes, 30 meters (90 meters nominal quality) in land areas +- [DMR1](http://gis.arso.gov.si/evode/profile.aspx?id=atlas_voda_Lidar%40Arso&initialExtent=499500.5%2C109841.5%2C264.58333) derived from pre-processed Lidar, which contains only the terrain at 1m resolution, 1200m-1400m over Slovenia ### Footprints database @@ -89,12 +90,12 @@ zoom | ocean | land **7** | `ETOPO1` | `SRTM`, `NRCAN` in Canada, with `GMTED` in high latitudes above 60° **8** | `ETOPO1` | `SRTM`, `NRCAN` in Canada, with `GMTED` in high latitudes above 60° **9** | `ETOPO1` | `SRTM`, `NRCAN` in Canada, `EUDEM` in Europe, with `GMTED` in high latitudes above 60° -**10** | `ETOPO1`, `NED Topobathy` in California | `SRTM`, `data.gov.at` in Austria, `NRCAN` in Canada, `SRTM`, `NED/3DEP` 1/3 arcsec, `data.gov.uk` in United Kingdom, `INEGI` in Mexico, `ArcticDEM` in latitudes above 60°, `LINZ` in New Zealand, `Kartverket` in Norway -**11** | `ETOPO1`, `NED Topobathy` in California | `SRTM`, `data.gov.at` in Austria, `NRCAN` in Canada, `SRTM`, `NED/3DEP` 1/3 arcsec and 1/9 arcsec, `data.gov.uk` in United Kingdom, `INEGI` in Mexico, `ArcticDEM` in latitudes above 60°, `LINZ` in New Zealand, `Kartverket` in Norway -**12** | `ETOPO1`, `NED Topobathy` in California | `SRTM`, `data.gov.at` in Austria, `NRCAN` in Canada, `SRTM`, `NED/3DEP` 1/3 arcsec and 1/9 arcsec, `data.gov.uk` in United Kingdom, `INEGI` in Mexico, `ArcticDEM` in latitudes above 60°, `LINZ` in New Zealand, `Kartverket` in Norway -**13** | `ETOPO1`, `NED Topobathy` in California | `SRTM`, `data.gov.at` in Austria, `NRCAN` in Canada, `SRTM`, `NED/3DEP` 1/3 arcsec and 1/9 arcsec, `data.gov.uk` in United Kingdom, `INEGI` in Mexico, `ArcticDEM` in latitudes above 60°, `LINZ` in New Zealand, `Kartverket` in Norway -**14** | `ETOPO1`, `NED Topobathy` in California | `SRTM`, `data.gov.at` in Austria, `NRCAN` in Canada, `SRTM`, `NED/3DEP` 1/3 arcsec and 1/9 arcsec, `data.gov.uk` in United Kingdom, `INEGI` in Mexico, `ArcticDEM` in latitudes above 60°, `LINZ` in New Zealand, `Kartverket` in Norway -**15** | `ETOPO1`, `NED Topobathy` in California | `SRTM`, `data.gov.at` in Austria, `NRCAN` in Canada, `SRTM`, `NED/3DEP` 1/3 arcsec and 1/9 arcsec, `data.gov.uk` in United Kingdom, `INEGI` in Mexico, `ArcticDEM` in latitudes above 60°, `LINZ` in New Zealand, `Kartverket` in Norway +**10** | `ETOPO1`, `NED Topobathy` in California | `SRTM`, `data.gov.at` in Austria, `NRCAN` in Canada, `SRTM`, `NED/3DEP` 1/3 arcsec, `data.gov.uk` in United Kingdom, `INEGI` in Mexico, `ArcticDEM` in latitudes above 60°, `LINZ` in New Zealand, `Kartverket` in Norway, `DMR1` in Slovenia +**11** | `ETOPO1`, `NED Topobathy` in California | `SRTM`, `data.gov.at` in Austria, `NRCAN` in Canada, `SRTM`, `NED/3DEP` 1/3 arcsec and 1/9 arcsec, `data.gov.uk` in United Kingdom, `INEGI` in Mexico, `ArcticDEM` in latitudes above 60°, `LINZ` in New Zealand, `Kartverket` in Norway, `DMR1` in Slovenia +**12** | `ETOPO1`, `NED Topobathy` in California | `SRTM`, `data.gov.at` in Austria, `NRCAN` in Canada, `SRTM`, `NED/3DEP` 1/3 arcsec and 1/9 arcsec, `data.gov.uk` in United Kingdom, `INEGI` in Mexico, `ArcticDEM` in latitudes above 60°, `LINZ` in New Zealand, `Kartverket` in Norway, `DMR1` in Slovenia +**13** | `ETOPO1`, `NED Topobathy` in California | `SRTM`, `data.gov.at` in Austria, `NRCAN` in Canada, `SRTM`, `NED/3DEP` 1/3 arcsec and 1/9 arcsec, `data.gov.uk` in United Kingdom, `INEGI` in Mexico, `ArcticDEM` in latitudes above 60°, `LINZ` in New Zealand, `Kartverket` in Norway, `DMR1` in Slovenia +**14** | `ETOPO1`, `NED Topobathy` in California | `SRTM`, `data.gov.at` in Austria, `NRCAN` in Canada, `SRTM`, `NED/3DEP` 1/3 arcsec and 1/9 arcsec, `data.gov.uk` in United Kingdom, `INEGI` in Mexico, `ArcticDEM` in latitudes above 60°, `LINZ` in New Zealand, `Kartverket` in Norway, `DMR1` in Slovenia +**15** | `ETOPO1`, `NED Topobathy` in California | `SRTM`, `data.gov.at` in Austria, `NRCAN` in Canada, `SRTM`, `NED/3DEP` 1/3 arcsec and 1/9 arcsec, `data.gov.uk` in United Kingdom, `INEGI` in Mexico, `ArcticDEM` in latitudes above 60°, `LINZ` in New Zealand, `Kartverket` in Norway, `DMR1` in Slovenia ## Sources native resolution diff --git a/joerd/download.py b/joerd/download.py index 358c780..7b39762 100644 --- a/joerd/download.py +++ b/joerd/download.py @@ -62,6 +62,15 @@ def get(url, options={}): # to restart from the beginning every time). accept_range = False + #auth username + username = options.get('username') + + #auth pass + password = options.get('password') + + #auth top level url + auth_url = options.get('auth-url') + # we need to download _something_ if the file position is less than the # known size, or the size is unknown. while filesize is None or filepos < filesize: @@ -77,6 +86,17 @@ def get(url, options={}): req = urllib2.Request(url) + # if the user provided a username and password, add Authorization + if username is not None: + redirectHandler = urllib2.HTTPRedirectHandler() + cookieProcessor = urllib2.HTTPCookieProcessor() + passwordManager = urllib2.HTTPPasswordMgrWithDefaultRealm() + passwordManager.add_password(None, auth_url, username, password) + authHandler = urllib2.HTTPBasicAuthHandler(passwordManager) + opener = urllib2.build_opener(redirectHandler,cookieProcessor,authHandler) + urllib2.install_opener(opener) + #req.add_header("Authorization", "Basic %s" % base64string) + # if the server supports accept range, and we have a partial # download then attemp to resume it. if accept_range and filepos > 0: @@ -93,7 +113,9 @@ def get(url, options={}): tmp.seek(0, os.SEEK_SET) tmp.truncate(0) + try: + f = urllib2.urlopen(req, timeout=timeout) # try to get the filesize, if the server reports it. @@ -183,4 +205,13 @@ def options(in_opts={}): tries = in_opts.get('tries', 10) out_opts['tries'] = int(tries) - return out_opts + username = in_opts.get('username', None) + out_opts['username'] = username + + password = in_opts.get('password', None) + out_opts['password'] = password + + auth_url = in_opts.get('auth-url', None) + out_opts['auth-url'] = auth_url + + return out_opts \ No newline at end of file diff --git a/joerd/index.py b/joerd/index.py index c247c3f..7b8f208 100644 --- a/joerd/index.py +++ b/joerd/index.py @@ -1,9 +1,11 @@ +from __future__ import print_function import pyqtree import yaml import logging import thread + # Create an index given a YAML file consisting of a list of strings and a # function to parse it. Extra, fixed arguments for the function can be also # be given. Each object returned from the `parse_fn` should have a member @@ -20,6 +22,7 @@ def create(index_file, bbox, parse_fn, *parse_args): if t: idx.insert(bbox=t.bbox.bounds, item=t) n += 1 + print("\r Loaded... %d objects" % n, end="") logger.info("Created index with %d objects." % n) return idx @@ -27,4 +30,4 @@ def create(index_file, bbox, parse_fn, *parse_args): # Returns a list of all objects intersecting the given bbox in the index. def intersections(idx, bbox): - return idx.intersect(bbox.bounds) + return idx.intersect(bbox.bounds) \ No newline at end of file diff --git a/joerd/source/dmr1.py b/joerd/source/dmr1.py new file mode 100644 index 0000000..b816e24 --- /dev/null +++ b/joerd/source/dmr1.py @@ -0,0 +1,285 @@ +from bs4 import BeautifulSoup +from joerd.util import BoundingBox +import joerd.download as download +import joerd.check as check +import joerd.srs as srs +import joerd.index as index +import joerd.mask as mask +import joerd.tmpdir as tmpdir +from joerd.mkdir_p import mkdir_p +from contextlib2 import closing, ExitStack +from shutil import copyfile, move +import os.path +import os +import io +import requests +import logging +import re +import tempfile +import sys +import zipfile +import traceback +import subprocess +import glob +from osgeo import gdal +from osgeo import ogr +from osgeo import osr +import yaml +import time +import subprocess +import mimetypes +from itertools import groupby + +class DMR1Tile(object): + def __init__(self, parent, tile, name, block, bbox): + self.uri = parent.uri + self.download_options = parent.download_options + self.base_dir = parent.base_dir + self.name = name + self.block = block + self.tile = tile + self.bbox = bbox + self.region_bbox = parent.region_bbox + + def __key(self): + return self.name + + def __eq__(a, b): + return isinstance(b, type(a)) and \ + a.__key() == b.__key() + + def __hash__(self): + return hash(self.__key()) + + def urls(self): + uri_list = ["%(uri)s/b_%(block)s/D96TM/TM1_%(name)s.txt" % dict(uri=self.uri,block=self.block,name=self.name)] + return uri_list + + def verifier(self): + return is_txt + + def options(self): + return self.download_options + + def output_file(self): + file_name = "TM1_" + self.name + ".tif" + return os.path.join(self.base_dir, file_name) + + def _tif_file(self): + # returns the name of the geotiff file within the distributed archive. + return "TM1_%(name)s.tif" % dict(name=self.name) + + def unpack(self, store, txt_file): + + with store.upload_dir() as target: + target_dir = os.path.join(target, self.base_dir) + mkdir_p(target_dir) + with tmpdir.tmpdir() as temp_dir: + tif_file = os.path.join(temp_dir, self._tif_file()) + + xyz_file = os.path.join(temp_dir, "TM1_%(name)s.xyz" % dict(name=self.name)) + + #Flip x any y, as they are fliped in dmr1 + data = [] + with open(txt_file.name, 'r') as dmr_file: + for l in dmr_file: + data.append(l.strip().split(";")) + + data = sorted(data, key = lambda x: (float(x[1]), float(x[0]))) + + xyzlist = [] + for l in data: + xyzlist.append('{x};{y};{z}\n'.format(x=l[0],y=l[1],z=l[2])) + with open(xyz_file, 'w') as xyz: + xyz.writelines(xyzlist) + + #Correct NS (North South) resolution and convert to xyz to tif + ds = gdal.Open(xyz_file) + ds = gdal.Warp(tif_file,ds,format="GTiff",dstSRS="EPSG:3794") + ds = None + + #subprocess.check_output("sort -k2 -n -t';' -k1 {tfile} -o {ofile}".format(tfile=txt_file.name, ofile=xyz_file), cwd=temp_dir, shell=True) + #subprocess.check_output("gdalwarp -t_srs EPSG:3794 -ts -overwrite {xfile} {tfile}".format(xfile=xyz_file,tfile=tif_file), cwd=temp_dir, shell=True) + + #Move the file to the store + output_file = os.path.join(target, self.output_file()) + move(tif_file, output_file) + + + def freeze_dry(self): + return dict(type='dmr1', tile=self.tile) + +#IS_DMR1_FILE = re.compile('^b_([0-9]{2})/D96TM/TM1_([0-9]{2,3})_([0-9]{2,3}).txt') +IS_DMR1_FILE = re.compile(r'^b_([0-9]{2})\|TM1_([0-9]{2,3})_([0-9]{2,3})\|([\d*\.?\d*]+)\|([\d*\.?\d*]+)\|([\d*\.?\d*]+)\|([\d*\.?\d*]+)') + +def is_txt(self): + def func(tmp): + if mimetypes.guess_type(tmp.name)[0] == 'text/plain': + return True + else: + return False + return func + + +def _parse_dmr1_tile(tile, parent): + m = IS_DMR1_FILE.match(tile) + + if not m: + return None + + d96_left = float(m.group(4)) + d96_bottom = float(m.group(5)) + d96_right = float(m.group(6)) + d96_top = float(m.group(7)) + + block = int(m.group(1)) + name = "%(x)s_%(y)s" % dict(x=str(m.group(2)),y=str(m.group(3))) + + #D96/TM (EPSG::3794) to WGS84 (EPSG::4326) + src = osr.SpatialReference() + tgt = osr.SpatialReference() + src.ImportFromEPSG(3794) + tgt.ImportFromEPSG(4326) + + transform = osr.CoordinateTransformation(src, tgt) + coords = transform.TransformPoint(d96_left, d96_bottom) + left,bottom = coords[0:2] + coords = transform.TransformPoint(d96_right, d96_top) + right,top= coords[0:2] + + bbox = BoundingBox(left, bottom, right, top) + + #This is so the memory doesn't overflow, due to the big size of the index + #and is to be removed on a production server + if bbox.intersects(parent.region_bbox): + return DMR1Tile(parent, tile, name, block, bbox) + + return None + +class DMR1(object): + def __init__(self, options={}): + self.base_dir = options.get('base_dir', 'dmr1') + self.uri = options['uri'] + self.fishnet_url = options['fishnet_url'] + box = options['bbox'] + self.region_bbox = BoundingBox(box['left'], box['bottom'], box['right'],box['top']) + self.download_options = options + self.tile_index = None + + def get_index(self): + if not os.path.isdir(self.base_dir): + os.makedirs(self.base_dir) + index_name = 'index.yaml' + index_file = os.path.join(self.base_dir, index_name) + # if index doesn't exist, or is more than 24h old + if not os.path.isfile(index_file) or \ + time.time() > os.path.getmtime(index_file) + 86400: + self.download_index(index_file) + + def download_index(self, index_file): + logger = logging.getLogger('dmr1') + logger.info('Fetcthing D96TM index...') + + """ + fishnet = 'LIDAR_FISHNET_D96.shp' + + #these tiles do not have dmr1 + blacklist = ['b_21/D96TM/TM1_393_35.txt', + 'b_21/D96TM/TM1_392_35.txt', + 'b_23/D96TM/TM1_509_169.txt' + 'b_24/D96TM/TM1_616_148.txt', + 'b_24/D96TM/TM1_617_148.txt', + 'b_24/D96TM/TM1_618_148.txt', + 'b_24/D96TM/TM1_622_147.txt', + 'b_24/D96TM/TM1_622_148.txt', + 'b_24/D96TM/TM1_623_148.txt', + 'b_24/D96TM/TM1_623_149.txt'] + + req = requests.get(self.fishnet_url, stream=True) + with tmpdir.tmpdir() as d: + with zipfile.ZipFile(io.BytesIO(req.content)) as zip_file: + zip_file.extractall(d) + + fishnet_file = os.path.join(d, fishnet) + driver = ogr.GetDriverByName("ESRI Shapefile") + dataSource = driver.Open(fishnet_file, 1) + layer = dataSource.GetLayer() + + for feature in layer: + link = feature.GetField("BLOK") + "/D96TM/TM1_" + feature.GetField("NAME") + ".txt"; + if link not in blacklist: + links.append(link) + layer.ResetReading() + """ + req = requests.get(self.fishnet_url, allow_redirects=True) + + with open(index_file, 'w') as file: + file.write(req.content) + #file.write(yaml.dump(whitelist)) + + def _ensure_tile_index(self): + if self.tile_index is None: + index_file = os.path.join(self.base_dir, 'index.yaml') + + bbox = (13.39027778,45.40750000,16.62694444,46.88333333) + self.tile_index = index.create(index_file, bbox, _parse_dmr1_tile, self) + + return self.tile_index + + def existing_files(self): + for base, dirs, files in os.walk(self.base_dir): + for f in files: + if f.endswith('tif'): + yield os.path.join(base, f) + + def rehydrate(self, data): + assert data.get('type') == 'dmr1', \ + "Unable to rehydrate %r from Slovenia." %data + return _parse_dmr1_tile(data['tile'], self) + + def downloads_for(self, tile): + tiles = set() + # if the tile scale is greater than 20x the D96TM scale, then there's no + # point in including D96TM, it'll be far too fine to make a difference. + # D96TM is 1m (Aboue same as 1/9th arc second). + + if tile.max_resolution() > 20 * 1.0 / (3600 * 9): + return tiles + + # buffer by 0.0075 degrees (81px) to grab neighbouring tiles and ensure + # some overlap to take care of boundary issues. + tile_bbox = tile.latlon_bbox().buffer(0.0075) + + tile_index = self._ensure_tile_index() + + for t in index.intersections(tile_index, tile_bbox): + tiles.add(t) + + return tiles + + def filter_type(self, src_res, dst_res): + return gdal.GRA_Lanczos if src_res > dst_res else gdal.GRA_Cubic + + def vrts_for(self, tile): + """ + Returns a list of sets of tiles, with each list element intended as a + separate VRT for use in GDAL. + + D96TM is overlapping. + """ + vrt = [] + tiles = self.downloads_for(tile) + + def func(tile): + return (tile.tile) + + for k, t in groupby(sorted(tiles, key=func), func): + vrt.append(set(t)) + + return vrt + + def srs(self): + return srs.d96() + +def create(options): + return DMR1(options) \ No newline at end of file diff --git a/joerd/source/srtm.py b/joerd/source/srtm.py index 7af8386..a15859b 100644 --- a/joerd/source/srtm.py +++ b/joerd/source/srtm.py @@ -9,6 +9,7 @@ from joerd.mkdir_p import mkdir_p from contextlib2 import closing, ExitStack from shutil import copyfile, move +from ftplib import FTP import os.path import os import requests @@ -82,12 +83,12 @@ def _unpack_hgt(self, zip_name, target_dir): zfile.extract(n, target_dir) if n != self.fname: move(os.path.join(target_dir, n), - os.path.join(target_dir, self.fname)) + os.path.join(target_dir, self.fname)) return raise LookupError("None of the alternative names %r were found " - "in the SRTM zipfile %r. Contents are: %r" % - (names, zip_name, exists)) + "in the SRTM zipfile %r. Contents are: %r" % + (names, zip_name, exists)) def unpack(self, store, data_zip, mask_zip=None): with store.upload_dir() as target: @@ -112,13 +113,18 @@ def unpack(self, store, data_zip, mask_zip=None): # mask off the water using the mask raster raw file output_file = os.path.join(target, self.output_file()) mask.raw(os.path.join(d, self.fname), mask_file, 255, - "SRTMHGT", output_file) + "SRTMHGT", output_file) def freeze_dry(self): return dict(type='srtm', link=self.link, is_masked=self.is_masked) def _parse_srtm_tile(link, parent, is_masked=None): + i = IS_SRTM_FILE.match(link) + + if not i: + return None + fname = link.replace(".SRTMGL1.hgt.zip", ".hgt") bbox = parent._parse_bbox(link) if is_masked is None: @@ -128,7 +134,6 @@ def _parse_srtm_tile(link, parent, is_masked=None): class SRTM(object): - def __init__(self, options={}): self.base_dir = options.get('base_dir', 'srtm') self.url = options['url'] @@ -155,7 +160,7 @@ def get_one_index(self, name): index_file = os.path.join(self.base_dir, fname) # if index doesn't exist, or is more than 24h old if not os.path.isfile(index_file) or \ - time.time() > os.path.getmtime(index_file) + 86400: + time.time() > os.path.getmtime(index_file) + 86400: self.download_index(index_file, name) def download_index(self, index_file, name): @@ -167,30 +172,40 @@ def download_index(self, index_file, name): url = None if name == 'tile': - url = self.url + url = self.url if name == 'mask': url = self.mask_url - r = requests.get(url) + url_r = url + "/" + + r = requests.get(url.replace('http','https')) soup = BeautifulSoup(r.text, 'html.parser') links = [] for a in soup.find_all('a'): link = a.get('href') if link is not None: - bbox = self._parse_bbox(link) - if bbox: - links.append(link) + if link.find(url.replace('http','https')) != -1: + req = requests.get(link) + isoup = BeautifulSoup(req.text, 'html.parser') + for aa in isoup.find_all('a'): + ilink = aa.get('href') + if ilink is not None: + fname = ilink.replace(url_r, '') + bbox = self._parse_bbox(fname) + if bbox: + links.append(fname) with open(index_file, 'w') as f: f.write(yaml.dump(links)) + + def _ensure_tile_index(self): if self.tile_index is None: index_file = os.path.join(self.base_dir, 'index_tile.yaml') bbox = (-180, -90, 180, 90) - self.tile_index = index.create(index_file, bbox, _parse_srtm_tile, - self) + self.tile_index = index.create(index_file, bbox, _parse_srtm_tile,self) return self.tile_index @@ -271,4 +286,4 @@ def _parse_bbox(self, link): def create(options): - return SRTM(options) + return SRTM(options) \ No newline at end of file diff --git a/joerd/srs.py b/joerd/srs.py index 45923ad..d3a3476 100644 --- a/joerd/srs.py +++ b/joerd/srs.py @@ -12,6 +12,15 @@ 'AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0],UNIT["degree",' \ '0.0174532925199433],AUTHORITY["EPSG","4269"]]' +D96_WKT = 'PROJCS["Slovenia 1996 / Slovene National Grid",GEOGCS["Slovenia 1996",' \ +'DATUM["Slovenia_Geodetic_Datum_1996",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],' \ +'TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6765"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],' \ +'UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4765"]],PROJECTION["Transverse_Mercator"],' \ +'PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",15],PARAMETER["scale_factor",0.9999],' \ +'PARAMETER["false_easting",500000],PARAMETER["false_northing",-5000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],' \ +'AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","3794"]]' + + def wgs84(): sr = osr.SpatialReference() @@ -23,3 +32,8 @@ def nad83(): sr = osr.SpatialReference() sr.ImportFromWkt(NAD83_WKT) return sr + +def d96(): + sr = osr.SpatialReference() + sr.ImportFromWkt(D96_WKT) + return sr \ No newline at end of file