This is a good question, and one that I tend to get asked from time to time. First, as you've pointed out, the equation for TWI = ln(a / tan(B)), where a is the 'specific' catchment area (i.e. the upslope inflowing area normalized for a measure of contour length) and B is the slope gradient, in radians, at the grid cell. As you correctly pointed out TWI will ...
With Python, you can access raster statistics using the Python GDAL/OGR API.
from osgeo import gdal
# open raster and choose band to find min, max
raster = r'C:\path\to\your\geotiff.tif'
gtif = gdal.Open(raster)
srcband = gtif.GetRasterBand(1)
# Get raster statistics
stats = srcband.GetStatistics(True, True)
# Print the min, max, mean, stdev based on ...
Defining ridges vs hill/mountain tops is pretty scale-dependent. Jeff Jenness covers conceptually how to model topographic landforms in his article Some Thoughts on Analyzing Topographic Habitat Characteristics. If you poke around on his website, you can find his poster on this as well, under ArcGIS tools > Land Facet Corridor Designer. (Link is here)
With bash alone, you can use :
gdalinfo -mm input.tif
It returns a range of infos within which is the string Computed Min/Max=-425.000,8771.000, for my Eurasian raster.
Some cleanup and you get your vertical min/max variables:
$zMin=`gdalinfo -mm ./input.tif | sed -ne 's/.*Computed Min\/Max=//p'| tr -d ' ' | cut -d "," -f 1 | cut -d . -f 1`
From a theoretical point of view depression filling only has one solution, although there can be numerous ways of coming to that solution, which is why there are so many different depression filling algorithms. Therefore, theoretically a DEM that is filled with either the Planchon and Darboux or the Wang and Liu, or any of the other depression filling ...
The best all round tool here is a raster calculator.
gdal_calc is a GDAL raster calculator implemented in Python here, with some examples here.
If you e.g. wants to keep values above +50:
gdal_calc.py -A input.tif --outfile=result.tif --calc="A*(A>50)" --NoDataValue=0
You can specify several files -A to -Z, where each of them get a corresponding ...
Creating watersheds should help you locate both ridges and hill top. Then, I would define a hill top as a local maximum, while a point on a ridge is not the maximum (there is one other point higher or equal to this point). You can identify local maxima using the focal statistic tool.
another way to look at the problem is to analyse at the opposite of your ...
Hypsography concerns the land's elevation, altitude or height above sea-level or some other reference surface. (Hypso is derived from the Greek Ύψος for height.)
Topography concerns physical and cultural features of the land and so includes hypsography, hydrology, the built environment, major boundaries, communication channels, etc. (Topo is derived from ...
It's super-easy in QGIS 3.0:
Run the "Set Z Value" Processing algorithm
Click the button on the right of "Z Value", and select Field -> "DYBDE".
Run the algorithm. The z values for the geometry's vertices will be set to the value from the DYBDE field.
In case you have the values of depths and you want to get elevation values with negative number for ...
1) you read the data (x,y,z) from text files, shapefiles, etc., with Python only or with different Python modules (pandas, csv, Python - Excel, Fiona, Pyshp, osgeo:OGR or ...)
2) you can plot the points
in 2D: with matplotlib
in 3D: with visvis or matplotlib for example
3) you can compute contour lines (and other things...)
4) you interpolate the ...
There is more than one source depending on the area of interest around the globe.
Google Maps is the terrain layer, which provides a shaded relief (aka
hillshade) view of the topography derived from a digital elevation
model. Google has done a nice job generating a visually pleasing
terrain layer, and we use it for all of our Google Maps-based
I believe it is actually used to generate a Magnetic north line, using the reference point P on the bottom edge of the map. At least that's how I interpreted the statement in the declination legend (shown below) which seems to refer to the scale you mentioned.
I picked a sheet a random. At the time the map was done the angle was 31.25 (+ 2.5 degrees?) ...
Given your crop.tif raster layer, you can filter its pixels whose elevation is above the threshold (50 m) using gdal_calc.py:
gdal_calc.py -A crop.tif --outfile=result.tif --calc="A>=50" --NoDataValue=0
The result.tif will be made of 1 where the condition is satisfied, 0 otherwise.
Then, it will be possible to vectorize result.tif using gdal_polygonize....
SRTM (Shuttle Radar Topography Mission) was a shuttle mission, no satellite involved. But essentially the satellites do not cross the poles.
In a sun-synchronous orbit, which most imaging satellites are in, you get a pattern like:
This is great because it means that the orbit can be timed and most parts of the Earth get covered at around noon, getting good ...
Is Python an option?
Use RasterIO (a Python GDAL/ numpy bridge) to load the raster to NumPy array, then use numpy.amax() to find the maximum value, followed by numpy.where() to find the row/column indices, then calculate the lat and lon from the raster extents.
In the world of hydrology and geomorphology, there is indeed a metric that we use to classify/quantify the "curviness" of a river......sinuosity.
Sinuosity is simply a measure of the actual path length of the river divided by the shortest path length (straight line distance).
So, you could measure the sinuosity of the river as a whole (actual path length ...
MGRS is based entirely on UTM (Universal Tranverse Mercator) projected Coordinate Reference Systems. The first two numbers of an 8 digit grid, for instance, are the same as the UTM zone (these divide the world into 59 strips, each running from the equator to one of the poles). The difference is in the lettering that MGRS uses - the letter after the UTM ...
You get this kind of picture because every file has a different range of gray values, and QGIS scales the colours between min and max seperately. To solve this:
Create a virtual raster on all your files using GDAL, or from the QGIS menue.
Load that instead of the individual files as one single layer.
1. Get the pixels value: gdal's gdallocationinfo allow to access a pixel's value.
The gdallocationinfo utility provide a mechanism to query information
about a pixel given it's location in one of a variety of coordinate
systems. Several reporting options are provided.
$ gdallocationinfo crop.tif 50 50
The question (as clarified in a comment) asks how to
remove local slope to calculate relative ruggedness.
There is a simple way to do this. It relies on computing the slope using the same local data as the ruggedness (which usually is a 3 by 3 square neighborhood). I recall verifying that ArcGIS computes slope (s) and aspect in exactly this manner: its ...
In QGIS you can use the Interpolation plugin (I think it is installed by default, if not it is a standard plugin installed through Plugin Manager) found in Raster -> Interpolation.
With that plugin you can take your vector layer of points and turn them into either a TIN or a surface model. You can either pull the elevation from a field in an attribute ...
One easy way of doing this would be to inverse your DEM by multiplying it by negative one (Raster Calculator) then running the Fill tool on the inverted DEM. Finally, subtract the filled DEM from the inverted and multiply by negative one again (putting it back to the original scale). This will effectively turn peaks into depressions and find the spill height ...
It is important to remember that when computing hillshading, you need to have an illumination source. Using the sun as an illumination source may mean that a cell is shaded at noon, when the sun is directly overhead, but not at 4:00 p.m. Without an illumination point, your example seems more like a color coded slope map.
ESRI calculates illumination of ...
I think the problem is related to how you change the projection of your data. If the DEM data was originally in geographic projection (WGS84), and not projected to UTM zone 10S, then you need to reproject the original DEM image (the WGS84) again in a proper way. When you reproject (Warp) the DEM from geographic to UTM projection, it is better that you choose ...
If you are using QGIS 3.10 or later, you have better control over the GDAL tools - See the visual change log.
It means, GDAL-Contour tool can take direct input from users to which elevation (height) it should create contour lines.
Run the tool (in the Processing Toolbox > GDAL > Raster Extraction) and find Advanced parameters > Additional command-...
OK- I figured out a fix for this. Susan Jones has a script http://arcscripts.esri.com/details.asp?dbid=16055 that works the way I was hoping the Bearing Distance to Line tool would work.
The output from the script were lines radiating at varying angles and distances from my base coordinate (datum).
Then I used Feature Vertices to Point to add an X,Y ...
I will attempt to answer my own question - dun dun dun.
I used SAGA GIS to examine the differences in filled watersheds using their Planchon and Darboux (PD) based filling tool ( and their Wang and Liu (WL) based filling tool for 6 different watersheds. (Here I only show case two sets of results - they were similar across all 6 watersheds) I say "based", ...
the difference between roughness and slope is a question of scale. I recommend that you think about the resolution of your raster at which you observe the slope but you don't see the roughness anymore, then you can smooth your surface (e.g. using a low pass, a mean filter or some spline) at this resolution. This will yield a new surface with zero roughness, ...