GAB geology manipulation workflow
Workflow description:
Work with the topography dataset
In [3]:
import os, sys
from osgeo import gdal
from gdalconst import *
import matplotlib.pyplot as plt
from IPython.core.display import Image
# GDAL does not use python exceptions by default
gdal.UseExceptions()
# Combine several tiles together topography into a unique Geotiff
!gdal_merge.py /Users/sgfm/GAB/srtm_62_19/srtm_62_19.tif /Users/sgfm/GAB/srtm_63_19/srtm_63_19.tif \
-o /Users/sgfm/GAB/data/merge_topo.tif
# For a full overview of image information use gdalinfo command
#!gdalinfo -mm /Users/sgfm/GAB/data/merge_topo.tif
# Use gdalwarp function to crop the region to Lecode simulation area
!gdalwarp -te 126, -35, 135, -30 /Users/sgfm/GAB/data/merge_topo.tif \
/Users/sgfm/GAB/data/merge_clip.tif -overwrite
# Define the new projection system and apply a new resolution (by defining the tr option)
!gdalwarp -t_srs EPSG:32752 -tr 1000 1000 -r bilinear -te 240000, 6200000, 1010000, 6600000 -ot Float32 \
/Users/sgfm/GAB/data/merge_clip.tif /Users/sgfm/GAB/data/topo_clip.tif -overwrite
#!gdalinfo -mm /Users/sgfm/GAB/data/topo_utm.tif
# Using gdal_calc.py function, define topographic values below 0 as nodata values
!gdal_calc.py -A /Users/sgfm/GAB/data/topo_utm.tif --outfile=/Users/sgfm/GAB/data/topo_calc_nodata.tif \
--calc="A*(A>0)" --NoDataValue=0 --overwrite
# Then transform nodata values as none
!gdal_translate -a_nodata none /Users/sgfm/GAB/data/topo_calc_nodata.tif \
/Users/sgfm/GAB/data/topo_calc.tif
# Open topography Geotiff
gtif = gdal.Open('/Users/sgfm/GAB/data/topo_calc.tif',GA_ReadOnly)
#!gdalinfo -mm /Users/sgfm/GAB/data/topo_calc.tif
# Create the Numpy array
arr = gtif.ReadAsArray()
# Find geographic references and extent
trans = gtif.GetGeoTransform()
extent = (trans[0], trans[0] + gtif.RasterXSize*trans[1],
trans[3] + gtif.RasterYSize*trans[5], trans[3])
# Display the array with Matplotlib
figure(figsize = (7,7))
plt.imshow(arr, extent=extent)
plt.show()
Work with the bathymetry dataset
In [4]:
# The bathemetry is obtained from ArcInfo binary coverages datasets
# First translate the ArcInfo binary coverages to Geotiff
!gdal_translate /Users/sgfm/GAB/71989/GA_bathymetry/bathy_01_flt /Users/sgfm/GAB/data/bathymetry.tif
#!gdalinfo -mm /Users/sgfm/GAB/data/bathymetry.tif
# Use gdalwarp function to crop the region to Lecode simulation area
!gdalwarp -te 126, -35, 135, -30 /Users/sgfm/GAB/data/bathymetry.tif \
/Users/sgfm/GAB/data/bathymetry_clip.tif -overwrite -dstnodata -9999999
#!gdalinfo -mm /Users/sgfm/GAB/data/bathymetry_clip.tif
# Define the new projection system and apply a new resolution (by defining the tr option)
!gdalwarp -t_srs EPSG:32752 -tr 1000 1000 -r bilinear -te 240000, 6200000, 1010000, 6600000 \
/Users/sgfm/GAB/data/bathymetry_clip.tif /Users/sgfm/GAB/data/bathymetry_clip2.tif -overwrite
# Using gdal_calc.py function, define bathymetric values above 0 as nodata values
!gdal_calc.py -A /Users/sgfm/GAB/data/bathy_utm.tif --outfile=/Users/sgfm/GAB/data/bathy_calc_nodata.tif \
--calc="A*(A<0.0)" --NoDataValue=0 --overwrite
# Then transform nodata values as none
!gdal_translate -a_nodata none /Users/sgfm/GAB/data/bathy_calc_nodata.tif \
/Users/sgfm/GAB/data/bathy_calc.tif
# Open bathymetry Geotiff
btif = gdal.Open('/Users/sgfm/GAB/data/bathy_calc.tif',GA_ReadOnly)
#!gdalinfo -mm /Users/sgfm/GAB/data/bathy_calc.tif
# Create the Numpy array
arr = btif.ReadAsArray()
# Find geographic references and extent
trans = btif.GetGeoTransform()
extent = (trans[0], trans[0] + btif.RasterXSize*trans[1],
trans[3] + btif.RasterYSize*trans[5], trans[3])
# Display the array with Matplotlib
figure(figsize = (7,7))
plt.imshow(arr, extent=extent)
plt.show()
Combined topography/bathymetry dataset
In [19]:
# Combine topography and bathymetry using gdal_calc.py function
!gdal_calc.py -A /Users/sgfm/GAB/data/bathy_calc.tif -B /Users/sgfm/GAB/data/topo_calc.tif \
--outfile=/Users/sgfm/GAB/data/gab_topo_nd.tif --calc="B+A" --overwrite
# Open bathymetry Geotiff
ttif = gdal.Open('/Users/sgfm/GAB/data/gab_topo_nd.tif',GA_ReadOnly)
# Create the Numpy array
arr = ttif.ReadAsArray()
# Find geographic references and extent
trans = ttif.GetGeoTransform()
extent = (trans[0], trans[0] + ttif.RasterXSize*trans[1],
trans[3] + ttif.RasterYSize*trans[5], trans[3])
# Display the array with Matplotlib
figure(figsize = (7,7))
plt.imshow(arr, extent=extent)
plt.show()
# Set the nodata values to none
!gdal_translate -a_nodata none /Users/sgfm/GAB/data/gab_topo_nd.tif \
/Users/sgfm/GAB/data/gab_topo.tif
Visualise DEM information
In [20]:
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 /Users/sgfm/GAB/data/gab_topo.tif ../SGFM-data/color-scale/color_gab.txt \
/Users/sgfm/GAB/data/region.png -of PNG
# From a DEM raster, output a 32-bit float raster with slope values (specifies in percent slope)
!gdaldem slope /Users/sgfm/GAB/data/gab_topo.tif /Users/sgfm/GAB/data/slope_utm.tif -p
# Based on elevation extent use an appropriate color scale to create the DEM slope map
!gdaldem color-relief /Users/sgfm/GAB/data/slope_utm.tif ../SGFM-data/color-scale/color_slope2.txt \
/Users/sgfm/GAB/data/slopeshade.png -of PNG
# Create an ASCII grid from the simulation area
!gdal_translate -of AAIGrid /Users/sgfm/GAB/data/gab_topo.tif /Users/sgfm/GAB/data/dem.asc
# Perform Planchon fill depression algorithm
#os.system('fill_depression /Users/sgfm/GAB/data/dem.asc /Users/sgfm/GAB/data/fill.asc 0.001')
# Run the multiple flow direction algorithm
os.system('flow_acc /Users/sgfm/GAB/data/dem.asc /Users/sgfm/GAB/data/facc.asc')
# Transform the MFD flow accumulation map in Geotiff
!gdal_translate -of "GTiff" /Users/sgfm/GAB/data/facc.asc /Users/sgfm/GAB/data/facc.tif
# Based on flow accumulation values use an appropriate color scale to create the DEM slope map
!gdaldem color-relief /Users/sgfm/GAB/data/facc.tif ../SGFM-data/color-scale/color_facc2.txt \
/Users/sgfm/GAB/data/facc.png -of PNG
# Show the relief map
Image(filename='/Users/sgfm/GAB/data/region.png')
Out[20]:
In [21]:
Image(filename='/Users/sgfm/GAB/data/slopeshade.png')
Out[21]:
In [22]:
Image(filename='/Users/sgfm/GAB/data/facc.png')
Out[22]:
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 [34]:
# Create Lecode node file and uniform deposit file
# that will be overwritten with the proper geology
import os, sys
# Create an ASCII grid from the GeoTiff grid
!gdal_translate -of AAIGrid /Users/sgfm/GAB/data/gab_topo.tif \
../SGFM-data/GAB/gab_topo.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/GAB/gab_topo.asc 50 6 1')
print "LECODE input have been created, check for .nodes & .dep files"
Export the carbonate map
In [28]:
# The carbonate geology is obtained from ArcInfo binary coverages datasets
# First translate the ArcInfo binary coverages to Geotiff
!gdal_translate /Users/sgfm/GAB/71991/GA_sediments/bulk_carb1 /Users/sgfm/GAB/data/sed_carb.tif
# Use gdalwarp function to crop the region to Lecode simulation area
!gdalwarp -te 126, -35, 135, -30 /Users/sgfm/GAB/data/sed_carb.tif \
/Users/sgfm/GAB/data/sed_carb_clip.tif -overwrite -dstnodata -9999999
#!gdalinfo -mm /Users/sgfm/GAB/data/bathymetry_clip.tif
# Define the new projection system and apply a new resolution (by defining the tr option)
!gdalwarp -t_srs EPSG:32752 -tr 1000 1000 -r bilinear -te 240000, 6200000, 1010000, 6600000 \
/Users/sgfm/GAB/data/sed_carb_clip.tif /Users/sgfm/GAB/data/sed_carb_clip2.tif -overwrite
# Using gdal_calc.py function, define carbonate values above 0 as nodata values
!gdal_calc.py -A /Users/sgfm/GAB/data/sed_carb_utm.tif --outfile=/Users/sgfm/GAB/data/sed_carb_nodata.tif \
--calc="A*(A>0.0)" --NoDataValue=0 --overwrite
# Then transform nodata values as none
!gdal_translate -a_nodata none /Users/sgfm/GAB/data/sed_carb_nodata.tif \
/Users/sgfm/GAB/data/carbonate.tif
#!gdalinfo -mm /Users/sgfm/GAB/data/carbonate.tif
# Open carbonate Geotiff
carbtif = gdal.Open('/Users/sgfm/GAB/data/carbonate.tif',GA_ReadOnly)
# Create the Numpy array
arr = carbtif.ReadAsArray()
# Find geographic references and extent
trans = carbtif.GetGeoTransform()
extent = (trans[0], trans[0] + carbtif.RasterXSize*trans[1],
trans[3] + carbtif.RasterYSize*trans[5], trans[3])
# Display the array with Matplotlib
figure(figsize = (7,7))
plt.imshow(arr, extent=extent)
plt.show()
Export the gravel map
In [30]:
# The gravel geology is obtained from ArcInfo binary coverages datasets
# First translate the ArcInfo binary coverages to Geotiff
!gdal_translate /Users/sgfm/GAB/71991/GA_sediments/gravel /Users/sgfm/GAB/data/sed_grav.tif
#!gdalinfo -mm /Users/sgfm/GAB/data/bathymetry.tif
# Use gdalwarp function to crop the region to Lecode simulation area
!gdalwarp -te 126, -35, 135, -30 /Users/sgfm/GAB/data/sed_grav.tif \
/Users/sgfm/GAB/data/sed_grav_clip.tif -overwrite -dstnodata -9999999
#!gdalinfo -mm /Users/sgfm/GAB/data/bathymetry_clip.tif
# Define the new projection system and apply a new resolution (by defining the tr option)
!gdalwarp -t_srs EPSG:32752 -tr 1000 1000 -r bilinear -te 240000, 6200000, 1010000, 6600000 \
/Users/sgfm/GAB/data/sed_grav_clip.tif /Users/sgfm/GAB/data/sed_grav_clip2.tif -overwrite
# Using gdal_calc.py function, define carbonate values above 0 as nodata values
!gdal_calc.py -A /Users/sgfm/GAB/data/sed_grav_utm.tif --outfile=/Users/sgfm/GAB/data/sed_grav_nodata.tif \
--calc="A*(A>0.0)" --NoDataValue=0 --overwrite
# Then transform nodata values as none
!gdal_translate -a_nodata none /Users/sgfm/GAB/data/sed_grav_nodata.tif \
/Users/sgfm/GAB/data/gravel.tif
#!gdalinfo -mm /Users/sgfm/GAB/data/gravel.tif
# Open gravel Geotiff
gravtif = gdal.Open('/Users/sgfm/GAB/data/gravel.tif',GA_ReadOnly)
# Create the Numpy array
arr = gravtif.ReadAsArray()
# Find geographic references and extent
trans = gravtif.GetGeoTransform()
extent = (trans[0], trans[0] + gravtif.RasterXSize*trans[1],
trans[3] + gravtif.RasterYSize*trans[5], trans[3])
# Display the array with Matplotlib
figure(figsize = (7,7))
plt.imshow(arr, extent=extent)
plt.show()
Export the sand map
In [31]:
# The sand geology is obtained from ArcInfo binary coverages datasets
# First translate the ArcInfo binary coverages to Geotiff
!gdal_translate /Users/sgfm/GAB/71991/GA_sediments/sand2 /Users/sgfm/GAB/data/sed_sand.tif
#!gdalinfo -mm /Users/sgfm/GAB/data/bathymetry.tif
# Use gdalwarp function to crop the region to Lecode simulation area
!gdalwarp -te 126, -35, 135, -30 /Users/sgfm/GAB/data/sed_sand.tif \
/Users/sgfm/GAB/data/sed_sand_clip.tif -overwrite -dstnodata -9999999
#!gdalinfo -mm /Users/sgfm/GAB/data/bathymetry_clip.tif
# Use gdalwarp function to crop the region to Lecode simulation area
!gdalwarp -t_srs EPSG:32752 -tr 1000 1000 -r bilinear -te 240000, 6200000, 1010000, 6600000 \
/Users/sgfm/GAB/data/sed_sand_clip.tif /Users/sgfm/GAB/data/sed_sand_clip2.tif -overwrite
# Using gdal_calc.py function, define carbonate values above 0 as nodata values
!gdal_calc.py -A /Users/sgfm/GAB/data/sed_sand_utm.tif --outfile=/Users/sgfm/GAB/data/sed_sand_nodata.tif \
--calc="A*(A>0.0)" --NoDataValue=0 --overwrite
# Then transform nodata values as none
!gdal_translate -a_nodata none /Users/sgfm/GAB/data/sed_sand_nodata.tif \
/Users/sgfm/GAB/data/sand.tif
#!gdalinfo -mm /Users/sgfm/GAB/data/sand.tif
# Open sand Geotiff
sandtif = gdal.Open('/Users/sgfm/GAB/data/sand.tif',GA_ReadOnly)
# Create the Numpy array
arr = sandtif.ReadAsArray()
# Find geographic references and extent
trans = sandtif.GetGeoTransform()
extent = (trans[0], trans[0] + sandtif.RasterXSize*trans[1],
trans[3] + sandtif.RasterYSize*trans[5], trans[3])
# Display the array with Matplotlib
figure(figsize = (7,7))
plt.imshow(arr, extent=extent)
plt.show()
Export the mud map
In [32]:
# The mud geology is obtained from ArcInfo binary coverages datasets
# First translate the ArcInfo binary coverages to Geotiff
!gdal_translate /Users/sgfm/GAB/71991/GA_sediments/mud /Users/sgfm/GAB/data/sed_mud.tif
#!gdalinfo -mm /Users/sgfm/GAB/data/bathymetry.tif
# Use gdalwarp function to crop the region to Lecode simulation area
!gdalwarp -te 126, -35, 135, -30 /Users/sgfm/GAB/data/sed_mud.tif \
/Users/sgfm/GAB/data/sed_mud_clip.tif -overwrite -dstnodata -9999999
#!gdalinfo -mm /Users/sgfm/GAB/data/bathymetry_clip.tif
!gdalwarp -t_srs EPSG:32752 -tr 1000 1000 -te 240000, 6200000, 1010000, 6600000 \
-r bilinear /Users/sgfm/GAB/data/sed_mud_clip.tif /Users/sgfm/GAB/data/sed_mud_clip2.tif -overwrite
# Using gdal_calc.py function, define carbonate values above 0 as nodata values
!gdal_calc.py -A /Users/sgfm/GAB/data/sed_mud_utm.tif --outfile=/Users/sgfm/GAB/data/sed_mud_nodata.tif \
--calc="A*(A>0.0)" --NoDataValue=0 --overwrite
# Then transform nodata values as none
!gdal_translate -a_nodata none /Users/sgfm/GAB/data/sed_mud_nodata.tif \
/Users/sgfm/GAB/data/mud.tif
#!gdalinfo -mm /Users/sgfm/GAB/data/mud.tif
# Open mud Geotiff
mudtif = gdal.Open('/Users/sgfm/GAB/data/mud.tif',GA_ReadOnly)
# Create the Numpy array
arr = mudtif.ReadAsArray()
# Find geographic references and extent
trans = mudtif.GetGeoTransform()
extent = (trans[0], trans[0] + mudtif.RasterXSize*trans[1],
trans[3] + mudtif.RasterYSize*trans[5], trans[3])
# Display the array with Matplotlib
figure(figsize = (7,7))
plt.imshow(arr, extent=extent)
plt.show()
In [33]:
# Transform Geotiff to ASCII grid for processing each sediment class
!gdal_translate -of AAIGrid /Users/sgfm/GAB/data/mud.tif ../SGFM-data/GAB/mud.asc
!gdal_translate -of AAIGrid /Users/sgfm/GAB/data/sand.tif ../SGFM-data/GAB/sand.asc
!gdal_translate -of AAIGrid /Users/sgfm/GAB/data/gravel.tif ../SGFM-data/GAB/gravel.asc
!gdal_translate -of AAIGrid /Users/sgfm/GAB/data/carbonate.tif ../SGFM-data/GAB/carbonate.asc
In [36]:
# Create the sediment declaration file that will be use later to
# allocate sediment class proportions on each grid nodes.
file = open("../SGFM-data/GAB/seddeclaration.txt ", "w")
# First write the number of sediment classes declared
file.write("6 \n")
'''
# Position of the sediment defined in the XmL file in relation to the
geological sedimentary maps above:
- Here the carbonate map is in position 1,
- the gravel map in position 2,
- the sand map in position 4 and
- the mud map in position 5
It means that there are two other parameters defined in Lecode at position 3 and 6 but
which are not defined with a geological map
'''
file.write("1 2 4 5 \n")
# Hardness of the sediment defined in the XmL file
file.write("5 3 1 1 1 1 \n")
file.close()
Perform the transformation to Lecode deposit input file.
In [37]:
import os, sys
'''
# Perform the transformation to LECODE deposit format using ocean_geol_2_dep script:
- argument 1: sediment declaration file,
- argument 2: number of declared sediment maps above,
- argument 3: position of the new sediment based on its location in the XmL input file,
- argument 4 to 4+nb_map: list of .asc files for each sediment maps,
- argument 5+nb_sed: thickness of the sedimentary layer (here 50 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('ocean_geol_2_dep ../SGFM-data/GAB/seddeclaration.txt 4 3 \
../SGFM-data/GAB/carbonate.asc ../SGFM-data/GAB/gravel.asc \
../SGFM-data/GAB/sand.asc ../SGFM-data/GAB/mud.asc \
50 ../SGFM-data/GAB/gab_6sed.dep')
Out[37]:
In []: