[gdal-dev] writing XYZ file

[denis cohen]

Salut,

Almost there, it looks like the data is not inverted, it still starts from
the upper left corner instead of lower left.

Thanks

Denis



On Sat, Nov 9, 2013 at 8:46 PM, Even Rouault
<even.rouault at mines-paris.org>wrote:

> Le samedi 09 novembre 2013 20:39:50, denis cohen a écrit :
> > Hi,
> >
> > Thanks for the script. It works fine except for the actual data
> (elevation
> > in my case, this is a DEM).
> > It does not write the correct values.
>
> Ah right, this version only worked for Byte data type. The below version
> should work for any GDAL data type
>
> ----- Begin of script invert.py ------
>
> # -*- coding: utf-8 -*-
> import sys
> from osgeo import gdal
>
> src_ds = gdal.Open(sys.argv[1])
> xsize = src_ds.RasterXSize
> ysize = src_ds.RasterYSize
> bandcount = src_ds.RasterCount
> datatype = src_ds.GetRasterBand(1).DataType
> out_drv = gdal.GetDriverByName('GTiff')
> out_ds = out_drv.Create(sys.argv[2], xsize, ysize, bandcount, datatype)
> out_ds.SetProjection(src_ds.GetProjectionRef())
> src_gt = src_ds.GetGeoTransform()
> out_gt = [ src_gt[i] for i in range(6) ]
> out_gt[3] = src_gt[3] + ysize * src_gt[5]
> out_gt[5] = -src_gt[5]
> out_ds.SetGeoTransform(out_gt)
> for j in range(src_ds.RasterYSize):
>     data = src_ds.ReadRaster(0, ysize - 1 - j, xsize, 1, buf_type =
> datatype)
>     out_ds.WriteRaster(0, j, xsize, 1, data, buf_type = datatype)
> out_ds = None
> src_ds = None
>
> ----- End of script ------
>
>
>
> >
> > Merci
> >
> > Denis
> >
> >
> >
> > On Sat, Nov 9, 2013 at 8:02 PM, Even Rouault
> >
> > <even.rouault at mines-paris.org>wrote:
> > > Le samedi 09 novembre 2013 10:46:20, denis cohen a écrit :
> > > > Hello,
> > > >
> > > > I am a new user of gdal.
> > > > I used gdal_translate to convert a netcdf file to xyz:
> > > > gdal_translate -of XYZ file.nc file.xyz
> > > >
> > > > However, the xyz file starts with the upperleft coordinate.
> > > > Is it possible to have it write the xyz file starting with the lower
> > > > left coordinate?
> > >
> > > Not directly although it wouldn't be complicated to add an option to do
> > > that.
> > >
> > > You can try the following GDAL Python script that will create an
> > > intermediary
> > > TIFF image that has both inverted georeferencing and imagery in a
> > > consistant
> > > way.
> > >
> > > python invert.py file.nc file_reserved.tif
> > > gdal_translate -of XYZ file_reversed.tif file.xyz
> > >
> > >
> > > ----- Begin of script invert.py ------
> > >
> > > # -*- coding: utf-8 -*-
> > > import sys
> > > from osgeo import gdal
> > >
> > > src_ds = gdal.Open(sys.argv[1])
> > > xsize = src_ds.RasterXSize
> > > ysize = src_ds.RasterYSize
> > > bandcount = src_ds.RasterCount
> > > out_drv = gdal.GetDriverByName('GTiff')
> > > out_ds = out_drv.Create(sys.argv[2], xsize, ysize, bandcount)
> > > out_ds.SetProjection(src_ds.GetProjectionRef())
> > > src_gt = src_ds.GetGeoTransform()
> > > out_gt = [ src_gt[i] for i in range(6) ]
> > > out_gt[3] = src_gt[3] + ysize * src_gt[5]
> > > out_gt[5] = -src_gt[5]
> > > out_ds.SetGeoTransform(out_gt)
> > >
> > > for j in range(src_ds.RasterYSize):
> > >     data = src_ds.ReadRaster(0, ysize - 1 - j, xsize, 1)
> > >     out_ds.WriteRaster(0, j, xsize, 1, data)
> > >
> > > out_ds = None
> > > src_ds = None
> > >
> > > ----- End of script ------
> > >
> > > Even
> > >
> > > --
> > > Geospatial professional services
> > > http://even.rouault.free.fr/services.html
>
> --
> Geospatial professional services
> http://even.rouault.free.fr/services.html
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20131109/c5f59c11/attachment-0001.html>