This document describes the HydrOffice BAG library and tools (0.2). For the source code, go here.

Source code for hydroffice.bag.elevation

from __future__ import absolute_import, division, print_function, unicode_literals

import os
import logging

log = logging.getLogger(__name__)

import numpy as np

from osgeo import gdal, osr
from .meta import Meta
from .helper import BAGError, Helper
from . import __version__
from .bag import BAGFile

gdal.UseExceptions()


[docs]class Elevation2Gdal(object): formats = { 'ascii': [b"AAIGrid", "bag.elevation.asc"], 'geotiff': [b"GTiff", "bag.elevation.tif"], 'xyz': [b"XYZ", "bag.elevation.xyz"], } def __init__(self, bag_elevation, bag_meta, fmt="geotiff", out_file=None, epsg=None): """Export the elevation layer in one of the listed formats""" assert isinstance(bag_elevation, np.ndarray) assert isinstance(bag_meta, Meta) self.bag_elv = bag_elevation self.bag_meta = bag_meta # get the IN-MEMORY ogr driver self.mem = gdal.GetDriverByName(b"MEM") if self.mem is None: raise BAGError("%s driver not available.\n" % self.formats[fmt][0]) log.debug("format: %s" % fmt) # set the output file self.out_file = out_file if self.out_file is None: self.out_file = os.path.abspath(self.formats[fmt][1]) log.debug("output: %s" % self.out_file) if os.path.exists(self.out_file): os.remove(self.out_file) log.debug("dtype: %s" % self.bag_elv.dtype) self.rst = self.mem.Create(utf8_path=self.out_file, xsize=self.bag_meta.cols, ysize=self.bag_meta.rows, bands=1, eType=gdal.GDT_Float32) self.rst.SetGeoTransform((self.bag_meta.sw[0], self.bag_meta.res_x, 0, self.bag_meta.ne[1], 0, -self.bag_meta.res_y)) self.bnd = self.rst.GetRasterBand(1) self.bnd.WriteArray(self.bag_elv[::-1]) self.bnd.SetNoDataValue(BAGFile.BAG_NAN) self.srs = osr.SpatialReference() if self.bag_meta.wkt_srs is not None: self.srs.ImportFromWkt(self.bag_meta.wkt_srs) else: log.warning("unable to recover valid spatial reference info") self.rst.SetProjection(self.srs.ExportToWkt()) self.bnd.FlushCache() # get the required ogr driver self.drv = gdal.GetDriverByName(self.formats[fmt][0]) # check if re-projection is required if not epsg: # if not, we just create a copy in the selected format dst_ds = self.drv.CreateCopy(self.out_file, self.rst) dst_ds = None self.rst = None return # we need to change projection: # - we create the output srs dst_srs = osr.SpatialReference() dst_srs.ImportFromEPSG(epsg) dst_wkt = dst_srs.ExportToWkt() # Call AutoCreateWarpedVRT() to fetch default values for target raster dimensions and geotransform tmp_ds = gdal.AutoCreateWarpedVRT(self.rst, None, # src_wkt : left to default value --> will use the one from source dst_wkt, gdal.GRA_NearestNeighbour, 0.125 # error threshold --> use same value as in gdalwarp ) # Create the final warped raster dst_ds = self.drv.CreateCopy(self.out_file, tmp_ds) dst_ds = None self.rst = None