Smathermather's Weblog

Remote Sensing, GIS, Ecology, and Oddball Techniques

Archive for the ‘UAS’ Category

Finding peace, finding ground: Drone flights for hydrologic modeling

Posted by smathermather on October 15, 2016

Š-L-M

Lately I’ve been helping with project to predict flooding in Dar es Salaam (“Home of Peace”), the capital former capital of Tanzania in East Africa. This is a really fun project for me, in part because most of the hard work is already done, but there still remain some tricky problems to solve.

This project is a part of Dar Ramani Huria, a “community-based mapping project… training university students and local community members to create highly accurate maps of the most flood-prone areas of the city using OpenStreetMap.”

Dar es Salaam (or Dar for short), is the fastest growing city in Africa, making the identification of flood prone zones during rainy seasons absolutely critical. The Ramani Huria crew has mapped the streams, ditches, and other waterways of Dar, as well as flown imagery (via drone) over much of the city critical parts of the city.

screen-shot-2016-10-15-at-3-01-04-pm

Problem Space:

The results are stunning, but using drone imagery and photogrammetric point clouds (instead of LiDAR) has it’s limitations.

One problem, that I’m not going to get into in any depth today, is the difficulty of conflating (vertically and horizontally) different datasets flown at different times with commodity GPS.

Another problem is the difficulty of turning photogrammetrically derived point clouds into Digital Terrain Models. There is proprietary software that does this well (e.g. LasTools and others), but we seek a free and open source alternative. Let’s visualize the problem. See below a photogrammetrically-derived point cloud (processed in Pix4D, visualized at https://plas.io):

screen-shot-2016-10-15-at-3-15-44-pm

We can directly turn this into a digital surface model (DSM):

screen-shot-2016-10-15-at-3-24-23-pm

5

We can see the underlying topography, but buildings and trees are much of the signal in the surface model. If we, for example, calculate height above the nearest stream as a proxy for flooding, we get a noisy result.

OSM to the Rescue

If you recall at the beginning of this post, I said this is “a really fun project for me” because “most of the hard work is already done”. Dar is a city with pretty much every building digitized. We can take advantage of this work to remove buildings from our point cloud. We’ll use the utilities associated with the Point Data Abstraction Layer (PDAL) library. Specifically, we’ll use PDAL’s ability to clip point cloud data with 2D OGR compatible geometries (a great tutorial on this here).

First thing we need to do is extract the buildings from OSM into shapefiles. For this, and easy way is to use Mapzen’s Metro Extracts:

screen-shot-2016-10-15-at-3-58-35-pm

I chose Datasets grouped into individual layers by OpenStreetMap tags (IMPOSM) as this gives me the capacity to just choose “Buildings”. I loaded those buildings into a PostGIS database for the next steps. What I seek is a shapefile of all the areas which are not buildings, essentially from this:

building

To this:

no_building

For this we’ll need to use ST_Difference + ST_ConvexHull. ST_ConvexHull will create a shape of the extent of our data, and we’ll use ST_Difference to cookie-cutter out all the building footprints:

DROP TABLE IF EXISTS "salaam-buff-diff-split1";
CREATE TABLE "salaam-buff-diff-split1" AS
WITH subset AS (
	SELECT ST_Transform(geom, 32737) AS geom FROM "salaam-buildings" LIMIT 20000
	),
extent AS (
	SELECT ST_ConvexHull(ST_Union(geom)) AS geom FROM subset
	),
unioned AS (
	SELECT ST_Union(geom) AS geom FROM subset
	),
differenced AS (
	SELECT ST_Difference(a.geom, b.geom) AS geom FROM
	extent a, unioned b
	)
SELECT 1 AS id, geom FROM differenced;

This gives us a reasonable result. But, when we use this in PDAL, (I assume) we want a subdivided shape to ensure we don’t have to access the entire point cloud at any given time to do our clipping. We’ll add to this process ST_SubDivide, which will subdivide our shape into parts not to exceed a certain number of nodes. In this case we’ll choose 500 nodes per shape:

DROP TABLE IF EXISTS "salaam-buff-diff-split1";
CREATE TABLE "salaam-buff-diff-split1" AS
WITH subset AS (
	SELECT ST_Transform(geom, 32737) AS geom FROM "salaam-buildings" LIMIT 20000
	),
extent AS (
	SELECT ST_ConvexHull(ST_Union(geom)) AS geom FROM subset
	),
unioned AS (
	SELECT ST_Union(geom) AS geom FROM subset
	),
differenced AS (
	SELECT ST_Difference(a.geom, b.geom) AS geom FROM
	extent a, unioned b
	)
SELECT 1 AS id, ST_Subdivide(geom, 500) AS geom FROM
differenced;

subdivided

Finally, if we want to be sure to remove the points from the edges of buildings (we can assume the digitized buildings won’t perfectly match our point clouds), then we should buffer our shapes:

DROP TABLE IF EXISTS "salaam-buff-diff-split1";
CREATE TABLE "salaam-buff-diff-split1" AS
WITH subset AS (
	SELECT ST_Transform(geom, 32737) AS geom FROM "salaam-buildings" LIMIT 20000
	),
buffered AS (
	SELECT ST_Buffer(geom, 2, 'join=mitre mitre_limit=5.0') AS geom FROM subset
	),
extent AS (
	SELECT ST_ConvexHull(ST_Union(geom)) AS geom FROM subset
	),
unioned AS (
	SELECT ST_Union(geom) AS geom FROM Buffered
	),
differenced AS (
	SELECT ST_Difference(a.geom, b.geom) AS geom FROM
	extent a, unioned b
	)
SELECT 1 AS id, ST_Subdivide(geom, 500) AS geom FROM
differenced;

If we recall our before point cloud:

screen-shot-2016-10-15-at-3-15-44-pm

screen-shot-2016-10-15-at-4-38-53-pm

Now we have filtered out most buildings:

screen-shot-2016-10-15-at-4-36-15-pm

screen-shot-2016-10-15-at-4-36-31-pm

In order to do this filtering in PDAL, we do two things. First we create a json file that defines the filter:

{
  "pipeline":[
    "/data/2015-05-20_tandale_merged_densified_point_cloud_part_1.las",
    {
      "type":"filters.attribute",
      "dimension":"Classification",
      "datasource":"/data/salaam-buff-diff-split.shp",
      "layer":"salaam-buff-diff-split",
      "column":"id"
    },
    {
      "type":"filters.range",
      "limits":"Classification[1:1]"
    },
    "/data/2015-05-20_tandale_merged_densified_point_cloud_part_1_nobuild_buff.las"
  ]
}

Then we use these json definitions to apply the filter:

nohup sudo docker run -v /home/gisuser/docker/:/data pdal/pdal:1.3 pdal pipeline /data/shape-clip.json

Problems continue

This looks pretty good, but as we interrogate the data, we can see the artifacts of trees and some buildings still linger in the dataset.

screen-shot-2016-10-15-at-4-37-38-pm

screen-shot-2016-10-15-at-4-38-03-pm

 

Enter the Progressive Morphological Filter

It is possible to use the shape of the surface model to filter out buildings and trees. To do so, we start with the assumption that the buildings and vegetation shapes distinctive from the shape of the underlying ground. We have already used the hard work of the OpenStreetMap community to filter most of the buildings, but we still have some buildings and plenty of trees. PDAL has another great tutorial for applying this filter which we’ll leverage.

Again we need a JSON file to define our filter:

{
  "pipeline": {
    "name": "Progressive Morphological Filter with Outlier Removal",
    "version": 1.0,
    "filters": [{
        "name": "StatisticalOutlierRemoval",
        "setMeanK": 8,
        "setStddevMulThresh": 3.0
      }, {
        "name": "ProgressiveMorphologicalFilter",
        "setCellSize": 1.5
    }]
  }
}

And then use that filter to remove all the morphology and statistical outliers we don’t want:

nohup sudo docker run -v /home/gisuser/docker/:/data pdal/pdal:1.3 pdal pcl -i /data/2015-05-20_tandale_merged_densified_point_cloud_part_1_nobuild_buff.las -o /data/nobuild_filtered2.las -p /data/sor-pmf.json

This command will remove our colorization for the points, so we’ll see the colorization according to height only:

screen-shot-2016-10-15-at-4-53-55-pm

screen-shot-2016-10-15-at-4-54-49-pm

 

What remains is to interpolate this into digital terrain model. That is a project for another day.

Posted in 3D, Drone, Other, PDAL, Photogrammetry, UAS | Tagged: , , , | 1 Comment »

OpenDroneMap — Paris Code Sprint

Posted by smathermather on February 29, 2016

I failed to make it to the Paris Code Sprint. It just wasn’t in the cards. But, my colleague Dakota and I sprinted anyway, with some help and feedback from the OpenDroneMap community.

So, what did we do? Dakota did most of the work. He hacked away at the cmake branch of ODM, a branch set up by Edgar Riba to substantially improve the installation process for ODM.

  • Fixed odm_orthophoto in the branch so that it produces geotiffs
  • Fixed PMVS so that it is multithreaded again
  • Added rerun-all and rerun-from function
  • Integrated @lupas78’s additions for an xyz point cloud output
  • Added benchmarking which is an important soft number for when we have code changes
  • (Technically before the sprint) wrote the first test for OpenDroneMap
  • Cleaned code
What did I do? Mostly, I got caught up with the project. I haven’t been very hands on since the python port, let alone the cmake branch, so I became a little more pythonistic by just trying to successfully modify the code.
  • I also added PDAL to the build processs
  • And I inserted PDAL into the point cloud translation process.

Currently, this means we’ve dropped support for LAZ output, as I haven’t successfully built PDAL with LAZ support, but it stages the project for LAZ support through PDAL, and allows us to tap into additional PDAL functionality in the future.

It was an intensive couple of days that would have been improved with French wine, but we were in Parma (Ohio). So, a shout out to the coders in Paris at the same time, and cheers to all.

Posted in 3D, Drone, OpenDroneMap, OpenDroneMap, Photogrammetry, PMVS, UAS | Tagged: , | Leave a Comment »

Moar kite flight pics

Posted by smathermather on April 27, 2015

        

Posted in 3D, Bundler, Image Processing, OpenDroneMap, OpenDroneMap, Optics, Photogrammetry, PMVS, UAS | Tagged: , , , , , , , , , | Leave a Comment »

Kite flight (too windy for balloons, ahem “aerostats”)

Posted by smathermather on April 20, 2015

Inflation of aerostat

The aerostat hanger.

The end of the string.

The 9-footer is just so stable. But not enough wind to lift the cameras this day.

And so we send up the 16-foot workhorse. See that little dot? That’s the camera array.

The 16-footer flew nice and vertical, but pulled really hard. Processed images to follow soon.

Canon S100s from Kaptery — the silver one is an NIR adapted one; the black one is RGB color.


Edit: forgot the camera array:

CIR image from balloon:

IR image from the flight.

IR image from the flight.

Posted in 3D, Bundler, Image Processing, OpenDroneMap, OpenDroneMap, Photogrammetry, PMVS, UAS | Tagged: , , , , , , , , , | Leave a Comment »

Airspace — A deep rabbit hole

Posted by smathermather on October 25, 2014

In previous maps we looked at Class B, C, and D airspace. Let’s add in Class E0 and E5… (not yet in 3D):

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

Map showing Class B, C, D, E0, and E5 airspace

Map showing Class B, C, D, E0, and E5 airspace

Previous posts:

https://smathermather.wordpress.com/2014/10/25/airspace-is-complicated-and-so-i-abuse-postgis-once-again/

and

https://smathermather.wordpress.com/2014/10/25/airspace-is-complicated-and-so-i-abuse-postgis-once-again-reprise/

Posted in 3D, Cartography, Database, Drone, PostGIS, PostgreSQL, SQL, UAS | Tagged: , , , , , | Leave a Comment »

Airspace is complicated — and so I abuse PostGIS once again — Reprise…

Posted by smathermather on October 25, 2014

In the previous post: https://smathermather.wordpress.com/2014/10/25/airspace-is-complicated-and-so-i-abuse-postgis-once-again/ we explore the 3D shape and complexity of controlled airspace.

Now here’s the rest of the code. We’ll add our affine transformation ala Seth Fitsimmons:

    SELECT 
        ST_Affine(
            ST_Rotate(geom, -pi() / 2),
            -- isometric
            cos(pi() / 6), -cos(pi() / 6), 0,
            sin(pi() / 6), sin(pi() / 6), 1,
            0, 0, 0,
            0, 0, 0
        )
    AS geom

And integrate that into our original function:

-- Inputs are a geometry and an elevation to move the geometry to
CREATE OR REPLACE FUNCTION threed_iso(footprint geometry, elevation numeric)
  RETURNS geometry AS
$BODY$

-- Force 3D, then translate to the input elevation
WITH floor AS
(
    SELECT ST_Translate( ST_Force3DZ(footprint), 0, 0, elevation ) AS geom
),
-- Now make isometric (and begin Seth Code)
iso AS
(
    SELECT 
        ST_Affine(
            ST_Rotate(geom, -pi() / 2),
            -- isometric
            cos(pi() / 6), -cos(pi() / 6), 0,
            sin(pi() / 6), sin(pi() / 6), 1,
            0, 0, 0,
            0, 0, 0
        )

    AS geom
    FROM floor
)
-- We'll force it back to 3D so QGIS is happy
SELECT ST_Force2D(geom) FROM iso
;
$BODY$
  LANGUAGE sql VOLATILE
  COST 100;

And voila!

DROP TABLE IF EXISTS class_c_isoc;

CREATE TABLE class_c_isoc AS
	SELECT gid, airspace, name, lowalt, highalt, threed_iso(geom, lowalt::numeric * 5) AS geom
	FROM class_c_subset;

Let’s take a look at Washington, DC and surrounds, another nice complicated example:

3D Figure of DC controlled airspace

3D Figure of DC controlled airspace

And again with map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL:

3D Figure of DC controlled airspace

3D Figure of DC controlled airspace

 

Posted in 3D, Cartography, Database, Drone, PostGIS, PostgreSQL, SQL, UAS | Tagged: , , , , , | 3 Comments »

Airspace is complicated — and so I abuse PostGIS once again

Posted by smathermather on October 25, 2014

Let’s ignore for a moment the drone hobbiest / enthusiast. What is the shape of airspace for airplanes and commercial and government unmanned aircraft flying under Certificates of Authorization, and how can we visualize it?

Thanks to Anita in the last post, we have the Class B,C,D,E Airspace Shape Files which helps us define the overall shape of controlled airspace.

Map of Detroit, Cleveland, and Pittsburgh Class B Airspace

Map of Detroit, Cleveland, and Pittsburgh Class B Airspace

But, these are 3D things. I want to visualize them thus. Let us put some constraints on the problem. Let’s do it all in PostGIS, that way we can see it in QGIS or on the web. Let’s not use full PostGIS 3D (i.e. the SFCGAL library), not because it isn’t awesome (it truly is) but because it can be hard to install at the moment (although see https://github.com/vpicavet/docker-pggis for an easy way with docker). Finally, true 3D with unconstrained viewing angles and 100% flexibility is… is… well it usually sucks. So, we’ll stick to isometric viewing (thanks to Seth Fitzsimmons of Stamen http://stamen.com/ for his PostGIS isometric code which will be released upon his word). (Update — all the code is there…):

-- Inputs are a geometry and an elevation to move the geometry to
CREATE OR REPLACE FUNCTION threed_iso(footprint geometry, elevation numeric)
  RETURNS geometry AS
$BODY$

-- Force 3D, then translate to the input elevation
WITH floor AS
(
    SELECT ST_Translate( ST_Force3DZ(footprint), 0, 0, elevation ) AS geom
),
-- Now make isometric (and begin Seth Code)
iso AS
(
    SELECT
        ST_Affine(
            ST_Rotate(geom, -pi() / 2),
            -- isometric
            cos(pi() / 6), -cos(pi() / 6), 0,
            sin(pi() / 6), sin(pi() / 6), 1,
            0, 0, 0,
            0, 0, 0
        )

    AS geom
    FROM floor
)
-- We'll force it back to 3D so QGIS is happy
SELECT ST_Force2D(geom) FROM iso
;
$BODY$
  LANGUAGE sql VOLATILE
  COST 100;

Ok, now let’s rock some geometries with this bad function:

DROP TABLE IF EXISTS class_c_isoc;

CREATE TABLE class_c_isoc AS
	SELECT gid, airspace, name, lowalt, highalt, threed_iso(geom, lowalt::numeric * 5) AS geom
	FROM class_c_subset;

And what do our controlled airspaces look like?

Isometric view of Cleveland controlled airspace

Isometric view of Cleveland controlled airspace

Kind of conical, in this case with some “wings” that serve as approaches. It makes sense, I guess. At the center, where the airport is, controlled airspace goes from the ground up. As we get further from the airport, a set of concentric rings starting at higher elevations provide a controlled space that allows for landing and taking off.

Image of Cleveland airspace with Stamen Toner map as backdrop

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

There are more complicated ones, naturally. We need look no further than Detroit for additional complexity:

Visualization of Detroit controlled airspace

Visualization of Detroit controlled airspace

airspace_detroit_toner

No wonder every time I fly through Detroit, no matter where I’m going, I end up flying over Lake Erie.

If we want really complicated, we need look no further than Cincinnati:
airspace_cinci

What is going on there? Is that due to shape of the hills, or city boundaries?

Finally, what does airspace look like over Ohio, West Virginia, and Pennsylvania (etc.), overall?
airspace_is_complicated_regional

And while the following map isn’t quite right, here is a figure including many of the small airports sans controlled airspace:

View of all controlled and uncontrolled airspace across Ohio and neighbors.

View of all controlled and uncontrolled airspace across Ohio and neighbors.

May an aeronautical pox be upon you!
The above line was not intended in bad taste, but just an homage to the “red dot fever” of early neo-geography (neo-geography which I’m informed definitionally doesn’t exist). Only a few minutes ago, it dawned on me that the deleted phrase could be misinterpreted these days… .

On a related note, if you want an interesting analysis of ebola, read Zeynep Tufekci’s analysis.

(for the record, all heights are exagerated by 5x, for clarity of reading).

Also, in case it wasn’t obvious: Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL.

Posted in 3D, Cartography, Database, Drone, PostGIS, PostgreSQL, SQL, UAS | Tagged: , , , , , | 4 Comments »

Someone is wrong on the internet…

Posted by smathermather on October 23, 2014

Ok @Mapbox, I’ve given you enough grace. I pulled the grumpy-old-man thing with Kenneth Field from ESRI a couple of weeks ago, and ended up apologizing. If i only give ESRI a few hours of grace, 3 months of grace for you is more than fair.

Now, to raise my foot up. We’ll see if it ends up in my mouth again.

Ok, so Mapbox and drones. Mapbox and drones. Those folks are excited about drones. They’ve got a toolchain their working on, and integration into their existing services, and lots of posts about drones. Downright giddiness I say:

https://duckduckgo.com/?q=site%3Amapbox.com+drones&t=canonical

And then this:

Really? 3 months ago, you launch “Don’t Fly Here”, create a repo for feedback and improvements to the data, but the biggest problem that you haven’t even tried to fix is the map really under represents the restrictions on where you should fly if you are trying to stay out of controlled air space, and the major update you do is adding temporary flight restrictions. Glitz and glory over getting the fundamentals right.

Here’s the current Mapbox no fly zones:

not_quite_right

But, we are missing a lot of smaller airports that also have controlled airspace. Let’s be considerate to our hobbiest drone friends and put them in with 3-mile buffers and refine our 5-mile buffers on medium and large airports to include the edges of the airport, just to be on the safe side:

a_little_righter

Well, that’s a bit less inviting. But hey, still plenty of places to fly, right? Uh, no. We forgot, according to the Aeronautical Information Manual, in certain busy airspace, there is a Mode C requirement for all craft flying in the airspace. This means you need to have a transponder to fly in this space between ground and 10,000ft above mean sea level. This is a transponder that weighs several times what your drone can lift — in otherwords a “non-starter”, effectively making these zones no-fly zones as well. What does Hopkins look like with it’s 30 nautical mile Mode C requirement?:

even_righter

Now, to be fair, this is murky legal territory at best and I am not a lawyer. How much applies to hobbiests given the 1981 Advisory Circular I can’t say I know. But, if we are to propose a map to clarify where we can and cannot fly as hobbiests, then we should be including as much information (in a simple and easy to use way) as we can. On these grounds, “Don’t Fly Here” fails.

It’s cool though. It’s an open source data project in an open source community. Community contributions and knowledge will fix all mistakes in time, so my pull request to get to the second map will be reviewed and rejected on solid grounds or integrated, right?

https://github.com/mapbox/drone-feedback/pull/40

Three months later, I’m still waiting… . Fix your map please. Engage your repo users. Do this right. Please.

Posted in Drone, UAS | Tagged: , , | 11 Comments »

Announcing OpenDroneMap — Software for civilian (and humanitarian?) UAS post processing

Posted by smathermather on September 15, 2014

OpenDroneMap logo

This past Friday at FOSS4G in Portland, I announced the (early) release of OpenDroneMap, a software toolchain for civilian (and humanitarian?) UAS/UAV image processing. The software is currently a simple fork of https://github.com/qwesda/BundlerTools, and will process from unreferenced overlapping photos to an unreferenced point cloud. Directions are included in the repo to create a mesh and UV textured mesh as the subsequent steps, but the aim is to have this all automated in a single work flow.


Projects like Google Streetview, Mapillary, PhotoSynth, and most small UAS (drone) postprocessing software, such as that offered by senseFly, share a commonality– they all use computer vision techniques to create spatial data from un-referenced photography.

Screenshot of drone image thumbnails

OpenDroneMap is an open source project to unlock and make easy-to-use related related computer vision techniques, so that whether from street level photos or from civilian drone imagery, the average user will be able to generate point clouds, 3D surfaces models, DEMs, and orthophotography from processing unreferenced photos.

odm_flow

 


 

Screen shot of textured mesh as viewed in MeshLab


To those who may be wondering — wow cool, but what happens to the data at the end of the day? How do we share it back to a common community? The aim is for the toolchain to also be able to optionally push to a variety of online data repositories, pushing hi-resolution aerials to OpenAerialMap, point clouds to OpenTopography, and pushing digital elevation models to an emerging global repository (yet to be named…). That leaves only digital surface model meshes and UV textured meshes with no global repository home. (If anyone is working on global storage of geographically referenced meshes and textured meshes, please get in touch…).


 

So, try it out: http://opendronemap.org will point you to the repo. Clone it, fork it, try it out. Let me know what you think.

Test data can be found here: https://github.com/OpenDroneMap/odm_data (credit Fred Judson, Ohio Department of Transportation)

Presentations on it can be found here: https://github.com/OpenDroneMap/presentations and eventually here:


 

PostScript: Re: meshes and pointclouds on the web, Howard Butler and others are working on some pretty cool tools for handling just this problem technologically. Check out, for example http://plas.io/ and https://github.com/hobu/greyhound

Posted in 3D, Bundler, Camera Calibration, Drone, Image Processing, OpenDroneMap, Optics, Photogrammetry, PMVS, UAS | Tagged: , , , , , , | 3 Comments »

Texturing surfaces

Posted by smathermather on September 8, 2014

Evaluating the application of texture to 3D models as we approach my presentation on OpenDroneMap. Promising stuff:

image

Posted in Drone, OpenDroneMap, UAS | Tagged: | Leave a Comment »