[Qgis-developer] Help Needed - How to get the projection parameters for a Raster file

[Germán Carrillo]

Hi Bishwarup,

some time ago I needed something like that. I got it with gdal and osr.

Here is the Python snippet:
---------------------------------------------------------------
# *Python terminal*: First, get the transform parameters, you'll use them
after in the code
dataset = gdal.Open("/ruta/archivo_original.tif", gdal.GA_ReadOnly)
parameters = ds.GetGeoTransform()
---------------------------------------------------------------
# *On a separate Python file*
import numpy
from osgeo import gdal
from osgeo import osr

# Cargue del archivo que contiene la matriz
array = numpy.loadtxt("/ruta/archivo_cp1",delimiter=';')

# Formato de salida
out_driver = gdal.GetDriverByName("GTiff")

# Creación del DataSet de salida
outds = out_driver.Create("/ruta/archivo_salida.tif", 1272, 787, 1,
    gdal.GDT_Float32, options=["COMPRESS=PACKBITS","TFW=YES"] )

# *Here you should put the parameters you got before*
outds.SetGeoTransform([552910.96327099996,28.5,0.0,422591.15634099999,0.0,-28.5])

# Definir el sistema de referencia de la imagen de salida
srs = osr.SpatialReference()
srs.ImportFromEPSG( 32618 )
outds.SetProjection( srs.ExportToWkt() )

# Escribir la matriz 'array' en la primera banda de la imagen de salida
outband = outds.GetRasterBand( 1 )
outband.WriteArray( array )
---------------------------------------------------------------

Hope it helps you.


Regards,

Germán


On Tue, Nov 2, 2010 at 9:21 AM, Barry Rowlingson <
b.rowlingson at lancaster.ac.uk> wrote:

> On Tue, Nov 2, 2010 at 1:53 PM, Bishwarup Banerjee
> <bishwarup.banerjee at gmail.com> wrote:
> > Dear All,
> >
> > Can any body help me to get and set projection parameters for raster
> files?
> >
> > I am doing some map calculation using PIL and the data is losing
> projection,
> > so want to re project the raster files.
> >
> > The language I am using is Python on windows platform..
>
>  Assuming this is a qgis question...
>
>  If you get the layer in 'l' then l.crs() gets you a coordinate
> reference system object:
>
> >>> l.crs()
> <qgis.core.QgsCoordinateReferenceSystem object at 0xb3b44ac>
>
> and hence:
>
> >>> l.crs().toProj4()
> PyQt4.QtCore.QString(u'+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
>
>  and to set it, just use setCrs with a QgsCoordinateReferenceSystem object.
>
> Barry
> _______________________________________________
> Qgis-developer mailing list
> Qgis-developer at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/qgis-developer
>



-- 
-----------
  |\__
(:>__)(
  |/

Soluciones Geoinformáticas Libres
http://geotux.tuxfamily.org/
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/qgis-developer/attachments/20101105/834671d9/attachment.html