[RNFdev] Not a pretty picture

[Dan Putler]

Frank, it works as advertised. It also made me realize that the  
needed scripting tools for this undertaking probably should be Python  
rather than R focused since Python is more commonly used by folks in  
the open source GIS community (as opposed to people who primarily do  
statistics using spatial data).

Dan

On 25-Sep-06, at 5:22 PM, Dan Putler wrote:

> Thanks Frank, I'll give this a shot. I've got the most recent  
> version installed on my Linux box. It turns out that the gaps  
> between my hand created lines and other lines was done on purpose.  
> The goal was to avoid having edges cross one another.
>
> Dan
>
> On 25-Sep-06, at 5:08 PM, Frank Warmerdam wrote:
>
>> Dan Putler wrote:
>>> Hi Frank,
>>> I assume that a shapefile is a reasonable format. To that end,  
>>> I've attached the shapefile with the boundary linear features for  
>>> the K2W postal code.
>>> The linear features (road segments) that are not ordered within  
>>> the layer. Specifically, two adjoining road segments for the same  
>>> street do not follow one another in the original StatsCan RNF.  
>>> They are, in fact, randomly ordered. Probably to make things  
>>> render in a more attractive fashion. However, the first point of  
>>> one road segment is the same as the last point of the preceding  
>>> road segment for that street (road segments run intersection to  
>>> intersection). I did not use this information in my script (which  
>>> simply grabs the vertices of the line features), but I can see  
>>> that it could be used order the road segments. Where things break  
>>> down is for linear features that I have digitalized. They  
>>> represent non-road boundaries that are not in the RNF itself.  
>>> These features can be identified since they have a NAME attribute  
>>> of "Fill", and have a value for the RD_SEG field that is greater  
>>> than 3000.
>>
>> Dan,
>>
>> I have attached a script, mkpoly.py, which uses OGR to assemble  
>> the polygon.
>> I see there are gaps between some of your line endings (likely the  
>> manually
>> created ones?) so I have turned on some tolerances when assembling.
>>
>> If you install FWTools (http://fwtools.maptools.org) then from an  
>> FWTools
>> shell you should be able to do something like:
>>
>>   python mkpoly.py K2W.shp K2W_poly.shp
>>
>> It works ok for the dataset you sent, but we may need to fiddle  
>> with it
>> a bit for other cases.  We might also want to append the new  
>> polygon to
>> an existing master file with the FSA as an attribute.  All easily  
>> done
>> if the geometry assembly works adequately.
>>
>> Let me know if it works for you.
>>
>> Best regards,
>> -- 
>> --------------------------------------- 
>> +--------------------------------------
>> I set the clouds in motion - turn up   | Frank Warmerdam,  
>> warmerdam at pobox.com
>> light and sound - activate the windows | http://pobox.com/~warmerdam
>> and watch the world go round - Rush    | President OSGeo, http:// 
>> osgeo.org
>>
>> #!/usr/bin/env python
>> #******************************************************************** 
>> **********
>> #  $Id: gdalchksum.py,v 1.3 2005/02/23 14:37:03 fwarmerdam Exp $
>> #
>> #  Project:  GDAL
>> #  Purpose:  Merge a bunch of linestrings in into a single polygon.
>> #  Author:   Frank Warmerdam, warmerdam at pobox.com
>> #
>> #******************************************************************** 
>> **********
>> #  Copyright (c) 2006, Frank Warmerdam <warmerdam at pobox.com>
>> #
>> #  Permission is hereby granted, free of charge, to any person  
>> obtaining a
>> #  copy of this software and associated documentation files (the  
>> "Software"),
>> #  to deal in the Software without restriction, including without  
>> limitation
>> #  the rights to use, copy, modify, merge, publish, distribute,  
>> sublicense,
>> #  and/or sell copies of the Software, and to permit persons to  
>> whom the
>> #  Software is furnished to do so, subject to the following  
>> conditions:
>> #
>> #  The above copyright notice and this permission notice shall be  
>> included
>> #  in all copies or substantial portions of the Software.
>> #
>> #  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,  
>> EXPRESS
>> #  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF  
>> MERCHANTABILITY,
>> #  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO  
>> EVENT SHALL
>> #  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,  
>> DAMAGES OR OTHER
>> #  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,  
>> ARISING
>> #  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR  
>> OTHER
>> #  DEALINGS IN THE SOFTWARE.
>> #******************************************************************** 
>> **********
>> #
>>
>> import gdal
>> import ogr
>> import osr
>> import sys
>>
>> def Usage():
>>     print 'Usage: mkpoly.py srcfile dstfile'
>>     sys.exit(1)
>>
>> #  
>> ===================================================================== 
>> ========
>> # 	Mainline
>> #  
>> ===================================================================== 
>> ========
>>
>> srcfile = None
>> dstfile = None
>>
>> gdal.AllRegister()
>> argv = gdal.GeneralCmdLineProcessor( sys.argv )
>> if argv is None:
>>     sys.exit( 0 )
>>
>> # Parse command line arguments.
>> i = 1
>> while i < len(argv):
>>     arg = argv[i]
>>
>>     if srcfile is None:
>>         srcfile = arg
>>     elif dstfile is None:
>>         dstfile = arg
>>     else:
>>         print 'unknown option: ', arg
>>         Usage()
>>
>>     i = i+1
>>
>> if dstfile is None:
>>     Usage()
>>
>> #  
>> --------------------------------------------------------------------- 
>> -------
>> #	Read all the source linestrings.
>> #  
>> --------------------------------------------------------------------- 
>> -------
>>
>> src_ds = ogr.Open( srcfile )
>> src_layer = src_ds.GetLayer(0)
>>
>> collection = ogr.Geometry( type = ogr.wkbGeometryCollection )
>>
>> feat = src_layer.GetNextFeature()
>> while feat is not None:
>>
>>     collection.AddGeometry( feat.GetGeometryRef() )
>>
>>     feat.Destroy()
>>     feat = src_layer.GetNextFeature()
>>
>> print 'Got %d lines.' % collection.GetGeometryCount()
>>
>> #  
>> --------------------------------------------------------------------- 
>> -------
>> #	Assemble into a polygon.
>> #  
>> --------------------------------------------------------------------- 
>> -------
>>
>> try:
>>     poly = ogr.BuildPolygonFromEdges( collection, Tolerance=0.001,
>>                                       bAutoClose=1, bBestEffort=1 )
>>     print poly.ExportToWkt()
>> except:
>>     print 'BuildPolygonFromEdges failed.'
>>
>> #  
>> --------------------------------------------------------------------- 
>> -------
>> #	Save results to new file.
>> #  
>> --------------------------------------------------------------------- 
>> -------
>>
>> nad83 = osr.SpatialReference()
>> nad83.SetFromUserInput('NAD83')
>>
>> shp_driver = ogr.GetDriverByName( 'ESRI Shapefile' )
>> shp_driver.DeleteDataSource( dstfile )
>>
>> shp_ds = shp_driver.CreateDataSource( dstfile )
>>
>> shp_layer = shp_ds.CreateLayer( 'out', geom_type = ogr.wkbPolygon,
>>                                 srs = nad83 )
>>
>> out_feat = ogr.Feature(feature_def = shp_layer.GetLayerDefn() )
>> out_feat.SetGeometry( poly )
>>
>> shp_layer.CreateFeature( out_feat )
>>
>> out_feat.Destroy()
>>
>> shp_layer = None
>> shp_ds.Destroy()
>>
>