[Gdal-dev] python - get projection from a file

[Brent Pedersen]

hi, thanks to all who replied, i'm familiar with the command-line tools from
gdal and ogr, just getting used to this API.  plus, i'm just switching from
perl to python to take advantage of the greater number of GIS libs.
for future reference (and in case anyone has pointers on how to improve this
code) here's a script that will read a point file translate it into long/lat
projection and print out the new points and one Attirbute column ('ATTRIB'
in this case).

#!/usr/bin/python -w

import gdal
import osr
import ogr

fd = ogr.Open('treetrail.shp');
lyr = fd.GetLayer(0);
osrs = lyr.GetSpatialRef();

latlong = osr.SpatialReference();
latlong.ImportFromProj4('+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs');
#print latlong.ExportToProj4();

transform = osr.CoordinateTransformation(osrs,latlong)

nlayers = lyr.GetFeatureCount();
i = 0

field_of_interest = 'ATTRIB';

for i in range(0,nlayers - 1):
    feat = lyr.GetNextFeature()

    info = feat.GetField(field_of_interest)
    fgeom = feat.GetGeometryRef()

    x = fgeom.GetX()
    y = fgeom.GetY()

    # must be a way to use TransformPoints() ...
    [tx,ty,z] = transform.TransformPoint(x,y);
    print "%f,%f,%s" % (tx,ty,info)




On 3/6/06, Matthew Perry <perrygeo at gmail.com> wrote:
>
> Brent,
>
> Once you open the dataset, you have to get the layer then the spatial
> reference for the layer. From there you can export the spatial
> reference to a variety of formats as needed (the "exportTo*" methods
> at http://gdal.maptools.org/ogr/classOGRSpatialReference.html)
>
> >>> import ogr
> >>> ds = ogr.Open('ca_counties.shp')
> >>> layer = ds.GetLayer()
> >>> sr = layer.GetSpatialRef()
> >>> sr
> <osr.SpatialReference instance at 0xb7b5510c>
> >>> sr.ExportToProj4()
> '+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997
> +units=m +no_defs '
> >>> sr.ExportToWkt()
> 'PROJCS["US National Atlas Equal Area",GEOGCS["Unspecified datum based
> upon the GRS 1980 Authalic
>
> Sphere",DATUM["Not_specified_based_on_GRS_1980_Authalic_Sphere",SPHEROID["GRS_1980_Authalic_Sphere",6370997,0]],PRIMEM["Greenwich",0],UNIT["Degree",
> 0.017453292519943295
> ]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",45],PARAMETER["longitude_of_center",-100],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]'
>
> - matt
>
> On 3/6/06, Brent Pedersen <bpederse at gmail.com> wrote:
> > hi, i'd like to be able to have a user upload a file (.shapefile for
> now). i
> > want to then use python/gdal to get the projection information
> associated
> > with that file. where do i start?
> >
> > i have ogr.Open(filename) and the simple stuff like that, but how do i
> get
> > the projection? i am having trouble seeing how ogr fits in with stuff
> like
> > osr.SpatialReference.ImportFromESRI ...
> >
> > thanks for any pointers,
> > -brent
> >
> > _______________________________________________
> > Gdal-dev mailing list
> > Gdal-dev at lists.maptools.org
> > http://lists.maptools.org/mailman/listinfo/gdal-dev
> >
> >
>
>
> --
> Matt Perry
> perrygeo at gmail.com
> http://www.perrygeo.net
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20060306/8386d5b7/attachment.html