Geology manipulation workflow
Workflow description:
In [26]:
import os, sys
os.system('unzip ../SGFM-data/2954/srtm_61_17.zip -d ../SGFM-data/2954/DEM')
print "DEM has been unziped in the SGFM Dropbox folder."
2. Test the scene projection with GDAL, make UTM transformation and warp it to simulation domain
These steps are similar to the ones detailed in the previous notebooks (TopoGen_basic and TopoGen).
In [27]:
import os, sys
from osgeo import gdal
from gdalconst import *
import matplotlib.pyplot as plt
# GDAL does not use python exceptions by default
gdal.UseExceptions()
# Open the Geotiff
#gtif = gdal.Open('../SGFM-data/2954/DEM/srtm_61_17.tif',GA_ReadOnly)
# For a full overview of image information use gdalinfo command
#!gdalinfo ../SGFM-data/2954/DEM/srtm_61_17.tif
# Display it
#arr = gtif.ReadAsArray()
#trans = gtif.GetGeoTransform()
#extent = (trans[0], trans[0] + gtif.RasterXSize*trans[1],
# trans[3] + gtif.RasterYSize*trans[5], trans[3])
#figure(figsize = (5,5))
#plt.imshow(arr, extent=extent)
#plt.show()
# Use gdalwarp function to crop the region to Lecode simulation area
!gdalwarp -te 120, -22, 120.5, -21.5 ../SGFM-data/2954/DEM/srtm_61_17.tif \
../SGFM-data/2954/DEM/clip_61_17.tif -overwrite
#!gdalinfo ../SGFM-data/2954/DEM/clip_61_17.tif
# Open the clipped Geotiff
cgtif = gdal.Open('../SGFM-data/2954/DEM/clip_61_17.tif',GA_ReadOnly)
# Create the Numpy array
arr = cgtif.ReadAsArray()
# Find geographic references and extent
trans = cgtif.GetGeoTransform()
extent = (trans[0], trans[0] + cgtif.RasterXSize*trans[1],
trans[3] + cgtif.RasterYSize*trans[5], trans[3])
# Display the array with Matplotlib
figure(figsize = (5,5))
plt.imshow(arr, extent=extent)
plt.show()
Clip the region to the extent of the simulation area
In [28]:
# Convert from WGS84 to UTM
!gdalwarp -t_srs EPSG:32751 ../SGFM-data/2954/DEM/clip_61_17.tif \
../SGFM-data/2954/DEM/utm_61_17.tif -overwrite
# For a full overview of image information use gdalinfo command
#!gdalinfo ../SGFM-data/2954/DEM/utm_61_17.tif
In [29]:
# Crop the region at the given resolution (by defining new resolution in the tr option)
!gdalwarp -te 200000 7570000 240000 7600000 -tr 100 100 -tap -r cubicspline \
../SGFM-data/2954/DEM/utm_61_17.tif ../SGFM-data/2954/DEM/utm_clip.tif -overwrite
# Open the clipped Geotiff
ugtif = gdal.Open('../SGFM-data/2954/DEM/utm_clip.tif',GA_ReadOnly)
# Create the Numpy array
arr = ugtif.ReadAsArray()
# Find geographic references and extent
trans = ugtif.GetGeoTransform()
extent = (trans[0], trans[0] + ugtif.RasterXSize*trans[1],
trans[3] + ugtif.RasterYSize*trans[5], trans[3])
# Display the array with Matplotlib
figure(figsize = (5,5))
plt.imshow(arr, extent=extent)
plt.show()
Analyse DEM for better definition of Lecode XmL input parameters
In [30]:
import os, sys
from osgeo import gdal
from gdalconst import *
import matplotlib.pyplot as plt
from IPython.core.display import Image
# Based on elevation extent use an appropriate color scale to create the DEM map
!gdaldem color-relief ../SGFM-data/2954/DEM/utm_clip.tif ../SGFM-data/color-scale/color_ramp.txt \
../SGFM-data/2954/DEM/region.png -of PNG
# From a DEM raster, output a 32-bit float raster with slope values (specifies in percent slope)
!gdaldem slope ../SGFM-data/2954/DEM/utm_clip.tif ../SGFM-data/2954/DEM/slope_utm.tif -p
# Based on elevation extent use an appropriate color scale to create the DEM slope map
!gdaldem color-relief ../SGFM-data/2954/DEM/slope_utm.tif ../SGFM-data/color-scale/color_slope2.txt \
../SGFM-data/2954/DEM/slopeshade.png -of PNG
# Create an ASCII grid from the simulation area
!gdal_translate -of AAIGrid ../SGFM-data/2954/DEM/utm_clip.tif ../SGFM-data/2954/DEM/dem.asc
# Perform Planchon fill depression algorithm
os.system('fill_depression ../SGFM-data/2954/DEM/dem.asc ../SGFM-data/2954/DEM/fill.asc 0.001')
# Run the multiple flow direction algorithm
os.system('flow_acc ../SGFM-data/2954/DEM/fill.asc ../SGFM-data/2954/DEM/facc.asc')
# Transform the MFD flow accumulation map in Geotiff
!gdal_translate -of "GTiff" ../SGFM-data/2954/DEM/facc.asc ../SGFM-data/2954/DEM/facc.tif
# Based on flow accumulation values use an appropriate color scale to create the DEM slope map
!gdaldem color-relief ../SGFM-data/2954/DEM/facc.tif ../SGFM-data/color-scale/color_facc2.txt \
../SGFM-data/2954/DEM/facc.png -of PNG
# Show the relief map
Image(filename='../SGFM-data/2954/DEM/region.png')
Out[30]:
In [31]:
Image(filename='../SGFM-data/2954/DEM/slopeshade.png')
Out[31]:
In [32]:
Image(filename='../SGFM-data/2954/DEM/facc.png')
Out[32]:
Create Lecode topographic file with an uniform initial sedimentary layer.
This layer will be change according to geological information later but the thickness you define here has to be consistent for the 2 transformation to ensure the correct definition of the topography.
In [42]:
import os, sys
# Create an ASCII grid from the GeoTiff grid
!gdal_translate -of AAIGrid ../SGFM-data/srtm_63_18/clip_utm.tif \
../SGFM-data/srtm_63_18/dem.asc
"""
Perform the transformation to LECODE input format using lecode_basic_input script
- first argument is the ASCII grid you've just generated
- second argument is the thickness of the initial sedimentary layer
- third argument is the number of sediment classes (they will be distribute evenly over the sedimentary layer)
- last argument is the layer hardness
"""
os.system('lecode_basic_input ../SGFM-data/2954/DEM/dem.asc 30 5 1')
print "LECODE input have been created, check for .nodes & .dep files"
In [33]:
import os, zipfile
import numpy as np
from osgeo import ogr
import shapefile, fiona
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from matplotlib.patches import Polygon
from collections import defaultdict
# Open the shapefile
g = ogr.Open( "../SGFM-data/2954/arcview/geology_simple_2954.shp" )
# Get the shapefile layers
layer = g.GetLayer( 0 )
layerDefinition = layer.GetLayerDefn()
# Read the layer attributes from the database file
print "Shapefile attribute"
for i in range(layerDefinition.GetFieldCount()):
print "- %d: %s" % ( i+1, layerDefinition.GetFieldDefn(i).GetName() )
Get the field descriptors from the database which correspond to the Geological records (usually 'CODE').
In [34]:
import fiona
# Using Basemap toolkit library define a 2D plot of the Geological data
fig, ax = subplots(figsize=(7,7))
# Extent of the region
west, east, south, north = 120, 120.5, -22, -21.5
# Define the Basemap projection
m = Basemap(projection='merc', llcrnrlat=south, urcrnrlat=north,
llcrnrlon=west, urcrnrlon=east, lat_ts=0, resolution='l')
colormap = defaultdict(lambda: np.random.random(3))
# Fiona library is then used to find the different rock types within the shapefile
# 'CODE' attribute.
with fiona.open('../SGFM-data/2954/arcview/geology_simple_2954.shp') as f:
for idx,rec in enumerate(f):
coords = rec['geometry']['coordinates'][0]
rocktype = rec['properties']['CODE']
x, y = m(*zip(*coords))
poly = Polygon(zip(x,y), facecolor=colormap[rocktype])
ax.add_patch(poly)
m.drawmapboundary()
m.drawcoastlines()
Out[34]:
In [35]:
g = ogr.Open( "../SGFM-data/2954/arcview/geology_simple_2954.shp" )
layer = g.GetLayer( 0 )
# Find the name of each sedimentary rock present in the database file
code_name = []
for feat in layer :
code_name.append(feat.GetFieldAsString ( "CODE"))
# Sort the list of name to remove duplicates
sorted_code_name = sorted(set(code_name))
print "Total number of geological units (CODE): %d" % (len(sorted_code_name))
print sorted_code_name
4. Define the sediment classes for Lecode simulation
Base on the shapefile geological information we will group the different units based on their names to define different lithologies that will be used in Lecode.
We then extract the information for the defined group in arrays for further processing.
In [36]:
# Grouping rock units by lithological classes
# Number of units after combination
rock_nb = 3
# Define empty list of combined rock units
rock1_name = []
rock2_name = []
rock3_name = []
# Append the list by grouping rock properties based on their names
for name in sorted_code_name :
if 'AD' in name or 'AF' in name or 'AG' in name :
rock1_name.append(name)
elif 'AW' in name or 'Ab' in name or 'Ag' in name :
rock2_name.append(name)
else:
rock3_name.append(name)
# Merge the combined rock units lists into a single rock list
rocks_list = []
rocks_list.append( rock1_name )
rocks_list.append( rock2_name )
rocks_list.append( rock3_name )
Loop through the shapefile rock property attributes and extract the information for the combined lithological classes and perform the projection transformation.
Each combined sedimentary unit is exported to an array defining the presence or absence of the specified lithological classes in the area.
In [39]:
from osgeo import ogr,osr
import gdal
# WGS 84 EPSG reference:
EPSG_shape = 4326
reference_filename = "../SGFM-data/2954/DEM/utm_clip.tif"
target_vector_file = "../SGFM-data/2954/arcview/geology_simple_2954.shp"
rocks_data = []
i = 0
for matnb in range(0,len(rocks_list)):
id = 0
for name in rocks_list[i]:
id += 1
#attribute_filter = "CODE = 'rocks_list[i]'"
attribute_filter = "CODE = '%s'" %( name )
burn_value = 1+i
# First, open the file that we'll be taking as a reference
# We will need to gleam the size in pixels, as well as projection
# and geotransform.
g = gdal.Open( reference_filename )
# We now create an in-memory raster, with the appropriate dimensions
drv = gdal.GetDriverByName('MEM')
target_ds = drv.Create('', g.RasterXSize, g.RasterYSize, 1, gdal.GDT_Byte)
target_ds.SetGeoTransform( g.GetGeoTransform() )
# This goes from WGS84 to the projection in the reference datasets
wgs84 = osr.SpatialReference( ) # Define a SpatialReference object
wgs84.ImportFromEPSG( EPSG_shape ) # And set it to WGS84 using the EPSG code
# Now for the target projection,
to_proj = osr.SpatialReference() # define the SpatialReference object
to_proj.ImportFromWkt( g.GetProjectionRef() )
target_ds.SetProjection ( to_proj.ExportToWkt() )
# Now, we define a coordinate transformtion object, *from* wgs84 *to* UTM
tx = osr.CoordinateTransformation( wgs84, to_proj )
# We define an output in-memory OGR dataset
# You could also do select a driver for an eg "ESRI Shapefile" here
# and give it a sexier name than out!
drv = ogr.GetDriverByName( 'Memory' )
dst_ds = drv.CreateDataSource( 'out' )
# This is a single layer dataset. The layer needs to be of polygons
# and needs to have the target files' projection
dst_layer = dst_ds.CreateLayer('', srs = to_proj, geom_type=ogr.wkbPolygon )
# Open the original shapefile, get the first layer, and filter by attribute
vector_ds = ogr.Open( target_vector_file )
layer = vector_ds.GetLayer ( 0 )
layer.SetAttributeFilter( attribute_filter )
# Get a field definition from the original vector file.
# We don't need much more detail here
feature = layer.GetFeature(0)
field = feature.GetFieldDefnRef( 0 )
# Apply the field definition from the original to the output
dst_layer.CreateField( field )
feature_defn = dst_layer.GetLayerDefn()
# Reset the original layer so we can read all features
layer.ResetReading()
for feat in layer:
# For each feature, get the geometry
geom = feat.GetGeometryRef()
# transform it to the reference projection
geom.Transform ( tx )
# Create an output feature
out_geom = ogr.Feature ( feature_defn )
# Set the geometry to be the reprojected/transformed geometry
out_geom.SetGeometry ( geom )
# Add the feature with its geometry to the output yaer
dst_layer.CreateFeature(out_geom )
# Clear things up
out_geom.Destroy
geom.Destroy
# Done adding geometries
# Reset the output layer to the 0th geometry
dst_layer.ResetReading()
# Now, we rastertize the output vector in-memory file
# into the in-memory output raster file
err = gdal.RasterizeLayer(target_ds, [1], dst_layer,
burn_values=[burn_value])
if err != 0:
print("error:", err)
# Read the data from the raster, this is your mask
data = target_ds.ReadAsArray()
if id == 1 :
ndata = data
if id > 1 :
ndata += data
rocks_data.append( ndata )
i += 1
# Define subplots for each sedimentary types defined previously
fig = plt.figure(figsize=(52, 38))
plt.subplot(1, 3, 1)
plt.imshow ( rocks_data[0], interpolation='nearest', cmap=plt.cm.get_cmap('rainbow', 50), vmin=0, vmax=i )
plt.title('Material Unit 1')
plt.subplot(1, 3, 2)
plt.imshow ( rocks_data[1], interpolation='nearest', cmap=plt.cm.get_cmap('rainbow', 50), vmin=0, vmax=i )
plt.title('Material Unit 2')
plt.subplot(1, 3, 3)
plt.imshow ( rocks_data[2], interpolation='nearest', cmap=plt.cm.get_cmap('rainbow', 50), vmin=0, vmax=i )
plt.title('Material Unit 3')
plt.show()
# Dump the data for further processing
savetxt("../SGFM-data/2954/DEM/rock1.csv", rocks_data[0], fmt="%d")
savetxt("../SGFM-data/2954/DEM/rock2.csv", rocks_data[1], fmt="%d")
savetxt("../SGFM-data/2954/DEM/rock3.csv", rocks_data[2], fmt="%d")
In [40]:
# For each combined rock units define the lithological classes and the hardness
# In this specific case we defined 3 combined rock units
# Hardness per material units:
hardness_unit1 = 1
hardness_unit2 = 5
hardness_unit3 = 10
# We choose to define 5 types of lithologies for our Lecode simulation:
# csand, msand, ssand, silt & mud
nb_sed = 5
rock_comp = []
'''
Distribute lithologies by comnined units base on rock properties
- each sediment types needs to be given a value and set to 0 if not present
- the sediments are expressed in percent (the sum needs to be equal to 100)
- organised sediment from coarser to finer
Note: look on the website for more information on how to set up the material in the XmL input file).
'''
# UNIT 1
sedcomp_unit1 = [0., 20., 30., 30., 20.]
# UNIT 2
sedcomp_unit2 = [20., 20., 30., 20., 10.]
# UNIT 3
sedcomp_unit3 = [40., 30., 20., 5., 5.]
# Check the lithological units to ensure proper declaration
if len( sedcomp_unit1 ) != nb_sed or array( sedcomp_unit1 ).sum() != 100.0 :
err = "sediments declaration number or sum does not match for unit 1"
print("error:", err, len(sedcomp_unit1), array( sedcomp_unit1 ).sum())
sedcomp_unit1.append(hardness_unit1)
rock_comp.append( sedcomp_unit1)
if len(sedcomp_unit2) != nb_sed or array( sedcomp_unit2 ).sum() != 100.0 :
err = "sediments declaration number or sum does not match for unit 2"
print("error:", err, len(sedcomp_unit2), array( sedcomp_unit2 ).sum())
sedcomp_unit2.append(hardness_unit2)
rock_comp.append( sedcomp_unit2)
if len(sedcomp_unit3) != nb_sed or array( sedcomp_unit3 ).sum() != 100.0 :
err = "sediments declaration number or sum does not match for unit 3"
print("error:", err, len(sedcomp_unit3), array( sedcomp_unit3 ).sum())
sedcomp_unit3.append(hardness_unit3)
rock_comp.append( sedcomp_unit3)
# Dump the data for further processing
savetxt("../SGFM-data/2954/DEM/rock_comp.csv", rock_comp, fmt="%.3f" )
print "Grid extension: numrows = %d, numcols = %d" % ( shape( rocks_data[0] ) )
Use geol_2_dep script to transform the geological information in a readable format for Lecode.
In [41]:
import os, sys
'''
# Perform the transformation to LECODE deposit format using geol_2_dep script:
- argument 1: number of lithologies (here nb_sed=5),
- argument 2: extent of the simulation area (numrows, numcols),
- argument 3: sediment distribution by units,
- argument 4 to 4+nb_sed: files containing array of presence or absence of a specific combined rock unit,
- argument 5+nb_sed: thickness of the sedimentary layer (here 30 metres)
(this has to be similar with the thickness you've defined when producing the topographic file),
- argument 6+nb_sed: name of the output file name.
'''
os.system('geol_2_dep 5 3 300 400 ../SGFM-data/2954/DEM/rock_comp.csv \
../SGFM-data/2954/DEM/rock1.csv ../SGFM-data/2954/DEM/rock2.csv \
../SGFM-data/2954/DEM/rock3.csv 30 ../SGFM-data/2954/DEM/layer.dep')
print "LECODE deposit file have been created, check for .dep files"
In []: