Source code for cartopy.io.srtm

# Copyright Cartopy Contributors
#
# This file is part of Cartopy and is released under the LGPL license.
# See COPYING and COPYING.LESSER in the root of the repository for full
# licensing details.

"""
The Shuttle Radar Topography Mission (SRTM) is an international research
effort that obtained digital elevation models on a near-global scale from
56S to 60N, to generate the most complete high-resolution digital topographic
database of Earth prior to the release of the ASTER GDEM in 2009.

   - Wikipedia (August 2012)

"""

import io
import os
import warnings

import numpy as np

from cartopy import config
import cartopy.crs as ccrs
from cartopy.io import fh_getter, Downloader, RasterSource, LocatedImage


class _SRTMSource(RasterSource):
    """
    A source of SRTM data, which implements Cartopy's :ref:`RasterSource
    interface <raster-source-interface>`.

    """
    def __init__(self, resolution, downloader, max_nx, max_ny):
        """
        Parameters
        ----------
        resolution
            The resolution of SRTM to download, in integer arc-seconds.
            Data is available at resolutions of 3 and 1 arc-seconds.
        downloader: :class:`cartopy.io.Downloader` instance or None
            The downloader to use for the SRTM dataset. If None, the
            downloader will be taken using
            :class:`cartopy.io.Downloader.from_config` with ('SRTM',
            'SRTM{resolution}') as the target.
        max_nx
            The maximum number of x tiles to be combined when producing a
            wider composite for this RasterSource.
        max_ny
            The maximum number of y tiles to be combined when producing a
            taller composite for this RasterSource.

        """
        if resolution == 3:
            self._shape = (1201, 1201)
        elif resolution == 1:
            self._shape = (3601, 3601)
        else:
            raise ValueError(
                f'Resolution is an unexpected value ({resolution}).')
        self._resolution = resolution

        #: The CRS of the underlying SRTM data.
        self.crs = ccrs.PlateCarree()

        #: The Cartopy Downloader which can handle SRTM data. Normally, this
        #: will be a :class:`SRTMDownloader` instance.
        self.downloader = downloader

        if self.downloader is None:
            self.downloader = Downloader.from_config(
                ('SRTM', f'SRTM{resolution}'))

        #: A tuple of (max_x_tiles, max_y_tiles).
        self._max_tiles = (max_nx, max_ny)

    def validate_projection(self, projection):
        return projection == self.crs

    def fetch_raster(self, projection, extent, target_resolution):
        """
        Fetch SRTM elevation for the given projection and approximate extent.

        """
        if not self.validate_projection(projection):
            raise ValueError(f'Unsupported projection for the '
                             f'SRTM{self._resolution} source.')

        min_x, max_x, min_y, max_y = extent
        min_x, min_y = np.floor([min_x, min_y])
        nx = int(np.ceil(max_x) - min_x)
        ny = int(np.ceil(max_y) - min_y)
        skip = False
        if nx > self._max_tiles[0]:
            warnings.warn(
                f'Required SRTM{self._resolution} tile count ({nx}) exceeds '
                f'maximum ({self._max_tiles[0]}). Increase max_nx limit.')
            skip = True
        if ny > self._max_tiles[1]:
            warnings.warn(
                f'Required SRTM{self._resolution} tile count ({ny}) exceeds '
                f'maximum ({self._max_tiles[1]}). Increase max_ny limit.')
            skip = True
        if skip:
            return []
        else:
            img, _, extent = self.combined(min_x, min_y, nx, ny)
            return [LocatedImage(np.flipud(img), extent)]

    def srtm_fname(self, lon, lat):
        """
        Return the filename for the given lon/lat SRTM tile (downloading if
        necessary), or None if no such tile exists (i.e. the tile would be
        entirely over water, or out of latitude range).

        """
        if int(lon) != lon or int(lat) != lat:
            raise ValueError('Integer longitude/latitude values required.')

        x = '%s%03d' % ('E' if lon >= 0 else 'W', abs(int(lon)))
        y = '%s%02d' % ('N' if lat >= 0 else 'S', abs(int(lat)))

        srtm_downloader = Downloader.from_config(
            ('SRTM', f'SRTM{self._resolution}'))
        params = {'config': config, 'resolution': self._resolution,
                  'x': x, 'y': y}

        # If the URL doesn't exist then we are over sea/north/south of the
        # limits of the SRTM data and we return None.
        if srtm_downloader.url(params) is None:
            return None
        else:
            return self.downloader.path(params)

    def combined(self, lon_min, lat_min, nx, ny):
        """
        Return an image and its extent for the group of nx by ny tiles
        starting at the given bottom left location.

        """
        bottom_left_ll = (lon_min, lat_min)
        shape = np.array(self._shape)
        img = np.zeros(shape * (ny, nx))

        for i, j in np.ndindex(nx, ny):
            x_img_slice = slice(i * shape[1], (i + 1) * shape[1])
            y_img_slice = slice(j * shape[0], (j + 1) * shape[0])

            try:
                tile_img, _, _ = self.single_tile(bottom_left_ll[0] + i,
                                                  bottom_left_ll[1] + j)
            except ValueError:
                img[y_img_slice, x_img_slice] = 0
            else:
                img[y_img_slice, x_img_slice] = tile_img

        extent = (bottom_left_ll[0], bottom_left_ll[0] + nx,
                  bottom_left_ll[1], bottom_left_ll[1] + ny)

        return img, self.crs, extent

    def single_tile(self, lon, lat):
        fname = self.srtm_fname(lon, lat)
        if fname is None:
            raise ValueError('No srtm tile found for those coordinates.')
        return read_SRTM(fname)


[docs]class SRTM3Source(_SRTMSource): """ A source of SRTM3 data, which implements Cartopy's :ref:`RasterSource interface <raster-source-interface>`. """
[docs] def __init__(self, downloader=None, max_nx=3, max_ny=3): """ Parameters ---------- downloader: :class:`cartopy.io.Downloader` instance or None The downloader to use for the SRTM3 dataset. If None, the downloader will be taken using :class:`cartopy.io.Downloader.from_config` with ('SRTM', 'SRTM3') as the target. max_nx The maximum number of x tiles to be combined when producing a wider composite for this RasterSource. max_ny The maximum number of y tiles to be combined when producing a taller composite for this RasterSource. """ super().__init__(resolution=3, downloader=downloader, max_nx=max_nx, max_ny=max_ny)
[docs]class SRTM1Source(_SRTMSource): """ A source of SRTM1 data, which implements Cartopy's :ref:`RasterSource interface <raster-source-interface>`. """
[docs] def __init__(self, downloader=None, max_nx=3, max_ny=3): """ Parameters ---------- downloader: :class:`cartopy.io.Downloader` instance or None The downloader to use for the SRTM1 dataset. If None, the downloader will be taken using :class:`cartopy.io.Downloader.from_config` with ('SRTM', 'SRTM1') as the target. max_nx The maximum number of x tiles to be combined when producing a wider composite for this RasterSource. max_ny The maximum number of y tiles to be combined when producing a taller composite for this RasterSource. """ super().__init__(resolution=1, downloader=downloader, max_nx=max_nx, max_ny=max_ny)
[docs]def srtm(lon, lat): """ Return (elevation, crs, extent) for the given longitude latitude. Elevation is in meters. """ warnings.warn("This method has been deprecated. " "See the \"What's new\" section for v0.12.", DeprecationWarning, stacklevel=2) return SRTM3Source().single_tile(lon, lat)
[docs]def add_shading(elevation, azimuth, altitude): """Add shading to SRTM elevation data, using azimuth and altitude of the sun. Parameters ---------- elevation SRTM elevation data (in meters) azimuth Azimuth of the Sun (in degrees) altitude Altitude of the Sun (in degrees) Return shaded SRTM relief map. """ azimuth = np.deg2rad(azimuth) altitude = np.deg2rad(altitude) x, y = np.gradient(elevation) slope = np.pi/2. - np.arctan(np.sqrt(x*x + y*y)) # -x here because of pixel orders in the SRTM tile aspect = np.arctan2(-x, y) shaded = np.sin(altitude) * np.sin(slope)\ + np.cos(altitude) * np.cos(slope)\ * np.cos((azimuth - np.pi/2.) - aspect) return shaded
[docs]def fill_gaps(elevation, max_distance=10): """Fill gaps in SRTM elevation data for which the distance from missing pixel to nearest existing one is smaller than `max_distance`. This function requires osgeo/gdal to work. Parameters ---------- elevation SRTM elevation data (in meters) max_distance Maximal distance (in pixels) between a missing point and the nearest valid one. Returns ------- elevation SRTM elevation data with filled gaps. """ warnings.warn("The fill_gaps function has been deprecated. " "See the \"What's new\" section for v0.14.", DeprecationWarning, stacklevel=2) # Lazily import osgeo - it is only an optional dependency for cartopy. from osgeo import gdal from osgeo import gdal_array src_ds = gdal_array.OpenArray(elevation) srcband = src_ds.GetRasterBand(1) dstband = srcband maskband = srcband smoothing_iterations = 0 options = [] gdal.FillNodata(dstband, maskband, max_distance, smoothing_iterations, options, callback=None) elevation = dstband.ReadAsArray() return elevation
[docs]def srtm_composite(lon_min, lat_min, nx, ny): warnings.warn("This method has been deprecated. " "See the \"What's new\" section for v0.12.", DeprecationWarning, stacklevel=2) return SRTM3Source().combined(lon_min, lat_min, nx, ny)
[docs]def read_SRTM(fh): """ Read the array of (y, x) elevation data from the given named file-handle. Parameters ---------- fh A named file-like as passed through to :func:`cartopy.io.fh_getter`. The filename is used to determine the extent of the resulting array. Returns ------- elevation The elevation values from the SRTM file. Data is flipped vertically such that the higher the y-index, the further north the data. Data shape is automatically determined by the size of data read from file, and is either (1201, 1201) for 3 arc-second data or (3601, 3601) for 1 arc-second data. crs: :class:`cartopy.crs.CRS` The coordinate reference system of the extents. extents: 4-tuple (x0, x1, y0, y1) The boundaries of the returned elevation array. """ fh, fname = fh_getter(fh, needs_filename=True) if fname.endswith('.zip'): from zipfile import ZipFile zfh = ZipFile(fh, 'rb') fh = zfh.open(os.path.basename(fname[:-4]), 'r') elev = np.fromfile(fh, dtype=np.dtype('>i2')) if elev.size == 12967201: elev.shape = (3601, 3601) elif elev.size == 1442401: elev.shape = (1201, 1201) else: raise ValueError( f'Shape of SRTM data ({elev.size}) is unexpected.') fname = os.path.basename(fname) y_dir, y, x_dir, x = fname[0], int(fname[1:3]), fname[3], int(fname[4:7]) if y_dir == 'S': y *= -1 if x_dir == 'W': x *= -1 return elev[::-1, ...], ccrs.PlateCarree(), (x, x + 1, y, y + 1)
read_SRTM3 = read_SRTM read_SRTM1 = read_SRTM
[docs]def SRTM3_retrieve(lon, lat): """ Return the path of a .hgt file for the given SRTM location. If no such .hgt file exists (because it is over the ocean) None will be returned. """ warnings.warn("This method has been deprecated. " "See the \"What's new\" section for v0.12.", DeprecationWarning, stacklevel=2) return SRTM3Source().srtm_fname(lon, lat)
[docs]class SRTMDownloader(Downloader): """ Provide a SRTM download mechanism. """ FORMAT_KEYS = ('config', 'resolution', 'x', 'y') _SRTM_BASE_URL = ('https://e4ftl01.cr.usgs.gov/MEASURES/' 'SRTMGL{resolution}.003/2000.02.11/') _SRTM_LOOKUP_CACHE = os.path.join(os.path.dirname(__file__), 'srtm.npz') _SRTM_LOOKUP_MASK = np.load(_SRTM_LOOKUP_CACHE)['mask'] """ The SRTM lookup mask determines whether keys such as 'N43E043' are available to download. """
[docs] def __init__(self, target_path_template, pre_downloaded_path_template='', ): # adds some SRTM defaults to the __init__ of a Downloader # namely, the URL is determined on the fly using the # ``SRTMDownloader._SRTM_LOOKUP_MASK`` array Downloader.__init__(self, None, target_path_template, pre_downloaded_path_template)
def url(self, format_dict): warnings.warn('SRTM requires an account set up and log in to access. ' 'Use of this Downloader is likely to fail with HTTP 401 ' 'errors.') # override the url method, looking up the url from the # ``SRTMDownloader._SRTM_LOOKUP_MASK`` array lat = int(format_dict['y'][1:]) # Change to co-latitude. if format_dict['y'][0] == 'N': colat = 90 - lat else: colat = 90 + lat lon = int(format_dict['x'][1:4]) # Ensure positive. if format_dict['x'][0] == 'W': lon = 360 - lon if SRTMDownloader._SRTM_LOOKUP_MASK[lon, colat]: return (SRTMDownloader._SRTM_BASE_URL + '{y}{x}.SRTMGL{resolution}.hgt.zip').format(**format_dict) else: return None def acquire_resource(self, target_path, format_dict): from zipfile import ZipFile target_dir = os.path.dirname(target_path) if not os.path.isdir(target_dir): os.makedirs(target_dir) url = self.url(format_dict) srtm_online = self._urlopen(url) zfh = ZipFile(io.BytesIO(srtm_online.read()), 'r') zip_member_path = '{y}{x}.hgt'.format(**format_dict) member = zfh.getinfo(zip_member_path) with open(target_path, 'wb') as fh: fh.write(zfh.open(member).read()) srtm_online.close() zfh.close() return target_path @staticmethod def _create_srtm_mask(resolution, filename=None): """ Return a NumPy mask of available lat/lon. This is slow as it must query the SRTM server to identify the continent from which the tile comes. Hence a NumPy file with this content exists in ``SRTMDownloader._SRTM_LOOKUP_CACHE``. The NumPy file was created with:: import cartopy.io.srtm as srtm import numpy as np np.savez_compressed(srtm.SRTMDownloader._SRTM_LOOKUP_CACHE, mask=srtm.SRTMDownloader._create_srtm_mask(3)) """ # lazy imports. In most situations, these are not # dependencies of cartopy. from bs4 import BeautifulSoup if filename is None: from urllib.request import urlopen url = SRTMDownloader._SRTM_BASE_URL.format(resolution=resolution) with urlopen(url) as f: html = f.read() else: with open(filename) as f: html = f.read() mask = np.zeros((360, 181), dtype=bool) soup = BeautifulSoup(html) for link in soup('a'): name = str(link.text).strip() if name[0] in 'NS' and name.endswith('.hgt.zip'): lat = int(name[1:3]) # Change to co-latitude. if name[0] == 'N': colat = 90 - lat else: colat = 90 + lat lon = int(name[4:7]) # Ensure positive. if name[3] == 'W': lon = 360 - lon mask[lon, colat] = True return mask @classmethod def default_downloader(cls): """ Return a typical downloader for this class. In general, this static method is used to create the default configuration in cartopy.config """ default_spec = ('SRTM', 'SRTMGL{resolution}', '{y}{x}.hgt') target_path_template = os.path.join('{config[data_dir]}', *default_spec) pre_path_template = os.path.join('{config[pre_existing_data_dir]}', *default_spec) return cls(target_path_template=target_path_template, pre_downloaded_path_template=pre_path_template)
# add a generic SRTM downloader to the config 'downloaders' section. config['downloaders'].setdefault(('SRTM', 'SRTM3'), SRTMDownloader.default_downloader()) config['downloaders'].setdefault(('SRTM', 'SRTM1'), SRTMDownloader.default_downloader())