Smathermather's Weblog

Remote Sensing, GIS, Ecology, and Oddball Techniques

Archive for the ‘Ecology’ Category

Gorilla Food Plants Biomass and Ranging

Posted by smathermather on July 4, 2017

The last two days, I have been working with Olivier Jean Leonce Manzi here at Karisoke Research Center on the question of the relationship between gorilla food plants biomass and the ranging patterns of gorillas outside Volcanoes National Park (VNP) in Rwanda.

A view from the edge of VNP looking toward Mounts Mgahinga and Muhabura

A view from the edge of VNP looking toward Mounts Mgahinga and Muhabura

Even though Volcanoes National Park is set aside and gorillas thrive inside, they don’t strictly stay within the bounds of the park — often leaving the park to forage in the adjacent farming communities. This tendency is an opportunity for human / wildlife conflict, and thus it is important to understand the relationship between gorilla browsing habits outside the park and their food resources in these farming communities.

Among the food plants in the study, eucalyptus, a non-native, is planted widely in Rwanda as a fast growing timber source. Gorillas like to climb up eucalyptus, stripping the outer bark, and use their teeth to scrape the inner bark from the tree. For larger trees, this seems to have little effect. For smaller trees, it can kill the tree easily, either by the weight of the gorilla on the tree or the damage to the vascular system of the tree.

 

Eucalyptus trees planted near VNP. Eating the inner bark of Eucalyptus is a favorite food of gorillas. For larger trees, this seems to have no effect. For smaller trees, this can maim or kill the non-native tree plantings.

Eucalyptus trees planted near VNP. Eating the inner bark of Eucalyptus is a favorite food of gorillas. For larger trees, this seems to have no effect. For smaller trees, this can maim or kill the non-native tree plantings.

Bamboo, native to the region, is widely planted outside VNP (and grows in forests inside VNP). As it sprouts during rainy season, it is a favorite food of gorillas both inside and outside the park.

Small bamboo stand outside VNP.

Small bamboo stand outside VNP.

Olivier’s study area is in the farming area along the base of Karisimbi and Bisoke Volcanoes just outside VNP.  (Quick aside: it was a combination of Karisimbi and Bisoke names that formed Karisoke Research Center’s name).

Map of study area near Karisimbi and Bisoke Volcanoes outside Volcanoes National Park.

Map of study area near Karisimbi and Bisoke Volcanoes outside Volcanoes National Park.

The two datasets we want to compare for the study are counts of gorillas ranging outside the park and transects of gorilla plant food biomass, also outside the park.

Gorilla from the Amahoro (Peace) group in a bamboo stand during rainy season. (Amahoro group is not part of this study.)

Gorilla from the Amahoro (Peace) group in a bamboo stand during rainy season. (Amahoro group is not part of this study.)

 

study_area_inset

Map of gorilla sightings outside VNP within our study area.

The measurements that Olivier made on gorilla plant food biomass is that of herbs (overall), trees (overall), eucalyptus, bamboo, and rubus (a group commonly known as raspberries, blackberries, dewberries, etc.) — favored foods for gorillas especially outside the park.

vegetation_plots.png

Vegetation plot locations along border of VNP

For the purposes of analysis, Olivier grouped the vegetation plots and gorilla counts into 18 approximately equal zones. Gorilla counts and vegetation plots were summarized per zone and compared.

vegetation_plot_zones.png

I remember his project as one of the more difficult ones to wrap my head around when I was here in December. In summary, here’s the problem that we wanted to solve: given a dataset of gorilla food biomass outside the park, and counts of gorillas outside the park, can we establish a relationship between the two. It should be an easy problem: on one side of the equation, we have count data, on the other ratio data. My first thought would be to apply a Poisson approach.

Gorilla.count ~ Herbs.biomass + Tree.biomass + Eucalyptus.biomass + Bamboo.biomass + Rubus.biomass

But we quickly run into a problem: our data are not independent in time or in space, and so a simple Poisson approach is likely not appropriate.

For as much work as I’ve done in the geospatial space, I have done precious little with spatial statistics, so this posed a conundrum to me. I can’t remember now if it was through much googling, or the great stewardship of Dr. Patrick Lorch, the research manager at my institution, that we settled upon using a Markov chain Monte Carlo (MCMC) approach to the problem. The advantage to this approach is that we can do the analysis using a Poisson distribution without regard to autocorrelation. In the R statistical package, this is an easy analysis to set up.

Finally, we’ll close with a picture of Olivier working on his literature review for the methods section:

olivier.jpg

Posted in Ecology, Gorillas, Karisoke, National Park, R | Tagged: , , , , , , | 1 Comment »

A little Gorilla Time

Posted by smathermather on June 12, 2017

I miss my mountain gorilla friends in Rwanda. Let’s write a little more code to support them. I’ll be visiting Karisoke again next week, so it seems timely to post a little more code (HT Jean Pierre Samedi Mucyo for working with me on this one).

The problem today is simple — given a time series of gorilla locations and dates, can we calculate rate of travel using PostgreSQL? Well, of course we can. We can do anything in Postgre.

We have two tricks here:

  1. The first is to order our data so we can just compare one row to the next.
  2. Once we do that, we need simply to use PostGIS to calculate distance, and ordinary time functions from Postgres to calculate time difference.

This is my first use of WITH RECURSIVE, and it’s probably unnecessary (could be replaced with windowing functions), but I was very proud to finally get over my fear of WITH RECURSIVE . (We actually use one windowing function in our prep of the data. But there we are… ).

snip

For the record, WITH RECURSIVE isn’t recursive, but it is useful here in allowing us to compare the current row with the previous.

Posted in Database, Ecology, Gorillas, Karisoke, National Park, PostgreSQL, SQL | Tagged: , , , , | Leave a Comment »

Uninformed botanical musings

Posted by smathermather on September 28, 2015

I had the good pleasure of attending FOSS4G Seoul. One of the organizers (Heegu Park) early on told me, in response to the workshop I planned, something to the effect of “Whatever you need, Steve, ask for it. Nothing is impossible.”  The organizers truly were capable of fulfilling any request.  More on that later.

Last time I was in Seoul, I took lots of pictures. This time, so few, I’m afraid. But I took a few. I need some help with botanical sleuthing.

 Jewelweed, Impatiens capensis, Geauga County, Ohio

Jewelweed, Impatiens capensis, Geauga County, Ohio

 

There is a flower native to the Eastern United States called jewelweed. Jewelweed is no exciting flower, but common in moist places, useful for treating poison ivy and other skin ailments. It’s medicinal and is among the first plants that I learned in walking in the woods in Southeast Michigan and Northwest Ohio. In the floodplains of Halfway Creek near my childhood home, it was common to the point of being weedy. (The above photo is from Geauga County in Northeast Ohio).
While hiking near Ongnyeobong Peak (옥녀봉 — the romanization is based on the sign at the peak, but I’m not sure it’s right) south of Seoul I saw this little impatiens in very similar habitat:

Unknown impatiens, Gwacheon City, South Korea

Unknown impatiens, Gwacheon City, South Korea

Anyone know it’s common or botanical name?

 

Location of unknown impatiens

Location of unknown impatiens

 

Posted in Conferences, Ecology, FOSS4G, FOSS4G2015 | Tagged: , , , | Leave a Comment »

Landscape Position using GDAL — PT 3

Posted by smathermather on November 25, 2014

More landscape position pictures — just showing riparianess. See also

https://smathermather.wordpress.com/2014/11/22/landscape-position-using-gdal/

and

https://smathermather.wordpress.com/2014/11/24/landscape-position-using-gdal-pt-2/

valleyz3 valleyz2 valleyz1 valleyz

Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL

Posted in Analysis, Ecology, GDAL, Landscape Position, Other, POV-Ray | Tagged: , , , , , , , , , , | Leave a Comment »

Landscape Position using GDAL — PT 2

Posted by smathermather on November 24, 2014

Just pretty pics today of estimated riparianess. If you prefer a bit of code, see previous post https://smathermather.wordpress.com/2014/11/22/landscape-position-using-gdal/

valleys_improved valleys_improved_1 valleys_improved_2 valleys_improved_3 valleys_improved_3.1 valleys_improved_4 valleys_improved_5

Posted in Analysis, Ecology, GDAL, Landscape Position, Other, POV-Ray | Tagged: , , , , , , , , , , | 4 Comments »

Landscape Position using GDAL

Posted by smathermather on November 22, 2014

Hat tip again to Seth Fitzsimmons. I’ve been looking for a good, easy to use smoothing algorithm for rasters. Preferably something so easy, I don’t even need to write a little python, and so efficient I can run it on 30GB+ datasets and have it complete before I get distracted again by the next shiny project (a few hours).

Seth’s solution? Downsample to a low resolution using GDAL, then sample back up to a higher resolution in order to smooth the raster. My innovation to his approach? Use Lanczos resampling to keep location static, and get a great smooth model:

Unsmoothed DEM

Unsmoothed DEM

Smoothed DEM

Smoothed DEM

Code to do this in gdal follows. “-tr” sets our resamping resolution, “-r lanczos” sets our resampling algorithm, and the “-co” flags are not strictly necessary, but I’ve got a 30GB dataset, so it helps to chop up the inside of the TIFF in little squares to optimize subsequent processing.

gdalwarp -tr 50 50 -srcnodata "0 -32767" -r lanczos  -co "BLOCKXSIZE=512" -co "BLOCKYSIZE=512" oh_leap_dem.tif oh_leap_dem_50.tif
gdalwarp -tr 10 50 -srcnodata "0 -32767" -r lanczos  -co "BLOCKXSIZE=512" -co "BLOCKYSIZE=512" oh_leap_dem_50.tif oh_leap_dem_10-50.tif

At first this excited me for cartographic reasons. We can use this to simplify contours, and then use simplified contours at different zoom levels for maps:

But, we can also use this for analyses. For example, if we difference these smoothed images with our original digital elevation model, we get a measurement of local elevation difference, the first step in establishing where valleys, ridges, and other land forms are.

# Resample to lower resolution
gdalwarp -tr 328.0523587211646 328.0523587211646 -srcnodata "0 -32767" -r lanczos  -co "BLOCKXSIZE=512" -co "BLOCKYSIZE=512" oh_leap_dem.tif oh_leap_dem_328.tif
# Upsample again to get nicely smoothed data
gdalwarp -tr 3.048293887897243 3.048293887897243 -srcnodata "0 -32767" -r lanczos  -co "BLOCKXSIZE=512" -co "BLOCKYSIZE=512" oh_leap_dem_328.tif oh_leap_dem_3-328.tif
# Merge two datasets together into single image as separate bands to ensure they are the same dimensions
# (gdal_calc, as a wrapper for numpy requires this)
gdal_merge -separate -o oh_leap_dem_3-328_m.tif oh_leap_dem.tif oh_leap_dem_3-328.tif
# And now we'll use gdal_calc to difference our elevation model with the smoothed one to get relative elevation 
gdal_calc -A oh_leap_dem_3-328_m.tif -B oh_leap_dem_3-328_m.tif --A_band=1 --B_band=2 --outfile=oh_leap_dem_lp_328.tif --calc="A-B"

So, if we want a good proxy for riparian zones, we can use a technique like this, instead of buffering our streams and rivers a fixed distance (in this case, I’ve used 4 different distances:

Map of landscape position estimated valleys in Cuyahoga County, Ohio

Map of landscape position estimated valleys in Cuyahoga County, Ohio

Pretty snazzy riparian finder. It seems to also find upland headwater wetlands (most of them historic and drained for Cuyahoga County). I am now running on 4 million acres of Ohio at a 10ft (~3 meter) resolution. It’s that efficient.

Addendum: It also finds escarpment edges, like the Portage Escarpment in the above, so it is a mix of a few landforms. Darn handy nonetheless.

Posted in Analysis, Ecology, GDAL, Landscape Position, Other, POV-Ray | Tagged: , , , , , , , , , , | 2 Comments »

Kicking the tires of PostGIS 2.0 — Testing ST_MakeValid

Posted by smathermather on April 20, 2012

The feature in PostGIS 2.0 that excited me most was not topology support, raster support, or 3D functions.  Ok, raster was near the top of my list.  But what I was really excited by was the ST_MakeValid function.  Sad, isn’t it?  Lack of vision probably– excited to try to solve recurring technical snafus in a computationally inexpensive way, rather than being more excited by the shiny new toys the PostGIS team has brought us for applying to new problems.

Never shy to try the impossible before trying the practical, I threw a really nasty 4 million acre, 3 meter vectorized landscape position dataset a colleague of mine put together.  It’s a nice dendritic, complicated mess, with polygon validity issues:

Always one for doing tests, I compared the speed with which ST_MakeValid fixed the polygons vs. Horst Duester’s cleanGeometry pgsql function, which I posted about earlier.  Sadly ST_MakeValid was much slower, and completed with error.  On the other hand, the cleanGeometry function dropped polygons, although I haven’t had time to look into why.  I don’t like my functions dropping data… .

Having posted to the PostGIS forum, Martin Davis poked a bit at the data and concluded that my problem was “self-touching polygons”, and suggested I try buffering the data a zero distance.  Ahh– it makes sense that this is the issue– it’s polygonized raster data in shapefile form.  Self-touching polygons are valid in the ESRI world but not in the simple features world (see Paul Ramsey’s presentation on PostGIS for Power Users from Foss4G 2011 for more on this).

So, ST_Buffer(the_geom, 0), and we’re good to go in half the time of the cleanGeometry function, without losing data.

Stay tuned for an ST_MakeValid wrapper which will do this trick for self-touching and self-intersecting polygons before trying to fix with ST_MakeValid.  Code below.

psql -U postgres -d test -f "C:\Program Files\PostgreSQL\9.1\share\contrib\postgis-2.0\legacy.sql"
CREATE TABLE tpi_clean AS
 SELECT gid, id, gridcode, "class name", cleanGeometry(geom) AS the_geom
 FROM tpi;

CREATE TABLE tpi_valid AS
 SELECT gid, id, gridcode, "class name", ST_MakeValid(geom) AS the_geom
 FROM tpi;
CREATE TABLE facepalm AS
  SELECT gid, id, gridcode, "class name", ST_Buffer(geom, 0) AS the_geom
  FROM tpi;

 

 

Posted in Database, Ecology, Landscape Position, PostGIS, PostgreSQL | Tagged: , , , , , , , , | Leave a Comment »

GDAL Slopes– Local Hydrologic Slope vs. the Standard Approach

Posted by smathermather on December 15, 2011

Open Source software is not, of course just about warm and fuzzies, great support, rapid development cycles, shared costs, etc., it’s also about getting your hands dirty with someone else’s code and implementing stuff more quickly and more intelligently because of it, and hopefully learning something in the process.  You don’t have to poke under the hood to drive the car, but sometimes it’s nice to anyway.

The two normal approaches to calculating slope in GIS raster tools (used in applications from ESRI’s Spatial Analyst and GDAL’s gdaldem) are Horn’s formula and Zevenbergen & Thorne’s formula, both integrated calculations of slope that do a bit of generalization in the process of calculating slope.  This is particularly useful if you are using slope for more than a localized calculations, for example in doing estimates of local light conditions.  That said, if you want to match local hydrologic slope (something which will more closely approximate a field measurement of slope), you are better off using a different formula altogether.  An article a colleague sent me: “When GIS Slope is Not What You Think” in this month’s The Forestry Source tipped me off to this, and I thought I’d do some quick and dirty digging in the GDAL code to see if I could implement something a little different.

Here is the standard Horn approach as currently implemented in GDAL:

float GDALSlopeHornAlg (float* afWin, float fDstNoDataValue, void* pData)
{
    const double radiansToDegrees = 180.0 / M_PI;
    GDALSlopeAlgData* psData = (GDALSlopeAlgData*)pData;
    double dx, dy, key;

    using std::max;
    dx = (max((afWin[0], afWin[3], afWin[3], afWin[6])) -
          max((afWin[2], afWin[5], afWin[5], afWin[8])))/psData->ewres;

    dy = ((afWin[6], afWin[7], afWin[7], afWin[8]) -
          max((afWin[0], afWin[1], afWin[1], afWin[2])))/psData->nsres;

    key = (dx * dx + dy * dy);

    if (psData->slopeFormat == 1)
        return (float) (atan(sqrt(key) / (8*psData->scale)) * radiansToDegrees);
    else
        return (float) (100*(sqrt(key) / (8*psData->scale)));
}

Borrowing from gdaldem roughness we implement this:

float GDALSlopeHydroAlg (float* afWin, float fDstNoDataValue, void* pData)

{
    // Hydrologic Slope is the max
    //  local slope btw center cell and adjacent cells

    const double radiansToDegrees = 180.0 / M_PI;
    GDALSlopeAlgData* psData = (GDALSlopeAlgData*)pData;
                double key;

    float pafLocalMin = afWin[0];
     float pafLocalMax = afWin[0];

    for ( int k = 1; k < 9; k++)
    {
        if (afWin[k] < pafLocalMin)
        {
            pafLocalMin=afWin[k];
        }

        if (afWin[k] > pafLocalMax)
        {
            pafLocalMax=afWin[k];
        }
    }

    key = afWin[4] - pafLocalMin;

    if (psData->slopeFormat == 1)
        return (float) (atan(sqrt(key) / (2*psData->scale)) * radiansToDegrees);
    else
        return (float) (100*(sqrt(key) / (2*psData->scale)));
}

Output is rougher than the original, as expected:

Local Hydrologic Slope

Zevenbergen & Thorne’s Slope

And looking at the stats we see a larger standard deviation for local hydrologic slope:


Z:\>gdalinfo Hydrologic_slope.tif
Driver: GTiff/GeoTIFF
Files: Hydrologic_slope.tif
...

Metadata:
STATISTICS_MINIMUM=0
STATISTICS_MAXIMUM=46.168109893799
STATISTICS_MEAN=10.206419514979
STATISTICS_MEDIAN=2.8041767686283e-266
STATISTICS_MODE=1.4463930749435e-307
STATISTICS_STDDEV=5.3794018830967
LAYER_TYPE=athematic

Z:\>gdalinfo Zevenbergen_Thorne.tif
Driver: GTiff/GeoTIFF
Files: Zevenbergen_Thorne.tif
...
Metadata:
STATISTICS_MINIMUM=0
STATISTICS_MAXIMUM=52.424499511719
STATISTICS_MEAN=3.0666069784504
STATISTICS_MEDIAN=4.331070843283e-184
STATISTICS_MODE=-4
STATISTICS_STDDEV=3.567882216098
LAYER_TYPE=athematic

A couple of additional notes.  First, this is 2.5ft resolution data.  At this resolution, the differences in algorithms may be somewhat academic.  Also, of note– the range for the Zevenbergen-Thorn slope algorithm is greater than that for the local hydrologic slope.  This is an outcome I did not anticipate.

 

Update:

Hmmm, looks like I have more revisions and questions from the GDAL dev board. May have to revisit this… .  Forgot about the corners of my 3×3 window… .

Posted in Analysis, Ecology, GDAL | Tagged: , , , , | Leave a Comment »

Landscape Position: Conclusion? (part 2)

Posted by smathermather on December 7, 2011

From earlier post:

“I’ve managed to pilot most of a fast high resolution landscape position workflow with PovRay as my magic tool. The final steps I hope to pipe through PostGIS Raster. In the meantime a screenshot and description: blues are riparian, raw ocre, etc upland categories, grey is mostly flat lake plain and mid slopes, all derived from just a high res DEM input (no hydro lines as input, so it works on areas where you don’t know where the streams may be). There will be more categories in the final, in the meantime, welcome to December.”

I didn’t use PostGIS Raster, but the incredibly useful gdal_calc.py to finish out my landscape position calculations.  I will endeavor to post my code to github soon, but the parts and pieces include using gdal and PovRay.  PovRay helps with the sampling (really!) of nearby neighbors in the raster DEM, and gdal does the averaging and differencing of those to get relative landscape position.  I spent some time yesterday scratching my head over how to show all the new landscape position information on a readable and useable map, and after discussion with collegues, decided to use it to divide the world into two categories– riparian & lowland + headwaters & upland (not reflected yet in the labels).  To find out more about landscape position, follow this link, or better yet this one.  (BTW, the green is park land, blue is riparian/lowland/stream networks, purple is the basin boundary).

Riparian map based on landscape position, calculated using PovRay and GDAL.

Posted in Analysis, Ecology, Landscape Position, Other, POV-Ray | Tagged: , , , , , , , , , , , | Leave a Comment »

Landscape Position: Conclusion?

Posted by smathermather on November 30, 2011

I’ve managed to pilot most of a fast high resolution landscape position workflow with PovRay as my magic tool. The final steps I hope to pipe through PostGIS Raster. In the meantime a screenshot and description: blues are riparian, raw ocre, etc upland categories, grey is mostly flat lake plain and mid slopes, all derived from just a high res DEM input (no hydro lines as input, so it works on areas where you don’t know where the streams may be). There will be more categories in the final, in the meantime, welcome to December.

20111130-221953.jpg

Posted in Analysis, Ecology, Landscape Position, Other, POV-Ray | Tagged: , , , , , , , , , , , | 1 Comment »