[gdal-dev] Re: Transform not working

[Jamie Lahowetz]

I figured it out... with all of your help. I had the coordtrans =
osr.CoordinateTransformation switched. I wanted to input lat/lon (WGS84) and
transform it to lcc (source) and I also had to set the source prj to lcc. I
think you all said this stuff already, just took me a little while. Thanks
for all the help.


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()
target = osr.SpatialReference()
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')
target.SetFromUserInput('WGS84')
coordtrans = osr.CoordinateTransformation(target, source)

#~ 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 float(value[3]),float(value[2])
        line.AddPoint(float(value[3]),float(value[2]))
        #~ set the field 'ID'
        feature.SetField('id', float(value[1]))

#~ trans from lat/lon input to lcc
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 12:54 PM, Jamie Lahowetz <deadpickle 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()
>



-- 
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/20100302/2d3e9b7e/attachment.html