[gdal-dev] writing XYZ file

[Even Rouault]

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