Smathermather's Weblog

Remote Sensing, GIS, Ecology, and Oddball Techniques

Posts Tagged ‘Technical’

Using PostGIS for Hydrologic Modeling (reblog)

Posted by smathermather on June 17, 2016

The Problem We have to filter out the roads and ditches without removing streams that cross roads or follow them closely. I’m going to use PostGIS to find the intersection of the streams lines data with a buffered roads polygon. If the intersected line is less than 50% of the length of the stream line, […]

via Filtering Roads from Extracted Streams Data — GeoKota

Posted in Database, PostGIS, PostgreSQL, SQL | Tagged: , , , , , | Leave a Comment »

Leaflet Crosshairs

Posted by smathermather on February 18, 2016

Ok, I’m not the first to come up with that name, but I like it so it’s staying for now.

The problem space is this: Your web development team is good with forms. They build forms like diurnal animals wake — daily: day in, day out, day in, day out. They want to build a form for a mobile web app and it just happens to use HTML5 to grab the geolocation from field deployed phones, tablets, etc.. Great! Geo problem solved — when we collect our form data, we’ll get geo-data for free.

Not so fast, form-building Valentine — what’s the quality of those data? How accurate is that phone GPS? Is it good enough? The answer is probably yes, most of the time.

But, good enough most of the time isn’t good enough, Faye. What I want is an embedded map where I can move the crosshairs to the actual location inside the form. It’s like all the corrections you do back at the office, but you can do them in the field while the site is still fresh in your mind. Simple. Useful. (Also, infrastructure like Fulcrum App does it because it’s so simple and useful). Hence, questions like this: http://gis.stackexchange.com/questions/90225/how-to-add-a-floating-crosshairs-icon-above-leaflet-map and pages like this: https://www.mapbox.com/blog/help-search-MH370/

I couldn’t get the solution on Stack Exchange to work for me. Besides, I think its the wrong solution. I don’t want to move the icon back to center on the map moving, I want the map to move and the icon to stay stationary. It’s a fine (and probably irrelevant) distinction, but it feels important to me.

So, we build a small page that has the following features:

  • A crosshairs that is stationary
  • A map that moves
  • When the map moves, our lat/lon values update in our form

Main code as follows (careful — careless use of jquery follows):

<!DOCTYPE html>
<html>
<head>
	<title>Leaflet Crosshairs</title>

	<!-- Include meta tag to ensure proper rendering and touch zooming -->
	<!--<meta name="viewport" content="width=device-width, initial-scale=1" />-->
	<meta name="viewport" content="width=device-width, initial-scale=1.0, maximum-scale=1.0, user-scalable=no">

	<!-- Include the jQuery library -->
	<script src="//ajax.googleapis.com/ajax/libs/jquery/1.11.3/jquery.min.js"></script>
	<!-- Include jQuery Mobile stylesheets -->
	<link rel="stylesheet" href="//ajax.googleapis.com/ajax/libs/jquerymobile/1.4.5/jquery.mobile.min.css">
	<!-- Include the jQuery Mobile library -->
	<script src="//ajax.googleapis.com/ajax/libs/jquerymobile/1.4.5/jquery.mobile.min.js"></script>

	<!-- Include leaflet css and js -->
	<link rel="stylesheet" href="//cdn.leafletjs.com/leaflet/v0.7.7/leaflet.css" />
	<script src="//cdn.leafletjs.com/leaflet/v0.7.7/leaflet.js"></script>

	<style>
	body {
		padding: 0;
		margin: 0;
	}
        html, body, #map {
            height: 100%;
            width:100%;
        }

        #metamap {
            width: 100%;
            height: 300px;
        }
        #crosshair {
            position: relative;
            z-index: 10;
            height: 200px;
            vertical-align: middle;
        }
	#crosshair img {
		position: absolute;
		margin: 0;
		top: 50%;
		left: 50%;
		margin-right: -50%;
		transform: translate(-50%, -50%);
	}
	</style>
</head>
<body>

    <div id="metamap">
        <div id="map">
            <div id="crosshair"><img class="crosshair" src=crosshair.png /></div>
        </div>
    </div>
    <br />
    <hr />
    Latitude: <input type="text" id="txtLatitude" />
    <br /><br />
    Longitude: <input type="text" id="txtLongitude" />
    
        <script>
            // Initiate map
            var map = L.map('map');

            // load map
            L.tileLayer('//api.tiles.mapbox.com/v4/{id}/{z}/{x}/{y}.png?access_token=pk.eyJ1IjoiY2xldmVsYW5kLW1ldHJvcGFya3MiLCJhIjoiWHRKaDhuRSJ9.FGqNSOHwiCr2dmTH2JTMAA', {
                maxZoom: 20,
                id: 'mapbox.satellite'
            }).addTo(map);

            // Now a function to populate our form with latitude and longitude values
            function onMapMove(e) {
                // txtLatitude.val(map.getCenter());
                var locale = map.getCenter();
                $('#txtLatitude').val(locale.lat);
                $('#txtLongitude').val(locale.lng);
            }
            
            // Boilerplate...
            function onLocationError(e) {
                alert(e.message);
            }

            // When the map moves we run our function up above
            map.on('move', onMapMove);

            // Boilerplate
            map.on('locationerror', onLocationError);
            
            // When we load the map, we should zoom to our current position using device geolocation
            map.locate({ setView: true, maxZoom: 20 });
        </script>
</body>
</html>

<body>


<div id="metamap">

<div id="map">

<div id="crosshair"><img class="crosshair" src=crosshair.png /></div>

        </div>

    </div>

    

<hr />

    Latitude: <input type="text" id="txtLatitude" />
    

    Longitude: <input type="text" id="txtLongitude" />
    
        <script>
            // Initiate map
            var map = L.map('map');

            // load map
            L.tileLayer('//api.tiles.mapbox.com/v4/{id}/{z}/{x}/{y}.png?access_token=pk.eyJ1IjoiY2xldmVsYW5kLW1ldHJvcGFya3MiLCJhIjoiWHRKaDhuRSJ9.FGqNSOHwiCr2dmTH2JTMAA', {
                maxZoom: 20,
                id: 'mapbox.satellite'
            }).addTo(map);

            // Now a function to populate our form with latitude and longitude values
            function onMapMove(e) {
                // txtLatitude.val(map.getCenter());
                var locale = map.getCenter();
                $('#txtLatitude').val(locale.lat);
                $('#txtLongitude').val(locale.lng);
            }
            
            // Boilerplate...
            function onLocationError(e) {
                alert(e.message);
            }

            // When the map moves we run our function up above
            map.on('move', onMapMove);

            // Boilerplate
            map.on('locationerror', onLocationError);
            
            // When we load the map, we should zoom to our current position using device geolocation
            map.locate({ setView: true, maxZoom: 20 });
        </script>
</body>
</html>

Screen shot of app in action

Things to fix:

  • Alignment of crosshairs so they are properly centered
  • Better looking crosshairs
  • Rounding for those coordinate values
  • Do we need jQuery? Pro’ly not

That was my fun for the day. Shout-out to Tanetta Jordan my brilliant paired programmer for the day. Without her, this would have taken a week… .

Oh ya, and git-repo here: https://github.com/cleveland-metroparks/leaflet-crosshairs

Look there for updates that should include the improvements above.

Posted in Javascript, Leaflet | Tagged: , , , | Leave a Comment »

GDAL tools of the trade: GDAL, MrSid, and nearblack revisited

Posted by smathermather on February 11, 2016

internet_archiveBack in 2011 I wrote a post about GDAL, translating MrSID files and trimming off the artifacts so mosaics work. I’m revisiting some similar work this week and realized I had lost my copy of FWTools that can deal with MrSIDs. The Internet Archive came to the rescue.

At some point LizardTech licensed a DLL that allowed for compilation compatible with FWTools licensing so we could do things like use gdal_translate to covert a MrSID file to a tiff. Later on, the licensing changed, but binaries created prior to the change hold that original license. The WayBack Machine archived one such copy here: https://web.archive.org/web/20130127013649/http://home.gdal.org/fwtools/FWTools202.exe .

Download, install, and now for some Batch scripting fun… .

W:\Aircraft\naip_ortho>for %I in (*.sid) do gdal_translate %I %~nI.tif

W:\Aircraft\naip_ortho>gdal_translate ortho_1-1_1n_s_oh035_2015_1.sid ortho_1-1_
1n_s_oh035_2015_1.tif
Input file size is 58222, 49140
0…10…20.

Posted in GDAL, Image Processing | Tagged: , , , | 2 Comments »

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 »

What is the center line of a complex polygon? Routing Stream and Rivers Part Two

Posted by smathermather on August 10, 2012

The series

and

one additional example:

Posted in Analysis, Database, PostGIS, PostgreSQL, SQL | Tagged: , , , | Leave a Comment »

What is the center line of a complex polygon? Routing Stream and Rivers

Posted by smathermather on August 9, 2012

Continuing my posts on the centerline of a complex polygon (you can work your way backward from here), here’s it’s application to a riverine system (legend: light blue river, black centerline derived from routing through the skeleton of voronoi polygons from a densified stream bank set + the skipped skeleton bits in pink– read through the series if you if you don’t grok that explaination):

Posted in Analysis, Database, PostGIS, PostgreSQL, SQL | Tagged: , , , | 8 Comments »

Proper (ab)use of a database, contour interpolation using #postgresql #postgis #arcgis

Posted by smathermather on August 6, 2012

Anyone who has been following along at home knows I don’t think much like a DBA.  Sometimes that’s good; mostly it’s probably bad.  In this post, I hope it will be interesting.

The problem of the day is how to take engineering contours derived from breaklines, a lidar point cloud, and all the lot, and do a good job interpolating that to a DEM.  This is an interesting problem, as it’s another one of those problems with so many sub-optimal solutions.  In part, it’s an impossible problem– going from a surface model to contours is inherently a loss of information, so going the other way is difficult to do believably, both with respect to believable smoothness and the requisite feature retention.  This is a little like the pugilistic pie puppets problem: ‘flaky vs. tender’ of making good pie crust.  Doing one well seems to preclude doing the other well.

Now, Chris Hormann, of PovRay fame has a solution, there are solutions in ArcGIS, and in GRASS, and plenty of other attempts.  Here’s a decent list, which I haven’t checked for timeliness:  http://vterrain.org/Elevation/contour.html.  Since we have a ArcGIS license, it handles the dataset I have OK, the results checked out OK, this is what we chose to use for now.  I tried GRASS, but it was too slow for the amount of data being processed.  The command in ArcGIS ends up looking something like this:


TopoToRaster_sa "N2260635.shp elevation Contour;N2260640.shp elevation Contour;N2260645.shp elevation Contour;N2265633.shp elevation Contour;N2265640.shp elevation Contour;N2265645.shp elevation Contour;N2270635.shp elevation Contour;N2270640.shp elevation Contour;N2270645.shp elevation Contour;N2275650.shp elevation Contour" "C:\Temp\TopoToR_N2269.tif" 2 "2265031 645070 2270065 639945" 20 # # ENFORCE CONTOUR 40 # 1 0 # # # # # #

And thus we go from this:

To this:

Now to do it 623 more times… .  And in walks PostgreSQL for some abuse as a spatial engine with string manipulation to boot.

Each tile of contours I run through this approach needs the adjacent 8 tiles of contours loaded in order to avoid dreaded edge effects:

Except that sometimes, there aren’t 8 tiles available, like on the edge of large water bodies:

So, we can’t just use a heuristic that loads the 8 adjacent tiles– we need to enumerate these tiles.  An ST_Touches should do the trick here:

SELECT p.tile, a.tile AS tiles
FROM cuy_contour_tiles AS p, cuy_contour_tiles As a
WHERE ST_Touches(p.geom, a.geom);

But, that does not arrange our table very nicely:

I want something that looks like this:

N2225655, N2225660, N2220650 …

N2225660, N2225665, N220665 … etc.

Thank goodness for Postgres’ String_agg, and moreover for Leo Hsu and Regina Obe’s blog.

Adapting their entry, we can reorganize those entries like this:

SELECT
STRING_AGG(a.tile, ',' ORDER BY a.tile)
As tiles
FROM cuy_contour_tiles AS p, cuy_contour_tiles As a
WHERE ST_Touches(p.geom, a.geom)
GROUP BY p.tile
ORDER BY p.tile;

Ah, now that’s better:

And if we want, we can use an arbitrary delimiter, and build our entire command in one swell foop:

SELECT p.tile,
'TopoToRaster_sa "'
|| 'W:\contours\cuy_contours\shapefiles\' || p.tile || '.shp elevation Contour;'
|| 'W:\contours\cuy_contours\shapefiles\'
|| STRING_AGG(a.tile, '.shp elevation Contour;W:\contours\cuy_contours\shapefiles\' ORDER BY a.tile)
|| '.shp elevation Contour;" '
|| '"v:\svm\dem\cuyahoga\' || p.tile || '.tif" 2'
|| ' "' || replace(replace(replace(ST_Extent(p.geom)::text, 'BOX(', ''), ')', ''), ',', ' ') || '"'
|| ' 20 # # ENFORCE CONTOUR 40 # 1 0 # # # # '
As tiles
FROM cuy_contour_tiles AS p, cuy_contour_tiles As a
WHERE ST_Touches(p.geom, a.geom)
GROUP BY p.tile
ORDER BY p.tile

Ok.  I have 623 more tiles to do.  Until next time… .

Posted in Analysis, ArcGIS, Database, PostGIS, PostgreSQL, SQL | Tagged: , , , , , , , , , , | Leave a Comment »

What is the center line of a complex polygon? Routing as the solution

Posted by smathermather on August 2, 2012

I’ve had a series of posts (including this one) on finding the center line of a complex polygon– especially with an interest in finding the center line of streams and river bodies.  We have to give credit to my colleague Tom Kraft for solving the polygon centerline problem.  Tom struck upon the very elegant solution– derive the medial axis, and then use a few end points to route between to extract just the important features.  It’s not fully automated, but the amount of user input is so nominal as to be as close to automatic as I care about.  So far, we have this implemented in a mix of ArcGIS and PostGIS, but I can easily picture a pretty simple web interface with a PostGIS back end for extracting the centerlines of polygons in a nominally interactive way.  As a teaser, a non-stream example:

Posted in Analysis, Database, PostGIS, PostgreSQL, SQL | Tagged: , , , | 1 Comment »