[gdal-dev] writing XYZ file

denis cohen dcohen at iastate.edu
Sat Nov 9 11:59:11 PST 2013


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>


More information about the gdal-dev mailing list