Not sure if that is an issue so feel free to close if it is just the way it is.
I am computing TWI over a pretty large area using 32 cores machine with a ton of RAM. My DEM are around 1.4gb and I am working on SSD. When I compute the Dinf flow direction, I get a value of 295392041 flats to resolve which even with a 32 cores machine, takes more than 24h.
I have the feeling that the algorithm see my no_data value as actual values as, if I crop my DEM but keep the same extent, I get even more flats to resolve.
Is there a preferred no_data value? I use -9999 as seen above. Here is a gdalinfo output of one of my raster file and my workflow at the bottom.
Am I doing something wrong?
J$ gdalinfo /Users/julienschroder/Downloads/CTIAnalysis_RasterLayers_26OCT2015/AKPCTR_DEM.tif
Driver: GTiff/GeoTIFF
Files: /Users/.../AKPCTR_DEM.tif
/Users/.../AKPCTR_DEM.tif.ovr
/Users/.../AKPCTR_DEM.tfw
/Users/j/Users/.../AKPCTR_DEM.tif.aux.xml
Size is 20332, 18577
Coordinate System is:
PROJCS["NAD_1983_StatePlane_Alaska_1_FIPS_5001",
GEOGCS["NAD83",
DATUM["North_American_Datum_1983",
SPHEROID["GRS 1980",6378137,298.2572221010002,
AUTHORITY["EPSG","7019"]],
AUTHORITY["EPSG","6269"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4269"]],
PROJECTION["Hotine_Oblique_Mercator"],
PARAMETER["latitude_of_center",57],
PARAMETER["longitude_of_center",-133.6666666666667],
PARAMETER["azimuth",-36.86989764583333],
PARAMETER["rectified_grid_angle",0],
PARAMETER["scale_factor",0.9999],
PARAMETER["false_easting",5000000],
PARAMETER["false_northing",-5000000],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]
Origin = (303571.061199599993415,1117618.294335020007566)
Pixel Size = (49.999999999999993,-50.000000000000000)
Metadata:
AREA_OR_POINT=Area
DataType=Generic
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 303571.061, 1117618.294) (162d14'22.18"W, 23d20'23.00"N)
Lower Left ( 303571.061, 188768.294) (156d53'39.64"W, 19d 8'34.25"N)
Upper Right ( 1320171.061, 1117618.294) (157d10'48.89"W, 29d17'48.43"N)
Lower Right ( 1320171.061, 188768.294) (151d15'29.26"W, 24d29'35.05"N)
Center ( 811871.061, 653193.294) (156d54'44.66"W, 24d 6'27.45"N)
Band 1 Block=20332x128 Type=Float32, ColorInterp=Gray
NoData Value=-9999
Overviews: 10166x9289, 5083x4645, 2542x2323, 1271x1162, 636x581, 318x291, 159x146
# some setup
import os, glob, rasterio
import numpy as np
#set path for Taudem EXE
os.chdir('/home/UA/truc/TauDEM/')
# Compute pit filling Taudem algorithm
os.system('mpirun -n 32 pitremove -z /workspace/Shared/Users/truc/CTI_ref/DEM_cropped.tif -fel /fast_scratch/truc/CTI_processing_cropped/AKPCTR_DEMfel.tif')
# Mask the glaciers
mask_rasters('/fast_scratch/truc/CTI_processing_cropped/AKPCTR_DEMfel.tif', ['/workspace/Shared/Users/truc/CTI_ref/RGI_v5_AKPCTR_NewExtent.tif', '/workspace/Shared/Users/truc/CTI_ref/Seam3000_NewExtent.tif'],1,-9999,'/fast_scratch/truc/CTI_processing_cropped/AKPCTR_NoIceSeam.tif')
#Run the flow direction taudem algorithm
os.system('mpirun -n 32 dinfflowdir -fel /fast_scratch/truc/CTI_processing_cropped/AKPCTR_NoIceSeam.tif -slp /fast_scratch/truc/CTI_processing_cropped/AKPCTR_NoIceSeam_slp.tif -ang /fast_scratch/truc/CTI_processing_cropped/AKPCTR_NoIceSeam_ang.tif')
#Reclass the slope raster from 0 to 0.0001
slope_reclasser('/fast_scratch/truc/CTI_processing_cropped/AKPCTR_NoIceSeam_slp.tif','/fast_scratch/truc/CTI_processing_cropped/AKPCTR_NoIceSeam_No0slp.tif')
#Run the contributing area taudem algorithm
os.system('mpirun -n 32 areadinf -ang /fast_scratch/truc/CTI_processing_cropped/AKPCTR_NoIceSeam_ang.tif -sca /fast_scratch/truc/CTI_processing_cropped/AKPCTR_NoIceSeam_sca.tif')
#Run the actual TWI calculation
os.system('mpirun -n 32 twi -sca /fast_scratch/truc/CTI_processing_cropped/AKPCTR_NoIceSeam_sca.tif -slp /fast_scratch/truc/CTI_processing_cropped/AKPCTR_NoIceSeam_No0Slp.tif -twi /fast_scratch/truc/CTI_processing_cropped/CTI_NoIceSeam_1.tif')
#Clip the CTI raster by land mask
mask_rasters('/fast_scratch/truc/CTI_processing_cropped/CTI_NoIceSeam_1.tif', '/fast_scratch/truc/CTI_processing_cropped/AKPCTR_WS_Outline_4kmBuff.tif',1,-9999,'/fast_scratch/truc/CTI_processing_cropped/CTI_NoIceSeam.tif')