deadpickle 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...@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.
Jamie,
My suggestion in the following script would be to create the file with
the "target" srs instead of WGS84, and to transform the geometry
before even assigning it to the feature that will be written.
A few notes:
* There is no way with OGR of changing the coordinate system of
an existing layer. The only point at which it can be set is
layer creation time.
* The CreateFeature() method on the layer creates the feature in
the datastore (ie. shapefile) but there is no meaningful
connection maintained between the ogr.Feature object and the
feature in the file after this. Alterations to the feature
or it's geometry will have no effect on the file without doing
another layer.SetFeature() call.
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-...@lists.osgeo.orghttp://lists.osgeo.org/mailman/listinfo/gdal-dev
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up | Frank Warmerdam, warmer...@pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush | Geospatial Programmer for Rent
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
