added topography class and plotting options with topography

gempy/plotting/plot.py

@@ -625,7 +625,7 @@ def plot_data(geo_data, direction="y", data_type = 'all', series="all", legend_f
     return plot.plot_data(direction=direction, data_type=data_type, series=series, legend_font_size=legend_font_size, **kwargs)
 
 
-def plot_section(geo_data, block, cell_number, direction="y", **kwargs):
+def plot_section(geo_data, block, cell_number, direction="y", topography=None,**kwargs):
     """
     Plot a section of the block model
 
@@ -644,7 +644,7 @@ def plot_section(geo_data, block, cell_number, direction="y", **kwargs):
         None
     """
     plot = PlotData2D(geo_data)
-    plot.plot_block_section(cell_number, block=block, direction=direction, **kwargs)
+    plot.plot_block_section(cell_number, block=block, direction=direction, topography=None, **kwargs)
     # TODO saving options
 
 

gempy/plotting/visualization.py

@@ -230,7 +230,7 @@ def _slice(self, direction, cell_number=25):
         return _a, _b, _c, extent_val, x, y, Gx, Gy
 
     def plot_block_section(self, cell_number=13, block=None, direction="y", interpolation='none',
-                           plot_data=False, ve=1, **kwargs):
+                           plot_data=False, topography = None, ve=1, **kwargs):
         """
         Plot a section of the block model
 
@@ -271,7 +271,7 @@ def plot_block_section(self, cell_number=13, block=None, direction="y", interpol
             self.plot_data(direction, 'all')
         # TODO: plot_topo option - need fault_block for that
 
-        # apply vertical exageration
+        # apply vertical exaggeration
         if direction == 'x' or direction == 'y':
             aspect = ve
         else:
@@ -287,6 +287,10 @@ def plot_block_section(self, cell_number=13, block=None, direction="y", interpol
                         interpolation=interpolation,
                         aspect=aspect,
                         **kwargs)
+        if topography:
+            # TODO: apply vertical exxageration to topography
+            topoline = topography._slice(direction = direction, extent = extent_val, cell_number = cell_number)
+            plt.fill(topoline[:, 0], topoline[:, 1], color='k')
 
         import matplotlib.patches as mpatches
         colors = [im.cmap(im.norm(value)) for value in self.formation_numbers]

gempy/utils/topography.py

@@ -0,0 +1,154 @@
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+import pandas as pn
+import gdal
+import skimage
+
+
+
+class topography():
+    '''Class to combine height elevation data (DEM, file format e.g. tif, asc) with the geological model '''
+
+    def __init__(self, path_dem, geodata=None, output_path=None):
+
+        if path_dem:
+            self.dem = gdal.Open(path_dem)
+        else:
+            print('Define path to file')
+
+        self.dem_zval = self.dem.ReadAsArray()
+        self.raster_extent, self.raster_resolution = self._get_raster_dimensions()
+
+        if geodata:
+            self.geo_data = geodata
+            self.extent_match = self.compare_extent()
+
+        if self.extent_match == False:  # crop dem and update values
+            if output_path:
+                self.dem = self.cropDEM2geodata(output_path)
+                self.dem_zval = self.dem.ReadAsArray()
+                self.raster_extent, self.raster_resolution = self._get_raster_dimensions()
+                self.extent_match = self.compare_extent()
+            else:
+                print('extents of dem and geodata do not match. please define output path for automatic cropping')
+
+        self.grid_info = self._get_grid_info()  # daraus tabelle machen für übersicht
+        self.dem_resized = skimage.transform.resize(self.dem_zval,
+                                                    (self.geo_data.resolution[0], self.geo_data.resolution[1]),
+                                                    preserve_range=True)
+
+        if output_path:  # create file only when extents match
+            self.surface_coordinates = self.convertDEM2xyz(output_path)  # [0]: shape(a,b,3), [1]: shape(a*b,3)
+        else:
+            print('for the full spectrum of plotting with topography, please define an output path')
+
+    def compare_extent(self):
+        if self.geo_data:
+            cornerpoints_geo = self._get_cornerpoints(self.geo_data.extent)
+            cornerpoints_dtm = self._get_cornerpoints(self.raster_extent)
+            if np.any(cornerpoints_geo[:2] - cornerpoints_dtm[:2]) != 0:
+                # print('Extent of geo_data and DEM do not match. Use function cropDEMtogeodata to crop')
+                return False
+            else:
+                # print('geodata and dem extents match')
+                return True
+
+    def show(self, compare=False):
+        print('Raster extent:', self.raster_extent,
+              '\nRaster resolution:', self.raster_resolution)
+        plt.imshow(self.dem_zval, extent=(self.raster_extent[:4]))
+        if compare == True:
+            if self.geo_data:
+                cornerpoints_geo = self._get_cornerpoints(self.geo_data.extent)
+                cornerpoints_dtm = self._get_cornerpoints(self.raster_extent)
+                if self.extent_match == False:
+                    plt.plot(cornerpoints_geo[:, 0], cornerpoints_geo[:, 1], 'ro', markersize=7,
+                             label='Geo_data extent')
+                    plt.plot(cornerpoints_dtm[:, 0], cornerpoints_dtm[:, 1], 'gX', markersize=11, label='DTM extent')
+                    plt.legend(loc=0, fancybox=True, shadow=True)
+                else:
+                    print('geodata and dem extents match')
+
+    def cropDEM2geodata(self, output_path):
+        '''returns new topography object'''
+        path_dest = output_path + '_cropped_DEM.tif'
+        print('Extents of geo_data and DEM do not match. DEM is cropped and stored as', path_dest)
+        new_bounds = (
+        self.geo_data.extent[0], self.geo_data.extent[2], self.geo_data.extent[1], self.geo_data.extent[3])
+        # destName = "C:\\Users\\elisa\\Documents\\git\\MSc\\GempyTopography\\cropped_DTM.tif"
+        gdal.Warp(path_dest, self.dem, options=gdal.WarpOptions(
+            options=['outputBounds'], outputBounds=new_bounds))
+        # cropped_dem = gdal.Open(path_dest)
+        # return topography(path_dest, self.geo_data, output_path)
+        return gdal.Open(path_dest)
+
+    def convertDEM2xyz(self, output_path):
+        '''returns array with the x,y,z coordinates of the topography.'''
+        path_dest = output_path + '_gempytopo.xyz'
+        shape = self.dem_zval.shape
+        gdal.Translate(path_dest, self.dem, options=gdal.TranslateOptions(options=['format'], format="XYZ"))
+        xyz = pn.read_csv(path_dest, header=None, sep=' ').as_matrix()
+        return xyz, np.dstack([xyz[:, 0].reshape(shape), xyz[:, 1].reshape(shape), xyz[:, 2].reshape(shape)])
+
+    def calculate_geomap(self, interpdata, plot=True):
+        geomap, fault = gp.compute_model_at(self.surface_coordinates[1], interpdata)
+        geomap = geomap[0].reshape(self.dem_zval.shape)  # resolution of topo gives much better map
+        if plot:
+            plt.imshow(geomap, origin="lower", cmap=gp.plotting.colors.cmap, norm=gp.plotting.colors.norm)  # set extent
+            plt.title("Geological map", fontsize=15)
+        return geomap
+
+    def _slice(self, direction, extent, cell_number=25):
+        xyzarray = self.surface_coordinates[1]
+        if direction == "x":
+            y = xyzarray[:, :, 1][:, cell_number]  # for y axis
+            z = xyzarray[:, :, 2][:, cell_number]  # for y axis
+            topoline = np.dstack((y, z)).reshape(-1, 2).astype(int)
+            upleft = np.array([extent[0], extent[3]])
+            upright = np.array([extent[1], extent[3]])
+            topolinebox = np.append(topoline, (upleft, upright), axis=0)
+
+        elif direction == "y":
+            y = xyzarray[:, :, 0][cell_number, :]
+            z = xyzarray[:, :, 2][cell_number, :]
+            topoline = np.dstack((y, z)).reshape(-1, 2).astype(int)
+            upleft = np.array([extent[0], extent[3]])
+            upright = np.array([extent[1], extent[3]])
+            topolinebox = np.append(topoline, (upright, upleft), axis=0)
+
+        elif direction == "z":
+            print('not implemented')
+
+        return topolinebox
+
+    def _get_raster_dimensions(self):
+        '''returns dtm.extent and dtm.resolution'''
+        ulx, xres, xskew, uly, yskew, yres = self.dem.GetGeoTransform()
+        z = self.dem_zval
+        if np.any(np.array([xskew, yskew])) != 0:
+            print('Obacht! DEM is not north-oriented.')
+        lrx = ulx + (self.dem.RasterXSize * xres)
+        lry = uly + (self.dem.RasterYSize * yres)
+        res = np.array([(uly - lry) / (-yres), (lrx - ulx) / xres]).astype(int)
+        return np.array([ulx, lrx, lry, uly, z.min(), z.max()]).astype(int), res
+
+    def _get_cornerpoints(self, extent):
+        upleft = ([extent[0], extent[3]])
+        lowleft = ([extent[0], extent[2]])
+        upright = ([extent[1], extent[3]])
+        lowright = ([extent[1], extent[2]])
+        return np.array([upleft, lowleft, upright, lowright])
+
+    def _get_grid_info(self):
+        ext = self.geo_data.extent
+        xres, yres, zres = self.geo_data.resolution[0], self.geo_data.resolution[1], self.geo_data.resolution[2]
+        # ext = self.raster_extent
+        # xres,yres,zres= self.raster_extent[0],self.raster_extent[1],25
+        dx = (ext[1] - ext[0]) / xres
+        dy = (ext[3] - ext[2]) / yres
+        dz = (ext[5] - ext[4]) / zres
+        # return np.array([ext[0],ext[1],xres,dx,ext[2],ext[3],yres,dy,ext[4],ext[5],zres,dz]).reshape(3,4).astype(int)
+        return np.array([ext[0], ext[1], xres, dx, ext[2], ext[3], yres, dy, ext[4], ext[5], zres, dz]).astype(int)
+