[RNFdev] Not a pretty picture

[Dan Putler]

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()
>