Smathermather's Weblog

Remote Sensing, GIS, Ecology, and Oddball Techniques

Archive for February, 2016

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 »

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 »

Valentine’s map for my daughter

Posted by smathermather on February 12, 2016

It’s fun every now and then to play with Stamen Design’s Map Stack. Tonight as I was playing with it in part to wind down from a long work week, my 3-year-old requested a pink and purple map.

So, a pink and purple map. I can do that:

Screen Shot 2016-02-12 at 10.18.44 PM

The cool thing about Map Stack, is I can share a link for my creation (click on the link above image), but I can also now embed this in a QGIS map.

http://mapstack.stamen.com/edit.html#watercolor%5Bsat=0%5D;color:$ff3030%5Bmask=!mapbox-water,comp=multiply,alpha=50%5D;toner-hybrid;color:$8629a6%5Bmask=!mapbox-water,comp=multiply,alpha=50%5D;color:$60c9fe%5Bmask=mapbox-water,comp=multiply,alpha=60%5D/2/30.1/-48.3

I can embed the link to that map in a specially crafted XML file that can be loaded in QGIS or many projects that use the GDAL library:

<GDAL_WMS>
    <Service name="TMS">
        <ServerUrl>http://a.sm.mapstack.stamen.com/((watercolor,$fff[hsl-saturation]),(mapbox-water,$ff3030[source-out@p])[multiply@50],toner-hybrid,(mapbox-water,$8629a6[source-out@p])[multiply@50],(mapbox-water,$60c9fe[source-in])[multiply@60])/${z}/${x}/${y}.png</ServerUrl>
    </Service>
    <DataWindow>
        <UpperLeftX>-20037508.34</UpperLeftX>
        <UpperLeftY>20037508.34</UpperLeftY>
        <LowerRightX>20037508.34</LowerRightX>
        <LowerRightY>-20037508.34</LowerRightY>
        <TileLevel>18</TileLevel>
        <TileCountX>1</TileCountX>
        <TileCountY>1</TileCountY>
        <YOrigin>top</YOrigin>
    </DataWindow>
    <Projection>EPSG:3857</Projection>
    <BlockSizeX>256</BlockSizeX>
    <BlockSizeY>256</BlockSizeY>
    <BandsCount>3</BandsCount>
    <Cache />
</GDAL_WMS>
Screen Shot 2016-02-12 at 10.34.16 PM

View of map in QGIS


Screen Shot 2016-02-12 at 10.33.37 PM

View of map in QGIS reprojected to World Robison

Posted in Other | 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 »

PostgreSQL table drop — a little fun with ouroboros

Posted by smathermather on February 11, 2016

My latests posts on PDAL have been fun. For the moment, a quick bit of code for dropping all the tables in your PostgreSQL database. BTW, the following is a bad idea. Many bad ideas are really useful. This is a really useful but bad idea:

echo "\dt" | psql | grep layer | awk '{print "DROP TABLE \"" $3"\";"}' | psql

How does this work? Let’s do a quick step through.

First we echo “\dt”.

$ echo "\dt"
\dt

This just prints literally (the literally kind of literally) “\dt”. \dt is a psql command for listing all the tables in your database. We “pipe” (this symbol is a pipe: |) that into psql. Essentially what this does is run this command within psql as follows:

$ echo "\dt" | psql
               List of relations
 Schema |        Name        | Type  |  Owner  
--------+--------------------+-------+---------
 public | layer_0-1.5        | table | user
 public | layer_1.5-3        | table | user
 public | layer_105-150      | table | user
 public | layer_15-30        | table | user
 public | layer_150-200      | table | user
 public | layer_3-6          | table | user
 public | layer_30-45        | table | user
 public | layer_45-60        | table | user
 public | layer_6-15         | table | user
 public | layer_60-105       | table | user
 public | paw_points         | table | user
 public | pointcloud_formats | table | user
 public | spatial_ref_sys    | table | user
(13 rows)

Now that we have a list of tables, we can manipulate that list in order to manipulate our tables. We’ll use the grep command to go line by line and only return the ones with the word “layer” in them. These are the table names that we are interested in dropping.

$ echo "\dt" | psql | grep layer
 public | layer_0-1.5        | table | user
 public | layer_1.5-3        | table | user
 public | layer_105-150      | table | user
 public | layer_15-30        | table | user
 public | layer_150-200      | table | user
 public | layer_3-6          | table | user
 public | layer_30-45        | table | user
 public | layer_45-60        | table | user
 public | layer_6-15         | table | user
 public | layer_60-105       | table | user

Now we use the awk command to just grab the column we want. By default, awk assumes spaces are our column delimiter, so we’ll grab the third column, however we could also use the pipes as our delimiter with the -F flag:

$ echo "\dt" | psql | grep layer | awk '{print $3}'
layer_0-1.5
layer_1.5-3
layer_105-150
layer_15-30
layer_150-200
layer_3-6
layer_30-45
layer_45-60
layer_6-15
layer_60-105

Let us next extend the awk command to print the commands we want to run against our database. Note I had to escape my double quotes with a “\”:

$ echo "\dt" | psql | grep layer | awk '{print "DROP TABLE \"" $3"\";"}' 
DROP TABLE "layer_0-1.5";
DROP TABLE "layer_1.5-3";
DROP TABLE "layer_105-150";
DROP TABLE "layer_15-30";
DROP TABLE "layer_150-200";
DROP TABLE "layer_3-6";
DROP TABLE "layer_30-45";
DROP TABLE "layer_45-60";
DROP TABLE "layer_6-15";
DROP TABLE "layer_60-105";

Now we have the commands we need to feed back into psql to run against the database. Here is where the dragon eats its own tail:

echo "\dt" | psql | grep layer | awk '{print "DROP TABLE \"" $3"\";"}' | psql
DROP TABLE
DROP TABLE
DROP TABLE
DROP TABLE
DROP TABLE
DROP TABLE
DROP TABLE
DROP TABLE
DROP TABLE
DROP TABLE

Posted in Database, PostGIS, PostgreSQL | Tagged: , | 3 Comments »

Taking Slices from LiDAR data: Part V

Posted by smathermather on February 10, 2016

For this post, let’s combine the work in the last 4 posts in order to get a single pipeline for doing the following:

  • Calculate relative height of LiDAR data
  • Slice that data into bands of heights
  • Load the data into a PostgreSQL/PostGIS/pgPointCloud database.
#!/bin/bash 

# readlink gets us the full path to the file. This is necessary for docker
readlinker=`readlink -f $1`
# returns just the directory name
pathname=`dirname $readlinker`
# basename will strip off the directory name and the extension
name=`basename $1 .las`

# Docker run allows us to leverage a pdal machine with pcl built in,
# thus allowing us to calculate height.
# See http://www.pdal.io/tutorial/calculating-normalized-heights.html
docker run -v $pathname:/data pdal/master pdal translate //data/"$name".las //data/"$name"_height.bpf height --writers.bpf.output_dims="X,Y,Z,Intensity,ReturnNumber,NumberOfReturns,ScanDirectionFlag,EdgeOfFlightLine,Classification,ScanAngleRank,UserData,PointSourceId,Height";

# Now we split the lidar data into slices of heights, from 0-1.5 ft, etc.
# on up to 200 feet. We're working in the Midwest, so we don't anticipate
# trees much taller than ~190 feet
for START in 0:1.5 1.5:3 3:6 6:15 15:30 30:45 45:60 60:105 105:150 150:200
 do
  # We'll use the height classes to name our output files and tablename.
  # A little cleanup is necessary, so we're removing the colon ":".
  nameend=`echo $START | sed s/:/-/g`

  # Name our output
  bpfname=$name"_"$nameend.bpf

  # Implement the height range filter
  pdal translate $name"_height".bpf $bpfname -f range --filters.range.limits="Height[$START)"

  # Now we put our data in the PostgreSQL database.
  pdal pipeline -i pipeline.xml --writers.pgpointcloud.table='layer_'$nameend --readers.bpf.filename=$bpfname --writers.pgpointcloud.overwrite='false'
done

Now, we can use parallel to make this run a little faster:

find . -name "*.las" | parallel -j6 ./pdal_processor.sh {}&

Sadly, we can run into issues in running this in parallel:

PDAL: ERROR:  duplicate key value violates unique constraint "pointcloud_formats_pkey"
DETAIL:  Key (pcid)=(1) already exists.


PDAL: ERROR:  duplicate key value violates unique constraint "pointcloud_formats_pkey"
DETAIL:  Key (pcid)=(1) already exists.

This issue is a one time issue, however — we just can’t parallelize table creation. Once the tables are created however, I believe we can parallelize without issue. I’ll report if I find otherwise.

Posted in 3D, Database, Docker, LiDAR, Other, PDAL, pointcloud, PostGIS, PostgreSQL | Tagged: , , , , , | Leave a Comment »

Taking Slices from LiDAR data: Part IV

Posted by smathermather on February 5, 2016

Trans-Alaska Pipeline System Luca Galuzzi 2005

In PDAL, a pipeline file can be used to do a variety of operations. Within the following context, I think of a pipeline file like an ad hoc preferences file, where I can use an external command to iterate through the things I want to change, while holding constant everything else in the pipeline file.

In my use case for this vignette, I’ll use the pipeline file to hold my database preferences for putting the point clouds into my PostGIS database. For the record, I’m using vpicavet’s docker pggis as the starting place for installing PostGIS with the pgpointcloud extension. I have adapted the following pipeline file from the PDAL writers.pgpointcloud example.

<?xml version="1.0" encoding="utf-8"?>
<Pipeline version="1.0">
  <Writer type="writers.pgpointcloud">
    <Option name="connection">
      host='localhost' dbname=‘lidar’ user='user' password=‘password’
    </Option>
    <Option name="table">54001640PAN_heightasz_0-1.5</Option>
    <Option name="compression">dimensional</Option>
<!--    <Option name="srid">3270</Option> -->
    <Filter type="filters.chipper">
      <Option name="capacity">400</Option>
      <Reader type="readers.las">
          <Option name="filename">54001640PAN_heightasz_0-1.5.las</Option>
<!--          <Option name="spatialreference">EPSG:3270</Option> -->
      </Reader>
    </Filter>
  </Writer>
</Pipeline>

Some things to note. I have commented out the SRID and readers.las.spatialreference in the XML above. We’ll rely on PDAL to discover this, and use the default output of epsg:4326 to store our point clouds for the time being.

Our wrapper script for the pipeline file is very simple. We will use the wrapper script to specify the table and the input file name.

#!/bin/bash

TABLE=`basename $1 .bpf`
INPUT=$1

#echo $TABLE
pdal pipeline -i pipeline.xml --writers.pgpointcloud.table=$TABLE --readers.bpf.filename=$INPUT

Now to use, we’ll use GNU Parallel, as much for it’s XArgs like functionality as scalability:

ls *.bpf | parallel -j6 ./pipeliner.sh {}

Now we can see what tables got loaded:

psql -d lidar
psql (9.5.0)
Type "help" for help.

lidar=# \dt
                    List of relations
 Schema |             Name              | Type  |  Owner  
--------+-------------------------------+-------+---------
 public | 54001640PAN_heightasz_0-1.5   | table | user
 public | 54001640PAN_heightasz_1.5-6   | table | user
 public | 54001640PAN_heightasz_100-200 | table | user
 public | 54001640PAN_heightasz_30-45   | table | user
 public | 54001640PAN_heightasz_45-55   | table | user
 public | 54001640PAN_heightasz_55-65   | table | user
 public | 54001640PAN_heightasz_6-30    | table | user
 public | 54001640PAN_heightasz_65-75   | table | user
 public | 54001640PAN_heightasz_75-85   | table | user
 public | 54001640PAN_heightasz_85-100  | table | user
 public | pointcloud_formats            | table | user
 public | spatial_ref_sys               | table | user
(12 rows)

W00t! We’ve got point clouds in our database! Next, we will visualize the data, and extract some analyses from it.

Posted in 3D, LiDAR, Other, PDAL, pointcloud | Tagged: , , | Leave a Comment »

Taking Slices from LiDAR data: Part III

Posted by smathermather on February 5, 2016

forest_structure
Borrowed from: http://irwantoshut.com

Continuing my series on slicing LiDAR data in order to analyze a forest, one of the objectives of the current project is to understand the habitats that particular species of birds prefer. This will be accomplished using field info from breeding bird surveys combined with LiDAR data of forest structure to help predict what habitats are necessary for particular species of breeding birds.

There are a number of studies doing just this for a variety of regions which I will detail later, but suffice it to say, structure of vegetation matters a lot to birds, so using LiDAR to map out structure can be an important predictive tool for mapping bird habitat. Being able to do that at scale across entire ecoregions—well, that’s just an exciting prospect.
Let’s get to some coding. I would like to take a few slices through our forest based on functional height groups:

Forest Canopy >15 meters
Sub-Canopy 10-15 meters
Tall Shrub 2-10 meters
Short Shrub 0.5-2 meters

(For the record, we don’t have to get these perfectly right. Sub-canopy, for example could be taller or shorter depending on how tall a forest it is in. This is just a convenience for dividing and reviewing the data for the time being. Also, we’ll split a little finer than the above numbers just for giggles.)

We’ll start by just echoing our pdal translate to make sure we like what we are getting for output:

for START in 0:0.5 0.5:1 1:2 2:5 5:10 10:15 15:20 20:35 35:50 50:75 75:100 100:200
 do
  nameend=`echo $START | sed s/:/-/g`
  echo pdal translate 54001640PAN_heightasz.bpf 54001640PAN_heightasz_$nameend.bpf -f range --filters.range.limits="Height[$START)"
done

Thus we get the following output:

pdal translate 54001640PAN_heightasz.bpf 54001640PAN_heightasz_0-0.5.bpf -f range --filters.range.limits=Height[0:0.5)
pdal translate 54001640PAN_heightasz.bpf 54001640PAN_heightasz_0.5-1.bpf -f range --filters.range.limits=Height[0.5:1)
pdal translate 54001640PAN_heightasz.bpf 54001640PAN_heightasz_1-2.bpf -f range --filters.range.limits=Height[1:2)
pdal translate 54001640PAN_heightasz.bpf 54001640PAN_heightasz_2-5.bpf -f range --filters.range.limits=Height[2:5)
... etc.

Let’s remove our echo statement so this actually runs:

for START in 0:0.5 0.5:1 1:2 2:5 5:10 10:15 15:20 20:35 35:50 50:75 75:100 100:200
 do
  nameend=`echo $START | sed s/:/-/g`
  pdal translate 54001640PAN_heightasz.bpf 54001640PAN_heightasz_$nameend.bpf -f range --filters.range.limits="Height[$START)"
done

We should generalize this a bit too. Let’s make this a script to which we can pass our filenames and ranges:

#!/bin/bash
namer=`basename $1 .bpf`
for START in $2
 do
  nameend=`echo $START | sed s/:/-/g`
  pdal translate $namer.bpf $namer"_"$nameend".bpf" -f range --filters.range.limits="Height[$START)"
done

Which to run, we use a statement as follows:

./tree_slicer.sh 54001640PAN_heightasz.bpf "0:0.5 0.5:1 1:2 2:5 5:10 10:15 15:20 20:35 35:50 50:75 75:100 100:200"

I forgot all my data is in feet, but my height classes in meters, so we’ll just multiply our input values by a factor of 3 to get in the same ballpark (future refinement likely):

./tree_slicer.sh 54001640PAN_heightasz.bpf "0:1.5 1.5:3 3:6 6:15 15:30 30:45 45:60 60:105 105:150 150:200"

We could alternatively stick to our original categories (short shrub, tall shrub, sub-canopy, canopy) as our break points:

./tree_slicer.sh 54001640PAN_heightasz.bpf "0:1.5 1.5:6 6:30 30:45 45:200"

Finally, we can convert to laz, and load all our slices of point clouds in plas.io an animate between the slices

for OUTPUT in $(ls *.bpf); do docker run -v /home/gisuser/test/test:/data pdal/master pdal translate //data/$OUTPUT //data/$OUTPUT.laz; done

If you look closely, you should be able to see where a tornado in the 80s knocked down much of the forest here. It’s signature is in the tall shrub / sub-canopy layer:

Posted in 3D, LiDAR, Other, PDAL, pointcloud | Tagged: , , | 2 Comments »

Point Clouds in Sketchup?

Posted by smathermather on February 4, 2016

Wow. Thanks to my colleague Brandon Garman for this tip — Point Clouds in Sketchup. Mind you, I don’t think this is an inexpensive proposition, but a pretty cool looking tool:

Posted in 3D, LiDAR, Other | Tagged: , | Leave a Comment »

Taking Slices from LiDAR data: Part II

Posted by smathermather on February 3, 2016

Ok, with a little help from Bradley Chambers on the PDAL mailing list, we are back in business. If we want to filter our newly calculated heights into a new PDAL output, we can do that easily, say all points 100-500 above ground level:


pdal translate 54001640PAN_heightasz.bpf 54001640PAN_heightasz_gt100.bpf -f range --filters.range.limits=&quot;Height[100:500]&quot;

A little sanity check to see if we are getting appropriate values:

pdal info --stats 54001640PAN_heightasz_gt100.bpf --dimensions &quot;Height&quot;
{
  &quot;filename&quot;: &quot;54001640PAN_heightasz_gt100.bpf&quot;,
  &quot;pdal_version&quot;: &quot;1.1.0 (git-version: 64c722)&quot;,
  &quot;stats&quot;:
  {
    &quot;statistic&quot;:
    [
      {
        &quot;average&quot;: 105.8738232,
        &quot;count&quot;: 179909,
        &quot;maximum&quot;: 194.5800018,
        &quot;minimum&quot;: 100,
        &quot;name&quot;: &quot;Height&quot;,
        &quot;position&quot;: 0
      }
    ]
  }
}

Ok, now I want to view this. I could convert to a *.laz file and view it with plas.io (as long as I use Chrome as my browser).
I have to switch back to docker, ’cause that’s where I have PDAL built with laszip:

docker run -v /home/gisuser/test:/data pdal/master pdal translate //data/54001640PAN_heightasz_gt100.bpf //data/54001640PAN_heightasz_gt100.laz

And now I can view in plas.io:

Screen Shot 2016-02-03 at 9.57.23 PM

Posted in 3D, Other, PDAL | Tagged: , , | Leave a Comment »