Saturday, April 8, 2017

A pair of short courses on "Geoscience Investigations of Point Clouds" and "Advancing understanding of geomorphology with topographic analysis": mid June, 2017 at Potsdam University, Golm

Two short courses are scheduled for mid June at Potsdam University. The short courses are independent of each other; however, the topics are related and probably address a similar audience.

Geoscience investigations of point clouds, June 7-9, 2017. Instructors B. Bookhagen, R. Arrowsmith, M. Isenburg, C. Crosby.
This course will explore the acquisition, post-processing, and classification of point clouds derived from airborne and terrestrial lidar scanners and structure from motion (SfM) photogrammetry from drones. The course will take place at campus Golm (UP) and includes one day of field-data collection and two days of data post-processing and analysis.
The application is here: The course is sponsored by Geo-X.
Here is the flyer: PDF for more details.

Advancing understanding of geomorphology with topographic analysis emphasizing high resolution topography, June 12-15, 2017. Instructors R. Arrowsmith, W. Schwanghart, C. Crosby, B. Bookhagen.
This course will focus on advanced understanding of geomorphology with topographic analysis emphasizing high-resolution topography. The course will take place at campus Golm (UP) and includes theoretical background and analysis of digital topography using TopoToolbox in a Matlab environment. The course is sponsored by StRATEGy.
Here is the flyer: PDF for more details.

Thursday, March 16, 2017

Mapping landforms with applications to geomorphology and earthquake geology EXERCISE

I wanted to share a project I have developed on and off for about 10 years. It is a classroom exercise for Mapping landforms with applications to geomorphology and earthquake geology. So far, it has an example for strike-slip faults (Wallace Creek along the San Andreas Fault), and blind thrust faults (Wheeler Ridge in Southern California). I have an intention to add a normal fault example but have not finished it yet.

The basic idea is that the exercises could be done "analog"--that is on paper in a classroom. They emphasize some simple morphologic and geomorphic mapping (but on high resolution topography base maps from lidar data collected by NCALM and available for download from OpenTopography), and then provide landform ages so that students can calculate slip rate or surface uplift rates. I am not sure they are so well explained so any feed back is welcome.

Exercise material:

Sunday, March 12, 2017

DEM grid size specification and learning a bit of TopoToolbox

I have been wanting to learn TopoToolbox for a while, and last week I had a chance to get started. Spending time at the Universität Potsdam, Institut für Erd- und Umweltwissenschaften, I have met with Wolfgang Schwanghart and Dirk Scherler occasionally. They are coauthors of TopoToolbox.

My ASU colleagues Adam Forte and Kelin Whipple have been using TopoToolbox increasingly over the last few years. Recently Kelin had a question as to why some DEMs grids were causing an error in TopoToolbox and why some were not. I dug into it and here is what I found. Projected DEMs should have EXACTLY the same x and y cell sizes and expressed in precise values (nearest meter or 10th or 100th of a meter).

DEMs (such as SRTM 30 meter) can be downloaded as geotiff formats and with geographic coordinate systems from sources such as OpenTopography: select Global Data tab.

For example, I have chosen a piece of data from Java along the Cimandiri fault (see Marliyani, et al., 2016 for example).

To be useful for most geomorphic analyses, the data should be projected to UTM so that the horizontal and vertical units are the same. Often, the projection is done in ArcGIS. And, usually it seems fine to let Arc determine the cell size automatically.

Notice how ArcGIS decided that the resolution for the cells should be: 30.8462728281739.
Here I set the cell size to exactly 30 m

Turning to TopoToolbox, as we compute drainage network properties, including contributing area, there is a check that the cell sizes are the same (inside the GRIDobj.m function):

if abs(abs(dx)-abs(dy))>1e-9;
    'The resolution in x- and y-direction must be the same');
The 1e-9 is a somewhat arbitrarily small number which one would think would not cause a problem.

I looked into this problem tracking along with the TopoToolbox processing and some of the MATLAB built in tools.

  1. Read the geotiff using MATLAB's geotiffread:
    [demdata, R] = geotiffread(filename);
    geotiffreaddifference = R.CellExtentInWorldX-R.CellExtentInWorldY;
    The R object has a number of values including the cell size of the geotiff:
    R.CellExtentInWorldX= 30.79776426908749800
    R.CellExtentInWorldY= 30.79776426908791000
    R.CellExtentInWorldX-R.CellExtentInWorldY= -4.12115e-13
  2. Read the geotiff into the DEM object in TopoToolbox and do a similar check as above (refmat is similar to R):
    DEM = GRIDobj(filename);
    refmatdifference = abs(DEM.refmat(2,1))-abs(DEM.refmat(1,2));
    Again we see the same very small grid size difference:
    abs(DEM.refmat(2,1))-abs(DEM.refmat(1,2))= -4.12115e-13
  3. The problem comes in a calculation that builds out the x and y vectors that includes a cumulative multiplication from the grid sizes inside a TopoToolbox function refmat2XY:
    nrrows = siz(1);
    nrcols = siz(2);
    x = [ones(nrcols,1) (1:nrcols)' ones(nrcols,1)]*R;
    x = x(:,1)';
    y = [(1:nrrows)' ones(nrrows,2)]*R;
    y = y(:,2);
    R is the refmat object. But, here we get the error:
    dx=x(1)-x(2)= -30.79776426916942000
    dy=y(1)-y(2)= 30.79776426777243600
    abs(dx)-abs(dy)= 1.39698e-09
    Which then would fail the GRIDobj test.
  4. If I run the same set of checks on the file I projected with exactly 30 m grid cell size (recall above), I don't get the problem:
    R.CellExtentInWorldX= 30.00000000000000000
    R.CellExtentInWorldY= 30.00000000000000000
    R.CellExtentInWorldX-R.CellExtentInWorldY= 0
    refmat difference as produced from GRIDobj:
    abs(DEM.refmat(2,1))-abs(DEM.refmat(1,2))= 0
    Difference after GRIDobj2mat:
    dx=x(1)-x(2)= -30.00000000000000000
    dy=y(1)-y(2)= 30.00000000000000000
    abs(dx)-abs(dy)= 0

I think what is happening is a bit of roundoff error (see this link for detailed discussion: What Every Computer Scientist Should Know About Floating-Point Arithmetic, by David Goldberg): "Therefore the result of a floating-point calculation must often be rounded in order to fit back into its finite representation." MATLAB does all of its calculations by default in double precision, so we should have access to 16 decimal digits (e.g., So, why we loose precision up to 10-9 at times is a little surprising to me, but it can obviously happen.

Thus there are a couple of workarounds:

  1. Comment out the dx and dy check in GRIDobj:
    %if abs(abs(dx)-abs(dy))>1e-9;
    %    error('TopoToolbox:GRIDobj',...
    %    'The resolution in x- and y-direction must be the same');
    %    end
    Or make the threshold larger.
  2. Explicitly set the grid resolution to a round number when projecting in ArcGIS or other software.
  3. Use the reproject2utm.m command to project dems before using other TopoToolbox operations. It sets the resolution to be exactly the same for x and y

Here is a script to run these calculations: DEM_cell_size_Script.m

Thanks to Chris Crosby, Benjamin Gross, Wolfgang Schwanghart, Kelin Whipple, and Mike Zoldak for discussions.

Sunday, March 5, 2017

Field trip in NW Himalaya (Himachal Pradesh)

I just returned from a very interesting and pleasant field trip in Himachal Pradesh (India) to examine the major elements of the frontal portion of the NW Himalaya. The trip was lead by Dr. Rasmus Thiede of the University of Potsdam (UP). This was my third trip to the Himalayas (the first in 2001 to the Sutlej and Spiti Rivers with UP colleagues including Rasmus and Bodo Bookhagen when they were starting their Ph.D.s under the supervision of Professor Manfred Strecker; the second was in 2010 with Dr. Wendy Bohon in the Ladakh region to contribute to her work along the Karakoram Fault in the Pangong Range area). We joined two UP students (Katharina Kretzschmar and Markus Nennewitz) who had been working in the area for two weeks already, focusing on developing a structural profile across the principal faults of the NW Himalaya in the region. I overlapped for a day with Dr. Saptarshi Dey who worked on his Ph.D. with Rasmus in the Kangra reentrant (see Dey, et al., 2016a,b). And, it was a pleasure to meet Professor Vikrant Jain from IIT Gandhinagar and his Ph.D. students. I did not (yet) get their last names but they were Ramindran, Sonam, and Pritha. Tashi Gyatseo is a long-time guide and friend of Rasmus who was chief of logistics for the trip. It was fun to meet everyone and to discuss the geology and geomorphology as well as learn a bit more of eachother's cultures.

Along with strengthening the collaborations among the groups, there were a couple of scientific targets of the trip. They build to a large degree from the Ph.D. work of Saptarshi Dey and focus on:

  1. The role of climate modulation of the sediment supply and transport capacity along the major river systems and the development (and removal) of fluvial fills and terraces.
  2. The activity of the major faults of the system, including the Himalayan Frontal Thrust and out of sequence faulting (activity higher in the tectonic wedge) possibly influenced by load variation from changes in sediment storage in the orogenic wedge.
So, we looked at a lot of interesting fluvial terraces mostly along the Ravi River system as well as reviewed the major structural elements of the system: HFT, Siwaliks, Main Boundary Thrust (MBT), and the Main Central Thrust (MCT).

Overview figure of the region from Gavillot, et al., 2016. The box shows the location of our trip and the map shown below.

Overview of the field trip with the orange-red line showing most of our tracks and the yellow points my main obversation locations. I am using the GaiaGPS app for my mapping. It seems ok but I might look for a more geologically oriented app for the next trip.

Panorama and annotated overlay showing the major structures and bounded rock units


View of the nice mountains (Dhaula Dhar Range) from the Kangra airport--a nice way to start the trip!
Katharina, Markus, and Rasmus with the impressive Dhaula Dhar Range behind as seen from the Kangra area.
Most of the group at a stop along the Main Boundary Thrust where a spring may indicate increased permeability along the fault zone (yellow-white material at left is travertine).
Upper Siwaliks at sunset
Katharina in the MCT mylonites.
Chambo Town--built on an important probably 10 ka terrace. We spent a pleasant cool night there.
Pretty agricultural terraces on a high terrace south of Chamba.
Nice terraces near the Chemera Reservoir.
Stair step terraces in a tributary a few km below the Chemera Reservoir.
Rasmus is a great teacher and it was wonderful to learn from him. Here he is near some fine grained schists discussing their reflection of the overall deformation in the MCT zone.

References cited
Dey, S., Thiede, R., Schildgen, T., Wittmann, H., Bookhagen, B., Scherler, D., and Strecker, M. Holocene internal shortening within northwest Sub-Himalaya: Out-of-sequence faulting of the Jwalamukhi Thrust, India: Tectonics, 2016a.

Dey, S., Thiede, R. C., Schildgen, T. F., Wittmann, H., Bookhagen, B., Scherler, D., … Strecker, M. R. (2016b). Climate-driven sediment aggradation and incision since the late Pleistocene in the NW Himalaya, India. Earth and Planetary Science Letters, 449, 321–331.

Gavillot, Y., Meigs, A., Yule, D., Heermance, R., Rittenour, T., Madugo, C., & Malik, M. (2016). Shortening rate and Holocene surface rupture on the Riasi fault system in the Kashmir Himalaya: Active thrusting within the Northwest Himalayan orogenic wedge. Bulletin of the Geological Society of America, 128(7), 1070–1094.

Tuesday, December 20, 2016

Landers earthquake fault scarp Structure from motion

I made a movie of structure from motion high resolution view of 1992 Landers California earthquake fault scarp. Video starts with 2012 hillshade (Johnson, et al., 2014; available from OpenTopography here: link) in Google Earth to show location and then to a ground based set of photographs (see blue rectangles as focal planes) visualized in Agisoft Photoscan.

I am pretty pleased that the ground-based model worked so well. Now we can move forward with fine scale alignment with earlier topographic point clouds and compute differences over the 25 years since the earthquake--a project I have worked on with Dallas Rhodes for many years (see Arrowsmith and Rhodes, 1994 and also Haddad, et al., 2012).

See also these posts:

  • 2015 Anniversary of 1992 Landers California earthquake
  • SfM mapping--also has an orthophoto kmz of the Johnson, et al., 2014 data
  • Structure from Motion micro documentary from Merri Lisa Trigilio
  • References:

    • Arrowsmith, J. R., & Rhodes, D. D. (1994). Original forms and initial modifications of the Galway Lake Road scarp formed along the Emerson Fault during the 28 June 1992 Landers, California, earthquake. Bulletin - Seismological Society of America, 84.
    • Haddad, D. E., Akciz, S. O., Arrowsmith, J. R., Rhodes, D. D., Oldow, J. S., Zielke, O., … Shilpakar, P. (2012). Applications of airborne and terrestrial laser scanning to paleoseismology. Geosphere, 8(4).
    • Johnson, K., Nissen, E., Saripalli, S., Arrowsmith, J. R., McGarey, P., Scharer, K., … Blisniuk, K. (2014). Rapid mapping of ultrafine fault zone topography with structure from motion. Geosphere, 10(5).

    Friday, December 9, 2016

    New report: NASA Challenges and Opportunities for Research in Earth Surface and Interiors

    The new report: NASA Challenges and Opportunities for Research in Earth Surface and Interiors has just been released officially. The main link to download is here: PDF.

    I was honored to be on the writing team and contributed to the surface process, human activities, topography, increasingly interconnected world, and professional development portions. We were charged with revisiting and updating the 2002 Solid Earth Science Working Group report “Living on a Restless Planet” (the SESWG Report). The update follows the SESWG framework and updates on many of the science and technology topics and will help to chart NASA Earth Sciences and Interiors priorities.

    Thanks to the rest of the committee, our co chairs James Davis and Louise Kellogg, and Ben Phillips from NASA.

    Sunday, November 13, 2016

    Accumulating links and interpretation for earthquake M7.8 - 53km NNE of Amberley, New Zealand

    I am accumulating some links and thoughts with respect to today's earthquake M7.8 - 53km NNE of Amberley, New Zealand. I hope that the damage won't be too severe and I am sending positive energy to the people there. Magnitude and depth were increasing as the seismologists reviewed the seismograms. It was felt widely across New Zealand (Felt reports:

    >15,000 felt reports as of 17:14 UTC

    Here is a sketch of the NZ tectonic setting (from

    March 6, 2017 update:
    GNS updates including offshore faulting from high resolution bathymetry

    Jan. 1-4 updates:
    GNS field work blog
    New Zealand Geographic article
    Before and after images

    Dec. 5 updates:
    Kekerengu Fault field work blog

    Nov. 28 updates:
    Kekerengu Fault and prior trenching (GNS)

    Nov. 27 updates:
    Repeat satellite imagery showing big shift! (Chris Milliner)
    Geospatial Information Authority of Japan interesting differencing
    EarthJay blog post

    Nov. 25 updates:
    Farm track ruptured by fault drone video
    The Kekerengu Fault rupture pictures from Julian's Rock and Ice Blog

    Nov. 23 updates:
    NPR report with some of the amazing drone videos and other commentary
    Temblor site with some of the amazing drone videos and other commentary

    Nov. 21 updates:
    Amazing drone video of the Kekerengu Fault rupture
    Kaikora coastal uplift
    Skateboarders making the most of the rupture....
    Rob Langridge GNS Science talks earthquakes with JR

    Nov. 19 updates:

    GeoNET: viewing the earthquake from space
    GeoNET: measuring the earthquake with GPS
    Fagereng, Complex Earthquake Raises Complex Questions (EOS)
    Trembling Earth blog entry 2

    Nov. 17 updates:
    Evolving magnitudes on quakestories
    Preliminary landslide mapping

    Nov. 15 and 16 updates:
    GNS Blog
    Trembling Earth blog entry 1
    COMET for interferogram and interpretation
    Kaikoura district faults report by GNS
    Temblor interpretation
    Commentary in Spinoff--GNS scientists interviewed
    Science commentary
    Duffy and Quigley commentary in the Conversation
    Tsunami damage Banks Peninsula
    Digital Globe images of surface rupture from Ryan Gold (USGS)
    Preliminary Sentinel-1 interferogram
    Yahoo News Seafloor uplift article

    Nov. 14 updates:
    Geonet what we know so far
    Geonet update
    IRIS Recent Earthquake Teachable moment
    IRIS Special Event Site: South Island, New Zealand
    Surface rupture, landslides, damage:

    Nov. 13 Interpretion: The magnitude and depth (as long as it does not get much deeper) could be consistent with a rupture on the Hope Fault. The location and to some degree the focal mechanism would be more consistent with one of the thrust faults further to the west. As information has flowed in, I get a sense that the rupture was deeper and mostly underneath the main crustal faults and approaching the subduction interface below.

    Reports are coming in that there is a tsunami that was generated and has hit the coast with heights of a few meters in places.

    Interesting to see aftershocks aligned along the faults to the northeast. USGS tectonic summary suggested some slip along the megathrust. I wonder if it is in a sort of accretionary complex above the megathrust and transitioning into to the shearing plate boundary (hence the depth and steeper dip and oblique slip focal mechanism). USGS finite fault and source time function show deeper slip almost 100 km northeast (and 60 seconds) from the hypocenter (making this a complex earthquake). Peak slip is ~ 4m (oblique) at about 25 km down dip distance from the surface along a 38 degree to the NW dipping surface. See also IRIS backprojection results. But that is necessarily a simple model. The Cape Campbell GPS station moved 2 m E, 1 m N, and 1 m up apparently (scroll down) and that would imply shallower or greater or more complex slip than the simple inversion implies. It makes me wonder what surface rupture will look like. USGS finite fault source now indicates strong ground motions and likely seafloor uplift zone (Kaikora to almost Wellington; see the image below).

    Here is a mashup of GNS faults with the USGS location, shaking estimates, additional earthquakes, and focal mechanism as of 16:23 UTC:

    Links (some coming from twitter feed)

    News links: