Monday, February 27, 2012

C# code snippet to determine if a point is in a polygon

A common test in GIS is to determine whether a point is inside a polygon or not.
The following C# code snippet can determine whether a point is inside a simple 2-D polygon. Just pass in an array of the polygon vertices and the point. The function returns true if the point is in the polygon and false if it is not.

private bool IsPointInPolygon(PointF[] polygon, PointF point)
bool isInside = false;
for (int i = 0, j = polygon.Length - 1; i < polygon.Length; j = i++)
if (((polygon[i].Y > point.Y) != (polygon[j].Y > point.Y)) &&
(point.X < (polygon[j].X - polygon[i].X) * (point.Y - polygon[i].Y) / (polygon[j].Y - polygon[i].Y) + polygon[i].X))
isInside = !isInside;
return isInside;

Monday, February 20, 2012

Simple method to count trees using Saga GIS

It is possible to make a rough estimation of the number of trees in an area from LiDAR derived digital surface (DSM) and digital terrain models (DTM). One method is to use some of the grid analysis modules algorithm in SAGA GIS, such as Gaussian Filter, and Watershed Segmentation. Then simply count the number of segmented table records with height greater than a value.

The example here counts the trees using the following general steps:

  1. Load the DSM and DTM datasets
  2. Calculate the canopy heights
  3. Smooth the canopy heights
  4. Segment the canopy heights
  5. Count the number of segments with canopy heights above a certain value
Load the source datasets
  1. Start SAGA GIS.
  2. Load and display the digital surface model (DSM) grid file, e.g. C:\data\dsm.asc.

  3. Load and display the digital terrain model (DTM) grid file, e.g. C:\data\dtm.asc.

Calculate the canopy heights
  1. Select Modules | Grid | Calculus | Grid Difference.

    The Grid Difference dialog box appears.
  2. In the Grid system field, choose the system of the source datasets, e.g. 683x 683y; 312480x 51952717y.
  3. In the A field, choose the digital surface model grid, e.g. dsm.
  4. In the B field, choose the digital terrain model grid, e.g. dtm.
  5. Click Okay.

    The canopy height model grid is created. The default name is Difference(A-B).
Smooth the canopy heights
  1. Select Modules | Grid | Filter | Gaussian Filter.

    The Gaussian Filter dialog box appears.
  2. In the Grid system field, choose the system of the source grids e.g.  683x 683y; 312480x 51952717y.
  3. In the Grid field, choose the canopy height grid, e.g. Difference (A-B).
  4. In the output Filtered Grid field, choose Create.
  5. In the Search Radius field, choose a value to approximate the tree radius e.g. 5.
  6. Click Okay.

    The output smoothed grid is created. Default name is Difference (A-B) [Gaussian Filter].
Segment the smoothed canopy height model
  1. Select Modules | Imagery | Segmentation | Watershed Segmentation.

    The Watershed Segmentation dialog box appears.
  2. In the Grid system field, choose the source grid system e.g.  683x 683y; 312480x 51952717y.
  3. In the Grid field, choose the smooth canopy height grid e.g. Difference (A-B) [Gaussian Filter].
  4. In the Output field, choose Seed Value.

    Note: since the segmentation will be run on the source grid of canopy height values, the resultant seed values will be the estimated heights of the tree tops.
  5. Ensure the Method is set to Maxima.

    Note: if set to Maxima, the watershed segmentation algorithm will start from the top (max) of the tree and flow down.
  6. Optional. Toggle on Borders.
  7. Click Okay.

    The segmented canopy height grid "Difference (A-B) [Gaussian Filter][Segments]" is created. 

    The canopy height point layer "Difference (A-B)[Gaussian Filter][Seed]" is created also.
Count the number of segments above a certain canopy height
  1. In the Data tree pane, right click on the newly created point layer node under the tree hierarchy - Shapes | Point | 01. Difference (A-B) [Gaussian Filter][Seed].

    A pop up menu appears.
  2. Choose Attributes | Show.

    The attributes table is displayed.
  3. Right click on the VALUE header.

    A pop up menu appears.
  4. Choose Sort Fields.

    The Sort Table dialog box appears.
  5. In the Sort first by field, choose VALUE.
  6. In the Direction field, choose descending.
  7. Click Okay.

    The table is sorted.
  8. Scroll down the table to the rows with the VALUE greater or equal to 5.

    The row number on the left indicates the rough estimate of the number of trees greater than 5 meters.

    Note: instead of using SAGA GIS to do the counting with the tabular data, it may be easier to export out the point layer into a database software and perform the analysis there

Monday, February 13, 2012

Simple canopy height derivation using Global Mapper

The forest canopy height is the height of the highest vegetation above the ground level. It is possible to calculate the canopy height from LiDAR derived digital terrain models (DTM) and digital surface models (DSM) gridded datasets by doing a simple subtraction in Global Mapper. The example below shows how to use the Combine Terrain Layers function in Global Mapper to calculate the canopy height.

  1. Start up Global Mapper. Load in a digital terrain model ESRI ArcGrid file e.g. dtm.asc.

  2. Load in a digital surface model ESRI ArcGrid file e.g. dsm.asc.

  3. Select File | Combine Terrain Layers.

    The Combine Terrain Options dialog box appears.
  4. In the Description field, type in relevant text e.g. CanopyHeight.
  5. In the Operation field, choose Subtraction (Difference) - Signed.
  6. In the Select First Elevation Layer to Combine field, choose the DSM layer e.g. dsm.asc.
  7. In the Select Elevation Layer to Combine the First Layer With field, choose the DTM layer e.g. dtm.asc.

  8. Click OK.
  9. Press CTRL+C. Toggle off the DTM and DSM layers.

    The canopy height is displayed as shown below where the pink indicates zero canopy heights and the blues are positive heights. 

Monday, February 6, 2012

Google Mapplet for State Plane Coordinate System readout

The State Plane Coordinate System (SPCS) is a map coordinate system commonly used in the United States. Most SPCS zones are based on the Transverse Mercator, Lambert Conformal Conic or Hotine Oblique Mercator projections. For more information about SPCS, visit this site at
This Google Mapplet can be used to display any State Plane Coordinate System coordinates. You do have to choose which coordinate system to use e.g. Alabama East Zone FIPS101 and the NAD83 datum. Then simply click anywhere on the map and the mapplet will display the SPCS easting, northing coordinates as shown below.
The datum transformation between the WGS84 datum used by Google Maps and the local NAD27 or NAD83 projection datum is done using the Standard Molodensky transformation. Parameters for the datums and ellipsoids are taken from the University of Colorado.

Click this link to run the mapplet.