[gdal-dev] Re: Transform not working

Jamie Lahowetz deadpickle at gmail.com
Mon Mar 1 20:05:03 EST 2010 

Thanks for the response. I'm sorry I'm really new to this and am not sure
where to transform the line geometry, I know it should be before
layer.CreateFeature(feature) but after I add the points to the geometry. Is
it before or after feature.SetGeometry(line)? And is the right way to call
the transform
target = osr.SpatialReference()
target.ImportFromProj4('+proj=lcc +lat_1=33.000000 +lat_2=45.000000
+lat_0=39.000000 +lon_0=-96.000000 +x_0=0.0 +y_0=0.0 +datum=NAD83')
coordtrans = osr.CoordinateTransformation(source, target)
line.Transform(coordtrans)?



import os
try:
    from osgeo import ogr
    from osgeo import osr

except:
    import ogr
    import osr

#~ start driver
driver = ogr.GetDriverByName('ESRI Shapefile')
#~ set workspace
ws = 'C:\\Users\\deadpickle\\Desktop\\case study\\verify\\'
os.chdir(ws)

#~ check if it exists
if os.path.exists('test.shp'):
    driver.DeleteDataSource('test.shp')

#~ create a new shp
ds = driver.CreateDataSource('test.shp')

#~ create spatref for new shp
source = osr.SpatialReference()
source.SetFromUserInput('WGS84')
#~ source.SetFromUserInput('+proj=lcc +lat_1=33.000000 +lat_2=45.000000
+lat_0=39.000000 +lon_0=-96.000000 +x_0=0.0 +y_0=0.0 +datum=NAD83')
#~ create a new data layer
layer = ds.CreateLayer('test',geom_type=ogr.wkbLineString,srs=source)

#~ add an id field
fieldDefn = ogr.FieldDefn('id', ogr.OFTInteger)
layer.CreateField(fieldDefn)

#~ create a geometry
line = ogr.Geometry(ogr.wkbLineString)

#~ get the def of the layer, scheme
featureDefn = layer.GetLayerDefn()

# create a new feature
feature = ogr.Feature(featureDefn)

#~ Read 2006track file
file = open("2006track.csv","r")
file.readline()
lnlist = file.readlines()
for lines in lnlist:
    value = lines.split(",")
    if float(value[1]) == 2:
        print value
        line.AddPoint(float(value[3]),float(value[2]))
        #~ set the field 'ID'
        feature.SetField('id', float(value[1]))

#~ feature = layer.GetFeature(0)
#~ geom = feature.GetGeometryRef()

#~ change projection to lcc
target = osr.SpatialReference()
target.ImportFromProj4('+proj=lcc +lat_1=33.000000 +lat_2=45.000000
+lat_0=39.000000 +lon_0=-96.000000 +x_0=0.0 +y_0=0.0 +datum=NAD83')
coordtrans = osr.CoordinateTransformation(source, target)
line.Transform(coordtrans)

#~ set the geometry for this shp
feature.SetGeometry(line)

# add the feature to the output layer
layer.CreateFeature(feature)

# destroy the geometry and feature and close the data source
ds.Destroy()
line.Destroy()
feature.Destroy()
file.close()




On Mon, Mar 1, 2010 at 6:49 PM, Chris Garrard
<garrard.chris+gis at gmail.com<garrard.chris%2Bgis at gmail.com>
> wrote:

> Jamie, Try transforming the line geometry before adding it to the feature.
> That should do it.
>
> chris
>
>
> On Mon, Mar 1, 2010 at 4:24 PM, deadpickle <deadpickle at gmail.com> wrote:
>
>> Not much of a response going on here. In order to change the shapefile
>> projection from geographic to projected do I need to reopen the file,
>> get each feature, and transform the points individually? If so what
>> calls do I need to invoke?
>>
>> On Mar 1, 12:54 pm, Jamie Lahowetz <deadpic... at gmail.com> wrote:
>> > I'm new to GDAL and am usingthe python bindings. I want to transform a
>> newly
>> > created shapefile from WGS84 to a Lambert Conical Conformal so that I
>> can
>> > calculate the length of the contained polylines. The script exits with
>> no
>> > errors but when I look at the shapefile in arcCatalog I see that the
>> > projection is still WGS84. Any ideas? Also, is there a way to calculate
>> the
>> > length directly from GDAL? Thanks.
>> >
>> > import os
>> > try:
>> >     from osgeo import ogr
>> >     from osgeo import osr
>> >
>> > except:
>> >     import ogr
>> >     import osr
>> >
>> > #~ start driver
>> > driver = ogr.GetDriverByName('ESRI Shapefile')
>> > #~ set workspace
>> > ws = 'C:\\Users\\deadpickle\\
>> > Desktop\\case study\\verify\\'
>> > os.chdir(ws)
>> >
>> > #~ check if it exists
>> > if os.path.exists('test.shp'):
>> >     driver.DeleteDataSource('test.shp')
>> >
>> > #~ create a new shp
>> > ds = driver.CreateDataSource('test.shp')
>> >
>> > #~ create spatref for new shp
>> > source = osr.SpatialReference()
>> > source.SetFromUserInput('WGS84')
>> >
>> > #~ create a new data layer
>> > layer = ds.CreateLayer('test',geom_type=ogr.wkbLineString,srs=source)
>> >
>> > #~ add an id field
>> > fieldDefn = ogr.FieldDefn('id', ogr.OFTInteger)
>> > layer.CreateField(fieldDefn)
>> >
>> > #~ create a geometry
>> > line = ogr.Geometry(ogr.wkbLineString)
>> >
>> > #~ get the def of the layer, scheme
>> > featureDefn = layer.GetLayerDefn()
>> >
>> > # create a new feature
>> > feature = ogr.Feature(featureDefn)
>> >
>> > #~ Read 2006track file
>> > file = open("2006track.csv","r")
>> > file.readline()
>> > lnlist = file.readlines()
>> > for lines in lnlist:
>> >     value = lines.split(",")
>> >     if float(value[1]) == 2:
>> >         print value
>> >         line.AddPoint(float(value[3]),float(value[2]))
>> >         #~ set the field 'ID'
>> >         feature.SetField('id', float(value[1]))
>> >
>> > #~ set the geometry for this shp
>> > feature.SetGeometry(line)
>> >
>> > # add the feature to the output layer
>> > layer.CreateFeature(feature)
>> >
>> > #~ feature = layer.GetFeature(0)
>> > geom = feature.GetGeometryRef()
>> >
>> > #~ change projection to lcc
>> > target = osr.SpatialReference()
>> > target.ImportFromProj4('+proj=lcc +lat_1= 33.000000 +lat_2=45.000000
>> > +lat_0=39.000000 +lon_0=-96.000000 +x_0=0.0 +y_0=0.0')
>> > coordtrans = osr.CoordinateTransformation(source, target)
>> >
>> > geom.Transform(coordtrans)
>> >
>> > # destroy the geometry and feature and close the data source
>> >
>> > ds.Destroy()
>> > line.Destroy()
>> > feature.Destroy()
>> > file.close()
>> >
>> > _______________________________________________
>> > gdal-dev mailing list
>> > gdal-... at lists.osgeo.orghttp://
>> lists.osgeo.org/mailman/listinfo/gdal-dev
>> _______________________________________________
>> gdal-dev mailing list
>> gdal-dev at lists.osgeo.org
>> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>>
>
>


-- 
Jamie Ryan Lahowetz
University of Nebraska - Lincoln
Graduate Student - Geosciences
402.304.0766
jlahowetz at gmail.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20100301/badb0244/attachment.html

More information about the gdal-dev mailing list