Earlier, I posted some simple pictures from this short video clip using LiDAR to determine vegetation shape, animated in PovRay. Now, I post the video itself (all 7 seconds):
I did what YouTube says not to do, and that is to upload a vignetted version of video. Such is. I’ll upload a better one shortly.
I won’t offer much prelude. Read this post, and this post first… .
Inadvertently I demonstrated the value of spatial indices, i.e. I meant to use them and didn’t. In attempting to sort my tables by their spatial index I got the following errors:
ALTER TABLE lidar_ground CLUSTER ON lidar_ground_the_geom_idx; ERROR: "lidar_ground_the_geom_idx" is not an index for table "lidar_ground" ALTER TABLE lidar_veg CLUSTER ON lidar_veg_the_geom_idx; ERROR: "lidar_veg_the_geom_idx" is not an index for table "lidar_veg"
Looking further, I saw there were no spatial indices on either table.
Returning to my original post, I had typos:
CREATE INDEX lidar_veg_the_geom_idx ON lidar USING gist(the_geom); CREATE INDEX lidar_ground_the_geom_idx ON lidar USING gist(the_geom);
Do you see it? It should be
CREATE INDEX lidar_veg_the_geom_idx ON lidar_veg USING gist(the_geom); CREATE INDEX lidar_ground_the_geom_idx ON lidar_ground USING gist(the_geom);
I created two identical indices on the lidar table, not on the two tables I was interested in. So, what was the relative performance hit?
Without the indices:
Vinton=# EXPLAIN ANALYZE SELECT pgis_fn_nn(b.the_geom, 10,1,1,'lidar_ground','true','vnum','the_geom') FROM lidar_veg AS b WHERE b.vnum > 576505 AND b.vnum < 576600; QUERY PLAN ---------------------------------------------------------------------------------------------------------------------------------------- Index Scan using lidar_veg_pkey on lidar_veg b (cost=0.00..108.97 rows=31 width=100) (actual time=571.154..28589.303 rows=50 loops=1) Index Cond: ((vnum > 576505) AND (vnum < 576600)) Total runtime: 28589.398 ms (3 rows)
With the indices:
EXPLAIN ANALYZE SELECT pgis_fn_nn(b.the_geom, 10,1,1,'lidar_ground','true','vnum','the_geom') FROM lidar_veg AS b WHERE b.vnum > 576505 AND b.vnum < 576600; QUERY PLAN ------------------------------------------------------------------------------------------------------------------------------ ------- Index Scan using lidar_veg_pkey on lidar_veg b (cost=0.00..108.97 rows=31 width=100) (actual time=86.489..138.579 rows=50 lo ops=1) Index Cond: ((vnum > 576505) AND (vnum < 576600)) Total runtime: 138.717 ms (3 rows)
Now with CLUSTERing:
ALTER TABLE lidar_veg CLUSTER ON lidar_veg_the_geom_idx; ALTER TABLE lidar_ground CLUSTER ON lidar_ground_the_geom_idx; EXPLAIN ANALYZE SELECT pgis_fn_nn(b.the_geom, 10,1,1,'lidar_ground','true','vnum','the_geom') FROM lidar_veg AS b WHERE b.vnum > 576505 AND b.vnum < 576600; QUERY PLAN ------------------------------------------------------------------------------------------------------------------------------ ----- Index Scan using lidar_veg_pkey on lidar_veg b (cost=0.00..108.97 rows=31 width=100) (actual time=1.929..56.285 rows=50 loop s=1) Index Cond: ((vnum > 576505) AND (vnum < 576600)) Total runtime: 56.378 ms (3 rows)
Indices sped up the query by a factor of 206 times. Adding the CLUSTER sped it up another 2.4 times over and above the indices for a total of a query that runs 507 times faster. Did I say the full query would take a few years? Ahem, I meant a couple days. Specifically about 2.
Something’s been nagging at me as far as the level of optimization (or lack there-of) that I did to my database before my other posts of using PostGIS to analyze LiDAR data (e.g. this post). It seemed my results were remarkably slow, but I couldn’t put my finger on why that was problematic.
Then, as I was testing a smaller lidar dataset in order to do a shading analysis for my future garden, I noticed that with those 1800 points (a high density point cloud considering I live on a quarter-acre), the nearest-neighbor algorithm, per point, was much faster than when applied to the larger dataset I was testing before. On one hand this seems reasonable. A bigger dataset takes longer to process. But, properly optimized, a large dataset should not take significantly longer per record. We want linear, or near linear, scaling of processing. When we don’t, in this case, it indicates that the data are not structured in a way that allows for fast access to records of interest.
Re-ordering records (PostgreSQL CLUSTER) to reflect the hierarchy of an index is a good way to optimize a database, assuming you know which index is your most important one to use to organize your data structure. In the case of large spatial datasets, this is usually the spatial index.
In this example, I will apply CLUSTERing to the PostGIS database, clustering on the spatial index in order to optimize spatial queries.
ALTER TABLE lidar_ground CLUSTER ON lidar_ground_the_geom_idx; ALTER TABLE lidar_veg CLUSTER ON lidar_veg_the_geom_idx;
I will report back with some results. Hopefully this helps.
I ran these tests a while back, but I’m just getting to post them now. These were run on a late model MacBook Pro, with a 2.4GHz Intel Core 2 Duo, with a couple of gigs of RAM on PostgreSQL. In Postgre, if I run
SELECT version();
I get:
PostgreSQL 8.4.1 on i386-apple-darwin10.0.0, compiled by GCC i686-apple-darwin10-gcc-4.0.1 (GCC) 4.0.1 (Apple Inc. build 5493), 64-bit
My hard drive is as follows:
Vendor: NVidia Product: MCP79 AHCI Speed: 3 Gigabit Description: AHCI Version 1.20 Supported
Ok. So enough prelude. Back in New Year’s day post, I used the nearest neighbor function posted on the Boston GIS page to find the nearest ground point to one of my vegetation points in an attempt to derived vegetation height within a LiDAR dataset loaded into PostGIS. I got as far as identifying the nearest neighbor, and then stopped. What next? Well, how fast (or slow) is that nearest neighbor search?
Vinton=# EXPLAIN ANALYZE SELECT pgis_fn_nn(b.the_geom, 10,1,1,'lidar_ground','true','vnum','the_geom') FROM lidar_veg AS b WHERE b.vnum > 576505 AND b.vnum < 576600; QUERY PLAN ---------------------------------------------------------------------------------------------------------------------------------------- Index Scan using lidar_veg_pkey on lidar_veg b (cost=0.00..108.97 rows=31 width=100) (actual time=571.154..28589.303 rows=50 loops=1) Index Cond: ((vnum > 576505) AND (vnum < 576600)) Total runtime: 28589.398 ms (3 rows)
I also timed it with a stop watch. No matter how I changed the search parameters for the nearest neighbor function, it took almost exactly 35 seconds. Extrapolating upward, that means that to do this for one of my blocks of data would take an excess of 50 hours. I have 559 of these blocks of data, so we’re looking at ~28,000 hours of computing time, or about 3.2 years. It’s my wife’s laptop, and I don’t think I should take it over for that long… .
So, with no real knowledge of how to optimize things, I’m thinking instead of a nearest neighbor approach, perhaps a point-in-polygon analysis with Voronoi polygons created from my LiDAR ground points. I can use R to generate my polygons, ala (yet again) Boston GIS, and then ST_Contains to do the spatial join. Hopefully by already having the geometry constructed with R, it will be a lot more efficient that expanding buffers from points.
I wanted to put together a quick and dirty map of a biological reserve in Ecuador, sort of a laptop exploration of a place quite distant. At first, I thought I’d use Shuttle Radar Topography Misssion data to get the elevation information. Then I discovered the ASTER Global DEM which is 30m resolution for the whole world. Wow. Cool cool data. (I used the Japanese site to download, as I had trouble getting the US site to load).
So, first I downloaded and extracted the zipped geotiffs that come from the ASTER Global DEM. I extracted the DEMs into the same directory, and mosaicked them into a single DEM:
gdalwarp AST*.tif ecuador_aster_dem.tif
Then I created a hillshade version of the DEM for mapping:
gdaldem hillshade ecuador_aster_dem.tif ecuador_aster_hillshade.tif
and then 100 and 50 meter contours for the site:
gdal_contour -a elev -3d -i 100 ecuador_aster_dem.tif ecuador_aster_contour_100.shp gdal_contour -a elev -3d -i 50 ecuador_aster_dem.tif ecuador_aster_contour_50.shp
I loaded it all into Quantum GIS, symbolized it, and added some OpenStreetMap data. To do this, I navigated to the location in OSM, chose the download tab, and chose OpenStreetMap XML Data as my output. Using the OSM plug-in in QGIS, I imported the file, and walla– a nearly functional map.
It finished, yay! A little Friday 5PM exuberance, I’m afraid. See previous post for explanation.
Here is the clipped version:
I loaded a set of contours into the database that are really nice, up-to-date LiDAR and breakline derived contours. The are engineering grade aerial contours, but they are a very big and complex dataset. So I’ve done what I had hoped was the hard part, using shp2pgsql, I’ve converted them to PostGIS insert statements and dumped them into the database only to find that some sections are corrupt. I’m hoping to use the function ST_Difference to clip out the bad areas, so I can insert them in again.
There are two problems in the contour dataset– some areas are duplicated records, and some areas are missing or otherwise corrupted areas. I’ve got a polygon file (in purple below) that delineates where the problem areas are, so I hope to use the function ST_Difference to remove these areas so I can repopulate with good data.
CREATE TABLE public.cuy_contours_2_clip
AS
SELECT elevation,
ST_Difference(a.the_geom, b.the_geom) AS the_geom
FROM base.cuy_contours_2_1 as a,
public.contour_clip_bounds as b;
Let’s see if PostGIS can handle this big of a dataset… .
I completed this project a long time ago and have not blogged on it in over a year. None-the-less, I realized that I hadn’t posted any final images for the viewshed analysis with Povray.
So here is the summer viewshed analysis. This is rendered with about (I think) 150,000 trees, each with 100,000+ leaves. The cell tower here is in red. The landscape colored in cyan where the tower can be seen from. A property boundary shows up in green:
And here is the same viewshed but a winter scene. Since this is a very flat wetland landscape, it’s important to understand the visual penetration of the cell tower with the vegetation, with and without leaves:
IMCORR is a package distributed by the National Snow and Ice Data Center (did you know we have one of those?), or NSIDC, that performs image cross correlation between two images using a comparison between a moving image chip in each image. It measures the displacement (in pixels) between the objects found in the two images, and writes that out to a text file. It is used to calculate how quickly glaciers are moving by calculating displacements. If you have an image of a glacier from two different time periods, this little program will calculate how far crevasse zones and other discernible features have moved between two images. I know I know. I have been talking about trees, and now ice. Bear with me.
IMCORR can also be used for simple image georeferencing, say when we don’t want to collect our own tie points. We can look for correlations between two images, and tabulate those as points, use them with a thin plate spline, in say gdalwarp, and thus correct for distortions in one dataset with another one.
Now, this is all in theory, and I haven’t tested it yet. But I’m thinking this may be a really useful technique, especially because we already have well referenced high-resolution aerials for so many areas. I have an effective 3″ color infrared (CIR) dataset for an entire county kicking around somewhere that has been referenced to frame on center, but with no terrain correction or other optical distortion corrections applied. I’m hoping to cross correlate these with a 6″ georeferenced aerial in order to make these data useable. 3″ CIR could be very useful for all sorts of fun applications… .
Just another shot of our digital surface model of the forest. This time, I rendered it with all of the trees at base elevation 0 (zero), placed the orthographic camera at height 150, so our 16-bit range we’re rendering to should scale 0-150 feet to 0-65535 values, giving us a precision of 0.002 feet. The original LiDAR data has a precision of 0.01 feet, so we’re not throwing data away with the integerization. By the way, we don’t have enough LiDAR points to really derive canopy shape– we’ve invented the detailed canopy shape by using a tree shaped interpolation technique (if we can quite liberally call it interpolation– I’d hate to see the actual function written out).
Anyway, the pretty pic– the darker the green, the taller the tree– we max out around 120 feet in most areas. We do find at least one 149 foot tall tree, but I’m guessing that’s on a steep slope, and is an over estimate due to a tree leaning of spreading over a cliff, or some such scenario. Of course, if it does exist, I should go look for it. I did measure one in the field last year that exceeded 130 feet, and I couldn’t spot the top, so it may have been much taller.
Oh, and Mac bragging rights here– I couldn’t render the full 30,000 x 30,000 pixel scene on Windows XP (although Vista may have handled it), but it rendered just fine in 35 minutes on my Mac laptop.
