Friday, October 30, 2009

Flagging intersection points between 2 features with FME

It seems to be easier to fix rather than flag unbroken intersection points between two feature layers. It took me a little while longer to flag the unbroken intersection points using FME transformers. The screenshot below shows two shape file layers - lakes (blue color fill) and rivers (red). I wanted to flag the intersecting points between the lake polygons and the river lines.

I did the following.


Define the two source datasets
  1. Start up FME Workbench. Choose Create a blank workspace.

  2. Select Source Data > Add Dataset.

    The Add Source Dataset dialog appears.
  3. In the Format field, choose ESRI Shape. In the Dataset field, click the browse [...] button.



  4. Select both the lake.shp and river.shp files. Click Open.

    The Add Source Dataset may look like this.


  5. Click OK.

    The Select Feature Types dialog box appears.


  6. Click OK.

    The two selected shape files are added to the workspace.

Define the destination dataset
  1. Select Destination Dataset > Add Dataset.

  2. In the Format field, select ESRI Shape. In the Dataset field, define the output folder e.g. C:\temp\out\.

    The Add Destination Dataset dialog box may look like this.
  3. Click OK.

    The message below appears.


  4. Click Yes.

    The Feature Type Properties dialog box appears.

  5. In the Feature Type Name field, type in intersect_flag. In the Allowed Geometries, choose shape_point.

    The Feature Type Properties dialog box may look like this.


  6. Click User Attributes tab.

  7. Define an attribute ID with Data Type number and Width 2.

    The Feature Type Properties dialog box may look like this.


    Note: id will simply be a counter storing running numbers of flags.

  8. Click OK.

    The source and destination datasets appear on the workspace.
Add in the transformers and define the flow
  1. Drag and drop the following transformers onto the workspace:
    1. Intersector
    2. Coordinate Fetcher (twice)
    3. Geometry Remover (twice)
    4. Aggregator
    5. Tester
    6. 2DPointAdder
    7. Counter

  2. Connect the source and destination datasets with the transformers as shown in the figure below.



  3. Open up the first Coordinate Fetcher transformer's properties.



    Note 1: Index 0 means fetch the first coordinate.
    Note 2: After the Intersector transformer has done its job, the output features are broken up at the point of intersections. The intersection points must therefore be at either ends.

  4. Open up the second Coordinate Fetcher transformer's properties.



  5. In the Index field, type in -1.

    Note: -1 means fetch the last coordinate.

  6. Open up the Aggregator transformer's parameters.



  7. In the Group By field, browse and select the fields or type in the string _x _y.

  8. In the Count Attribute (Optional) field, type in _count.

  9. Open up the Tester transformer's parameters.



  10. Define the test _count > 2 as shown above.

  11. In the Pass Criteria field, choose One test must pass.

    Note: At the intersection points, there must be more than two points. Anything less is either free end points or singly connected points. 

  12. Open the 2DPointAdder transformer's parameters.

  13. In the X Value field, choose _x. In the Y Value field, choose _y.

  14. Open up the Counter transformer's parameters.



  15. In the Count Output Attribute field, type in id. In the Count Start field, type in 1.

    Note: the Count Output Attribute field should be the same name as the destination dataset's attribute field.

  16. Run the translation.

    The intersection points are flagged as numbered points as shown in the screen shot below.

Note that this Intersector transformer method works based on line work intersection only. If the line work sits entirely within the area polygons, there would be no line intersection and therefore no flags.

Thursday, October 22, 2009

Fixing linework intersections using FME's Intersector transformer



I am looking into using FME for performing data validation and quality assurance of some digitized line work. I found a  rich set of transformers for this purpose. One of the data cleaning task that I want to do is to break the line work at the point of intersection. The transformers that can do this are the Intersector and the MRFClean transformers. A more specialized variant of the Intersector is the SelfIntersector transformer.

I gave the Intersector transformer a try on a test case shape file dataset.. I put in a few cases: (1) a self-intersecting line, (2) two intersecting lines, and (3) two overlapping or duplicate lines as shown below.

The Intersector transformer has a few options. For self-intersecting line work, you have the option of getting the transformer to break the line work for you or you can leave it as it is. I would prefer to delegate the breaking of the self-intersecting line work to the SelfIntersector transformer as I would be able to differentiate it from the other cases.

For line work with overlapping segments, the Intersector transformer has the option to output as a merged line work or as separate segments. I would prefer the merge option.

To fix the digitized line work using the Intersector transformer, I used the following steps.
  1. Start up the FME Workbench and choose to use the Workspace Wizard.

  2. In the Wizard, specify ESRI Shape as the source format and choose the input shape file.

  3. Specify ESRI Shape as the destination format

  4. Drag the Intersector transformer from the Transformers pane onto the Main pane. Connect the transformer with the source dataset and the dataset as shown in the screen shot below.



  5. Open up the Destination dataset's Feature Type Properties.Click User Attributes tab. Add in 2 number attributes - _overlaps and _segments as shown below. Click OK.


  6. Open up the Intersector's Parameters. Review the parameters. Click OK.



    Note 1: If Split Self-Intersecting Features is set to No, then self-intersecting line work will not be fixed.
    Note 2: If Separate Collinear Segments is set to No, then overlapping segments in the line work will be fixed and merged
    .

  7. Run the translation.

    The dataset is transformed and fixed.
The results are shown in the screenshots below. 
For the overlapping segments, the Intersector has merged into a single geometry. The _overlaps attribute is set to 2 since 2 overlapping segments were merged into 1. 


For intersecting line work, the Intersector transformer has broken them up into separate segments at the point of intersection. For each segment, the _segments attribute is set to 2 for the number of segments created.



The self-intersecting geometry is left as it is - unfixed because the relevant Intersector parameter was set to No.
The Intersector transformer works quite well for my purpose. It also works on polygon line work as well as 3-D datasets. For 3-D line work, it will optionally create one or two nodes at the point of intersection. The only thing left for me to do is to create some mapping file and scripts to generalize the source and destination dataset definitions.

Wednesday, October 14, 2009

Find free end points with GeoMedia


I wanted to find free end points of some linear geometry features using only GeoMedia or GeoMedia Professional's built-in tools without resorting to using third party applications like MRF Clean. An example linear feature set is shown below. Note that the line work is already broken up at the intersection.

I tried GeoMedia Professional's Validate Connectivity command but it doesn't quite give me a completely satisfying result. After some trial and error, I finally worked out a solution using GeoMedia's Functional Attributes, Union, Spatial Query, Attribute Query, Spatial Intersection, and Spatial Difference commands. In a nutshell, the idea can be broken down into the following parts:
  • Get a list of all the starting and ending points of the line geometries
  • Use Spatial Intersection to determine which points are intersecting with which points
  • Use Spatial Difference command to subtract away the points connected to 2 or more line segments to get the free end points
Here are the steps.

Find all the end points
  1. In GeoMedia, select Analysis > Functional Attributes.

    The Functional Attributes dialog box appears.


  2. In the Add functional attributes for combo box, choose the linear feature e.g. rcl. Click New.

    The Functional Attribute dialog box appears.

  3. In the Categories list, click Geometry. Then in the Functions list, double click STARTPOINT. In the Attributes list, double click the linear feature's geometry field e.g. Input.Geometry.

    In the Expression field, STARTPOINT(Input.Geometry) should be displayed.

  4. In the Functional attribute name field, type in EndpointGeometry.

    The Functional Attribute dialog box should be similar to the figure below.


  5. Click Add.

    The EndpointGeometry functional attribute is added.

  6. Click Close.

    The Functional Attribute dialog box is closed.

  7. In the Query name field, type in Start of rcl. Toggle off Display functional attributes in data window as shown below..



  8. Click OK.

    The Start of rcl query is displayed in the map window.


  9. Now repeat the previous steps to find the linear geometries' end points. Select Analysis > Functional Attributes.

    The Functional Attributes dialog box appears. The previous linear feature rcl and its EndpointGeometry functional attribute should still be selected.


  10. Click Properties.

    The Functional Attribute dialog box appears.


  11. Clear the Expression field.

  12. In the Functions list, double click ENDPOINT. In the Attributes list, double click the linear feature's geometry field, e.g. Input.Geometry.

    The string ENDPOINT(Input.Geometry) should be displayed in the Expression field.


  13. Click OK.

  14. In the Query name field, type in End of rcl. Click OK.

    The End of rcl query is displayed in the map window.
  15. Now select Analysis > Union.

    The Union dialog box appears.

  16. In the Union features in tree view, expand the Queries node and toggle on Start of rcl and End of rcl. In the Query name field, type in Endpoints of rcl.

    The Union dialog box should look similar to the figure below.


  17. Click OK.

    The Endpoints of rcl query is displayed.


We have now created a list of all the end points for the linear geometry feature. From this list we can use the Spatial Intersection command to find out which points are intersecting with which points. That's described below.

Find out which end points are intersecting with which end points
  1. Select Analysis > Spatial Intersection.

    The Spatial Intersection dialog box appears.

  2. In the Generate intersections for Features in combo box, choose the Endpoints of rcl query.

  3. In the That combo box, choose are spatially equal.

  4. In the second Features in combo box, choose the Endpoints of rcl query again.

    The Spatial Intersection dialog box may look like this now.


  5. Click OK.

    The spatial intersection query is created.


    Note 1: The Data Window shows the attributes of the spatial intersection query result. The example linear feature rcl has the field 'id' as the primary key.
    Note 2: Since both input features to the Spatial Intersection command are the same feature and hence the same primary key 'id', GeoMedia will rename the second input feature's primary key by appending a number '1' to it resulting in 'id1'.

A spatial intersection record is the result from two overlapping geometry features. In this example, GeoMedia will append the records from both feature (the same rcl feature) into one intersection record. Therefore, we can say for sure that at the free end points, the spatial intersection record's id and id1 values would be the same. We can also say for sure that at the line intersections, the spatial intersection records' id and id1 values may be the same or may not be the same. Given these two facts, to find out the free end points, all we have to do is to spatially subtract the points at the line intersections from the list of all end points. That is described below.

Use Spatial Difference to subtract away the points at the line intersections
  1. Select Analysis > Spatial Difference.

    The Spatial Difference dialog box appears.


  2. In the From features in combo box, choose the query Spatial Intersection of Endpoints of rcl and Endpoints of rcl.

  3. In the Subtract features in combo box, choose the same query Spatial Intersection of Endpoints of rcl and Endpoints of rcl.

  4. In the Subtract features in section, click Filter.

    The Spatial Intersection of Endpoints of rcl and Endpoints of rcl filter dialog box appears.

  5. Either type in or use the controls to form the string "id <> id1" in the Filter field as shown below.

  6. Click OK.

  7. In the Output difference as Query name field, type in "Free Endpoints of rcl".

    The Spatial Difference dialog box may look similar to this.


  8. Click OK.

    The free end points of the linear geometry feature are displayed in the map window.


Tuesday, October 6, 2009

GeoMedia 6.1's new Export Layout to GeoTiff command

There is a new hot fix 6.1.6.19 out.for GeoMedia and GeoMedia Professional 6.1. One of the new enhancement available in the hot fix is the ability to export the Layout window into a GeoTiff file. I tried out the new Export Layout to GeoTiff command with one of the sample GeoWorkspace files that has some Layouts defined. Here is how it works.
  1. Startup GeoMedia or GeoMedia Professional and open up a GeoWorkspace e.g. BatchPlotting.gws.

    The workspace is displayed in GeoMedia.

    Note: The red rectangle is approximately the area covered by the layout window.

  2. Select Window > Show Layout Window.

    The Layout Window appears.


  3. Select Sheets > Export Layout.

    The Export Layout Sheet As dialog box appears.

    Note: if you click the Save as type drop down list, you can see the new GeoTiff option.

  4. In the File name field, type in the name of the new file, e.g. south_hsv.tif. Choose Georeferenced TIFF - GeoTIFF (*.tif) in the Save as type field. Click Save.

    The GeoTIFF Export Options dialog box appears.
  5. Change the File resolution to your desired value. Click OK.

    The layout is exported to a GeoTIFF file.
I displayed the resultant layout GeoTIFF file in Global Mapper 10 successfully. The results can be seen in the screenshot below. I overlay the County boundary for reference and we can see that the layout GeoTIFF image sits correctly with respect to the County boundary; compare that with the red rectangle location in the first screenshot at the top.