[Gdal-dev] Append two shapefiles

[Dylan Beaudette]

Hi Alex,

You might want to talk to Howard Butler about that- he seems to be doing a lot 
of work on Python-GDAL stuff in windows.

Good luck,

Dylan

On Monday 05 December 2005 05:49 pm, Alex Mandel wrote:
> Thanks Dylan -
> I was putting full paths for the layers and using the -f, and for this
> time I also switched my slashes \ to / and it works great now.
>
> The next questions is how can I get just the ogr2ogr component into my
> program. I need to distribute a very thin and lite python program and
> can't expect my users to have FWtools.
>
> Ironically I think I'm one error away from writing a python module with
> gdal to do what I want, which would save me a ton of MB and give me
> control to copy only certain points.
>
> There's something wrong with the geometry line though. boa says that
> feat2.SetGeometry(pnt) has invalid syntax
> Anyone see my mistake right away, I figure it's something obvious that I
> just haven't learned yet. Keep in mind I'm only dealing with points, and
> will never need to handle any other geometry type.
>
> Thanks Everyone - Alex
> ---------------------------------------
>
> import gdal
> import gdal.ogr as ogr
>
> #Set ESRI Shapefile driver just in case
> driver = ogr.GetDriverByName('ESRI Shapefile')
> #Bring in 1st Shapefile, src
> ds1 = driver.Open("C:/Boxes/sandbox/one/one.shp")
> pLayer1 = ds1.GetLayer(0)
> #initialize first feature read
> feat = pLayer1.GetNextFeature()
>
> #Bring in 2nd Shapefile, dest
> ds2 = driver.Open("C:/Boxes/sandbox/two/two.shp")
> pLayer2 = ds2.GetLayer(0)
>
> #Count the numnber of fields
> layerDef2 = pLayer2.GetLayerDefn()
> field_count = layerDef.GetFieldCount()
>
> #Create index of fields to copy
> fld_list = {}
>
> for i in range(layerDef2.GetFieldCount()):
>      field = layerDef.GetFieldDefn(i)
>      fld_list[field.GetName()] = i
>
> while feat is not None:
>
>      multi_geom = feat.GetGeometryRef()
>
>      for i in range(multi_geom.GetGeometryCount()):
>          pnt = multi_geom.GetGeometryRef(i)
>
>          feat2 = ogr.Feature(feature_def=layerDef2)
> #Copy Feature table data
>          for fld_index in range(field_count):
>              fld_name = layerDef.GetFieldDefn(fld_index)
>              #if fld_list.has_key(fld_name.GetName()):
>              feat2.SetField(fld_index, \
> feat.GetField(fld_list.get(fld_name))
>
> #Paste into destination
>          feat2.SetGeometry(pnt)
>          pLayer2.CreateFeature(feat2)
>          feat2.Destroy()
>
> #Clean up and iterate
> 	feat.Destroy()
> 	feat = pLayer1.GetNextFeature()
>
>
> ----------------------------------------
>
> Dylan Beaudette wrote:
>  > hi Alex!
>  >
>  > didn't know you were lurking here!
>  >
>  > welcome the GDAL land --
>  >
>  > lots of good stuff here.
>  >
>  >
>  > here is how i do it:
>  >
>  > ogr2ogr -update -append $merge_file $this_file_to_append -nln
>
> $merge_layer
>
>  > $this_file_to_append_layer
>  >
>  > where the 2nd set of 'layers' are missing the '.shp' suffix.
>  >
>  > cheers,
>  >
>  > Dylan

-- 
Dylan Beaudette
Soils and Biogeochemistry Graduate Group
University of California at Davis
530.754.7341