Monday, July 29, 2013

Use Saga GIS to create Voronoi or Thiessen polygons

The new version of Saga GIS 2.1 has a functional Thiessen polygon command. Previously, a bug in the function creates erroneous polygons. Here is an example of using the command to create Voronoi or Thiessen polygons from discrete points with some polygon clipping thrown in for aesthetic effects. The example shows how to create polygons of cities in the great State of Texas.

  1. Start SAGA GIS. Load and display the cities and the state boundary polygon in a map view.

  2. Select Modules | Shapes | Points | Thiessen Polygons.

    The Thiessen Polygons dialog box appears.

  3. In the input Points field, choose the point layer, e.g. 01. cities_in_texas.
  4. In the output Polygons field, select create. Click Okay.

    The Thiessen polygons layer (02. cities_in_texas [Thiessen Polygons] is created. The polygons are spread out over some arbitrary four points.

  5. Select Modules | Shapes | Polygons | Polygon Clipping.

    The Polygon Clipping dialog box appears.

  6. In the Clip Features field, select the State boundary polygon layer, e.g. 01. texas.
  7. Select the Input Features field and click the [...] button.

    The Input Features dialog box appears.

  8. In the left list, select the Thiessen Polygons layer, e.g. 02. cities_in_texas [Thiessen Polygons]. Click the > button. Click Okay.

    The Input Features dialog box is closed.
  9. Click the Okay button.

    The Thiessen polygons are nicely clipped to the State polygon

Monday, July 22, 2013

Create cloud masks for aerial photos in batch

Clouds in aerial photos can interfere in the feature identification and matching stages of the ortho-photo rectification process. Depending on the software algorithm, it may be necessary to either remove the entire photo or mask out the cloud portions from the photo. If there are hundreds or thousands of photos with clouds, then it will be very tedious to manually digitize the cloud masks with a photo editing software. It would be nice to have some tool to automate the task somewhat.

Here is a simple procedure of using the free and open source software Imagemagick and the threshold option to generate the cloud masks from the photos. The assumption is that clouds in aerial photos would be brighter than the terrain underneath. An example aerial photo with clouds is shown below.
Photo credit: Kevin Payravi, Wikimedia Commons

  1. In the Windows Command Prompt, type in the following:

    C:/> convert Overview_head_of_cloud_with_shadow.jpg -threshold 40% cloud_mask.png

    Note: the threshold option tells Imagemagick to set pixel values below 40% to 0 (black) and above to white.

    The mask is created.

  2. Optional. In some software, black is the mask so it may be necessary to invert the mask with the following command.

    C:/> convert cloud_mask.png -channel rgba -transparent black -fill black -opaque white -fill white -opaque none -alpha off -define png:color-type=6 cloud_invert_mask.png

    The mask is inverted as RGB file type
Note: if there are bright objects on the ground e.g. river, stream, road, metallic rooftops, then the mask may contain these objects too. It may be necessary to edit the mask to remove these objects.

Monday, July 15, 2013

Use Notepad++ for comparing lists

When processing hundreds or thousands of LiDAR LAS files or aerial photographs, sometimes one or more files may not get processed and missed out from the processing flow. It is useful to be able to identify which files did not get processed so that they can be reprocessed again. There are many ways to compare and one such method is using Notepad++, a free text editor. There is a free plug-in Compare for Notepad++,  that can be used for comparing lists to identify the differences, i.e. which file(s) were not created. The following illustrates how to use Notepad++ for identifying missing file(s) using before and after lists.

  1. Open up a Windows Command Prompt and use the DIR command to generate a list of files before and after processing, e.g. list.txt and newlist.txt.

    C:/> DIR /b *.las > list.txt
  2. Run Notepad++. Open up the lists.

  3. Select Plugins | Compare | Compare.

    The differences are identified.
  4. Scroll up or scroll down to see the differences. Or use the scroll buttons.

Monday, July 8, 2013

Google Gadget for monitoring 911 incident calls in Portland, Oregon

This gadget allows users to monitor the City of Portland, Oregon's Police 911 Dispatch Incidents on a Google Maps backdrop by subscribing to the geoRSS feed of the 100 most recent, closed, non-confidential calls for service received by the City of Portland's 911 system. 

The incidents are shown as clickable and color coded icons with tool tips. Incidents that happened an hour or less ago will be fully opaque; at more than an hour but less than three hours, the icons will be displayed at 50% opacity; at greater than three hours, the icons will be shown at 25% opacity. 

Leave the gadget on the screen and it will update itself every few minutes with the latest information. Clicking on the icon will bring up more details about the incident. Together with Google Street View, users can truly immerse themselves into the incident environment

Monday, July 1, 2013

C# QuadTreeLib example for indexing and querying point items

In GIS, quad trees are usually used for spatial indexing of geographical features. I found this great C# QuadTree Implementation (QuadTreeLib) library project at The library project comes with an example of indexing rectangles but not for points. So I thought of posting a simple example for points. In this example, the points in a LiDAR LAS file is read and spatially indexed using this library.

//include these name spaces
using System.Drawing;
using QuadTreeLib;
//Implement the QuadTreeLib's IHasPoint interface as my own class e.g. PointItem.
class PointItem : QuadTreeLib.IHasPoint
private PointF _point;
private uint _seq;
public PointF Point
get { return _point; }
public uint Seq
get { return _seq; }
public PointItem(double x, double y,uint seq)
PointF p = new PointF((float)x, (float)y);
_point = p;
_seq = seq;

//Declare and initialize the PointQuadTree for storing PointItem classes
private PointQuadTree<PointItem> quadTree = null;
static void Main(string[] args)


LASReader reader = new LASReader(@"C:\Temp\myfile.las");
LASHeader header = reader.GetHeader().Copy();

//Create a rectangle encompassing the dataset bounds
System.Drawing.RectangleF lasRect = new System.Drawing.RectangleF((float)header.GetMinX(), (float)header.GetMinY(), (float)(header.GetMaxX()-header.GetMinX()),(float)(header.GetMaxY()-header.GetMinY()));
//Create the PointQuadTree class
quadTree = new PointQuadTree<PointItem>(lasRect);

//Populate the quad tree with the points from the LiDAR file
uint counter = 0;
while (reader.GetNextPoint())
LASPoint lasPnt = reader.GetPoint().Copy();
PointItem p = new PointItem(lasPnt.X, lasPnt.Y, counter);
//Define a query window bounds
double xlo = 0;
double xhi = 100000;
double ylo = 0;
double yhi = 200000;

//Query the quad tree and get the points within the window
ArrayList candidates = GetCandidates(xlo, ylo, xhi, yhi);

public ArrayList GetCandidates(double xlo, double ylo, double xhi, double yhi)
ArrayList candidates = new ArrayList();
RectangleF window = new RectangleF((float)xlo, (float)ylo, (float) (xhi - xlo), (float) (yhi - ylo));

//Call the PointQuadTree's Query method to find all the points within the window
List<PointItem> list = quadTree.Query(window);
for (int i = 0; i < list.Count; i++)
PointItem p = (PointItem) list[i];
return candidates;