Hillup now supports retrieving Canada Digital Elevation Data

Hillup/data/CDED50k.py

@@ -0,0 +1,152 @@
+'''
+Created on 2012-07-06
+
+@author: Edward Sargisson (edward@trailhunger.com)
+'''
+
+from sys import stderr
+from math import floor, ceil
+from ca_national_topographic_system.map_sheet import MapSheetFactory, MapSheet
+from os import unlink, close, write, mkdir, chmod, makedirs
+from os.path import basename, exists, isdir, join
+from httplib import HTTPConnection
+from urlparse import urlparse
+from ftplib import FTP, error_perm
+from osgeo import gdal, osr
+from ca_national_topographic_system import CDEDCell
+from tempfile import mkstemp
+from zipfile import ZipFile
+from hashlib import md5
+import re
+import os
+
+osr.UseExceptions() # <-- otherwise errors will be silent and useless.
+
+sref = osr.SpatialReference()
+sref.ImportFromProj4('+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs')
+
+def quads(minlon, minlat, maxlon, maxlat):
+    """ Generate a list of southeast (lon, lat) for 0.5 * 0.25-degree quads of CDED 1:50000 data.
+    """
+    mapSheetFactory = MapSheetFactory();
+    # CDED 1:50k longitude is divided into map sheets 0.5 deg wide
+    lon = floor(minlon * 2.0) / 2.0
+    while lon <= maxlon:
+        lat = ceil(maxlat * 4.0) / 4.0
+        while lat >= minlat:
+            # We get the northwest point in the loop so we send the south east
+            # point to get the map sheet
+            print 'getting map sheet for lat={0}, lon={1}'.format(lat - .25, lon + .5)
+            mapSheet = mapSheetFactory.getMapSheetBySouthEastCorner(lat - .25, lon + .5)
+            for i in range(2):
+                if i == 0:
+                    cdedCell = CDEDCell(mapSheet, 'EAST')
+                else:
+                    cdedCell = CDEDCell(mapSheet, 'WEST')
+                yield cdedCell
+            #print lon, lat
+            lat -= 0.25
+        lon += 0.5
+
+def datasource(cdedCell, source_dir):
+    mapSheet = cdedCell.mapSheet
+    mapSeriesDir = 'pub/geobase/official/cded/50k_dem/{0}'.format(mapSheet.series)
+    mapSheetFileName = '{0}.zip'.format(mapSheet)
+    mapSheetRetrCommand = 'RETR {0}'.format(mapSheetFileName)
+    mapSheetFullPath = '{0}/{1}'.format(mapSeriesDir, mapSheetFileName)
+
+    #
+    # Create a local filepath
+    #
+    #s, host, path, p, q, f = urlparse(url)
+
+    dem_dir = md5(mapSheetFullPath).hexdigest()[:3]
+    dem_dir = join(source_dir, dem_dir)
+
+    if (cdedCell.half == 'EAST'):
+        dem_file = '{0}_deme.dem'.format(cdedCell.mapSheet)
+    else:
+        dem_file = '{0}_demw.dem'.format(cdedCell.mapSheet)
+    dem_path = join(dem_dir, dem_file)
+    dem_none = dem_path[:-4]+'.404'
+
+    #
+    # Check if the file exists locally
+    #
+    if exists(dem_path):
+        return gdal.Open(dem_path, gdal.GA_ReadOnly)
+
+    if exists(dem_none):
+        return None
+
+    if not exists(dem_dir):
+        makedirs(dem_dir)
+        chmod(dem_dir, 0777)
+
+    assert isdir(dem_dir)
+
+    #
+    # Grab a fresh remote copy
+    #
+    print >> stderr, 'Retrieving', mapSheetFullPath, 'in DEM.CDED50k.datasource().'
+
+
+    ftp = FTP('ftp2.cits.rncan.gc.ca')   # connect to host, default port
+    ftp.login()               # user anonymous, passwd anonymous@
+    ftp.cwd(mapSeriesDir)
+    handle, zip_path = mkstemp(prefix='cded-', suffix='.zip')
+
+    def write_to_handle(data):
+        write(handle, data)
+
+    try:
+        ftp.retrbinary(mapSheetRetrCommand, write_to_handle)
+
+        close(handle)
+
+        #
+        # Get the DEM out of the zip file
+        #
+        zipfile = ZipFile(zip_path, 'r')
+
+        # The file names include a version so we have to look for it.
+        matcher = re.compile('^.{6}_[0-9]{4}_dem[ew].dem$')
+        extractList = filter(matcher.match, zipfile.namelist())
+
+        #
+        # Write the actual DEM
+        #
+        for extractFile in extractList:
+            # We remove the edition and version identification from the file name
+            # so that the cache lookup above will work.
+            extract_path = join(dem_dir, extractFile[0:6] + extractFile[11:])
+            dem_file = open(extract_path, 'w')
+            dem_file.write(zipfile.read(extractFile))
+            dem_file.close()
+
+        chmod(dem_path, 0666)
+    except error_perm:
+        print >> stderr, 'action=datasourceFailed, mapSheetFullPath={0}, reason=ftpPermanentError'.format(mapSheetFullPath)
+        # we're probably outside the coverage area
+        print >> open(dem_none, 'w'), mapSheetFullPath
+        return None
+    # todo put this back
+    #finally:
+        #unlink(zip_path)
+
+    #
+    # The file better exist locally now
+    #
+    return gdal.Open(dem_path, gdal.GA_ReadOnly)
+
+def datasources(minlon, minlat, maxlon, maxlat, source_dir):
+    print 'minlon={0}, minlat={1}, maxlon={2}, maxlat={3}, source_dir={4}'.format(minlon, minlat, maxlon, maxlat, source_dir)
+    """ Retrieve a list of CDED datasources overlapping the tile coordinate.
+    """
+    cdedCells = quads(minlon, minlat, maxlon, maxlat)
+    #cellList = list(cdedCells)
+    #print cellList
+    sources = [datasource(cdedCell, source_dir) for (cdedCell) in cdedCells]
+    return [ds for ds in sources if ds]
+
+

Hillup/data/__init__.py

@@ -5,7 +5,7 @@
 from itertools import product
 from tempfile import mkstemp
 
-import NED10m, NED100m, NED1km, SRTM1, SRTM3
+import NED10m, NED100m, NED1km, SRTM1, SRTM3, CDED50k
 
 from ModestMaps.Core import Coordinate
 from TileStache.Geography import SphericalMercator
@@ -53,7 +53,7 @@ def name(self):
 class Provider:
     """ TileStache provider for generating tiles of DEM slope and aspect data.
 
-        Source parameter can be "srtm-ned" (default) or "ned-only".
+        Source parameter can be "srtm-ned" (default), "ned-only" or "srtm-cded".
 
         See http://tilestache.org/doc/#custom-providers for information
         on how the Provider object interacts with TileStache.
@@ -79,6 +79,10 @@ def renderArea(self, width, height, srs, xmin, ymin, xmax, ymax, zoom):
 
         elif self.source == 'ned-only':
             providers = choose_providers_ned(zoom)
+
+        elif self.source == 'srtm-cded':
+            print "You are using data from Canada Geobase. \n You need to adhere to the GeoBase Unrestricted Use Licence at http://www.geobase.ca/geobase/en/licence.jsp. Specifically, you need to acknowledge attribution."
+            providers = choose_providers_cded(zoom)
 
         else:
             raise Exception('Unknown source "%s"' % source)
@@ -219,6 +223,14 @@ def choose_providers_srtm(zoom):
 
     return [(bottom, proportion), (top, 1 - proportion)]
 
+def choose_providers_cded(zoom):
+    """ Return a list of data sources and proportions for given zoom level.
+
+        Each data source is a module such as SRTM1 or SRTM3, and the proportions
+        must all add up to one. Return list has either one or two items.
+    """
+    return [(CDED50k, 1)]
+
 def choose_providers_ned(zoom):
     """ Return a list of data sources and proportions for given zoom level.
 

ca_national_topographic_system/__init__.py

@@ -0,0 +1,93 @@
+'''Methods for dealing with Canada's National Topographic System (NTS) - 
+specifically finding the NTS Map Sheet number for a latitude, longitude pair.
+'''
+
+from math import floor, ceil
+
+class MapSheetFactory(object):
+    '''
+    Creates map_sheet objects
+    '''
+    map_area_lookup = [
+                ['a', 'b', 'c', 'd'],
+                ['h', 'g', 'f', 'e'],
+                ['i', 'j', 'k', 'l'],
+                ['p', 'o', 'n', 'm']
+                ]
+    map_sheet_lookup = [
+                ['01', '02', '03', '04'],
+                ['08', '07', '06', '05'],
+                ['09', '10', '11', '12'],
+                ['16', '15', '14', '13']
+                ]
+
+    '''
+    Returns a Map Sheet representation for the map sheet specified by its
+    south east corner. 
+    '''
+    def getMapSheetBySouthEastCorner(self, lat, lon):
+
+        # the left most part of the series is based on the longitude
+        series_west = floor((lon + 48.0) / 8.0 * -1)
+        # the right most part of the series is based on the latitude
+        series_north = abs(ceil((lat - 44.0) / 4.0))
+
+        # A Series is divided into 16 Map Areas, marked sinusoidally from A to M starting from the south-east corner of the Series. Ex 96-M.
+        # series orgin is the south east corner
+        series_origin_lon = series_west * -8 - 48
+        series_origin_lat = (series_north  - 1) * 4 + 44
+
+        series_remainder_lon = lon - series_origin_lon
+        series_remainder_lat = lat - series_origin_lat
+
+        # series is 8 deg wide divided into 4
+        map_area_x = int(floor(series_remainder_lon / 2 * -1))
+        map_area_y = int(floor(series_remainder_lat))
+
+#        print 'lat={0}, lon={1}, series_west={2}, series_north={3}, s_o_lat={4}, s_o_lon={5}, s_r_lat={6}, s_r_lon={7}, map_area_x={8}, map_area_y={9}'.format(lat, lon, series_west, series_north, series_origin_lat, series_origin_lon, series_remainder_lat, series_remainder_lon, map_area_x, map_area_y)
+
+        map_area = self.map_area_lookup[map_area_y][map_area_x]
+
+        map_area_origin_lon = map_area_x * 2
+        map_area_origin_lat = map_area_y
+
+        # A Map Area is divided into 16 Map Sheets, marked sinusoidally from 1 to 15 starting from the south-east corner of the Map Area. Ex 96-M-16.
+        map_area_remainder_lon = series_remainder_lon + map_area_origin_lon 
+        map_area_remainder_lat = series_remainder_lat - map_area_origin_lat
+
+        # a map sheet is 0.5 deg wide and 0.25 deg high
+        map_sheet_x = int((map_area_remainder_lon * -1.0) / 0.5)
+        map_sheet_y = int(map_area_remainder_lat / 0.25)
+        map_sheet = self.map_sheet_lookup[map_sheet_y][map_sheet_x]
+
+        # print lat, lon, series_west, series_north, series_origin_lat, series_origin_lon, series_remainder_lat, series_remainder_lon, map_area_x, map_area_y, map_area, map_area_remainder_lat, map_area_remainder_lon, map_sheet_x, map_sheet_y, map_sheet
+        # print  '{0:02.0f}{1:.0f}{2}{3}'.format(series_west, series_north, map_area, map_sheet)
+        series_str = '{0:02.0f}{1:.0f}'.format(series_west, series_north)
+        return MapSheet(series_str, map_area, map_sheet)
+
+class MapSheet(object):
+    def __init__(self, series, mapArea, mapSheet):
+        self.series = series
+        self.mapArea = mapArea
+        self.mapSheet = mapSheet
+
+    def __eq__(self, other):
+        return self.series == other.series and self.mapArea == other.mapArea and self.mapSheet == other.mapSheet
+
+    def __repr__(self):
+        return '{0}{1}{2}'.format(self.series, self.mapArea, self.mapSheet)
+
+class CDEDCell:
+    ''' Represents a Cell from the Canada Digital Elevation Data. Every Map Sheet
+    in the National Topographic System is divided into two halves: a western half and
+    and eastern half. '''
+
+    def __init__(self, mapSheet, half):
+        self.mapSheet = mapSheet
+        self.half = half
+
+    def __eq__(self, other):
+        return self.mapSheet == other.mapSheet and self.half == other.half
+
+    def __repr__(self):
+        return '{0}-{1}'.format(self.mapSheet, self.half)

hillup-seed.py

@@ -33,8 +33,8 @@
                   help='Directory for generated slope/aspect tiles, default "%(tiledir)s". This directory will be used as the "source_dir" for Hillup.tiles:Provider shaded renderings.' % defaults)
 
 parser.add_option('-s', '--source', dest='source',
-                  help='Data source for elevations. One of "srtm-ned" for SRTM and NED data or "ned-only" for US-only downsample NED, default "%(source)s".' % defaults,
-                  choices=('srtm-ned', 'ned-only'))
+                  help='Data source for elevations. One of "srtm-ned" for SRTM and NED data, "ned-only" for US-only downsample NED, "srtm-cded" for SRTM and CDED (Canada Digital Elevation Data). Note that the CDED source does not currently support merging with NED or SRTM at the borders nor the map sheets near the North Pole which do not follow the standard form. default "%(source)s".' % defaults,
+                  choices=('srtm-ned', 'ned-only', 'srtm-cded'))
 
 parser.add_option('--tmp-directory', dest='tmpdir',
                   help='Optional working directory for temporary files. Consider a ram disk for this.')

tests/ca_national_topographic_system_test.py

@@ -0,0 +1,35 @@
+'''
+Created on 2012-07-08
+
+@author: Edward
+'''
+import unittest
+from ca_national_topographic_system import MapSheetFactory, MapSheet
+
+class testMapSheetFactory(unittest.TestCase):
+
+
+    def testShouldGet092j02WhenNearWhistlerBC(self):
+        mapSheetFactory = MapSheetFactory();
+        mapSheet = mapSheetFactory.getMapSheetBySouthEastCorner(50, -122.5)
+        self.assertEquals(mapSheet, MapSheet('092', 'j', '02'))
+
+    def testShouldGet092j03WhenNearWestsideWhistlerBC(self):
+        mapSheetFactory = MapSheetFactory();
+        mapSheet = mapSheetFactory.getMapSheetBySouthEastCorner(50, -123.0)
+        self.assertEquals(mapSheet, MapSheet('092', 'j', '03'))
+
+    def testShouldGet001n08WhenNearStJohnsNL(self):
+        mapSheetFactory = MapSheetFactory();
+        mapSheet = mapSheetFactory.getMapSheetBySouthEastCorner(47.25, -52)
+        self.assertEquals(mapSheet, MapSheet('001', 'n', '08'))
+
+    def testShouldGet030m11WhenNearTorontoOn(self):
+        mapSheetFactory = MapSheetFactory();
+        mapSheet = mapSheetFactory.getMapSheetBySouthEastCorner(43.5, -79)
+        self.assertEquals(mapSheet, MapSheet('030', 'm', '11'))
+
+
+if __name__ == "__main__":
+    #import sys;sys.argv = ['', 'Test.testName']
+    unittest.main()