[gdal-dev] writing XYZ file

Even Rouault even.rouault at mines-paris.org
Sat Nov 9 14:01:40 PST 2013 

Le samedi 09 novembre 2013 21:17:07, denis cohen a écrit :
> Hi,
> 
> May be I wasn't clear. The x,y coordinates start from the lower left but
> not the elevation (z coord).

Well, it seems to work as expected with my test dataset. Are you sure about 
that ? Perhaps there's something particular with your dataset that I can't see 
at my end, although I don't see what could go wrong for you and not for me.

> 
> Denis
> 
> 
> 
> On Sat, Nov 9, 2013 at 9:12 PM, Even Rouault
> 
> <even.rouault at mines-paris.org>wrote:
> > Le samedi 09 novembre 2013 20:59:11, denis cohen a écrit :
> > > Salut,
> > > 
> > > Almost there, it looks like the data is not inverted, it still starts
> > 
> > from
> > 
> > > the upper left corner instead of lower left.
> > 
> > Really ? Not in my testing.
> > 
> > $ gdalinfo ../autotest/gdrivers/data/small_world.tif
> > Driver: GTiff/GeoTIFF
> > Files: ../autotest/gdrivers/data/small_world.tif
> > Size is 400, 200
> > Coordinate System is:
> > GEOGCS["WGS 84",
> > 
> >     DATUM["WGS_1984",
> >     
> >         SPHEROID["WGS 84",6378137,298.257223563,
> >         
> >             AUTHORITY["EPSG","7030"]],
> >         
> >         AUTHORITY["EPSG","6326"]],
> >     
> >     PRIMEM["Greenwich",0],
> >     UNIT["degree",0.0174532925199433],
> >     AUTHORITY["EPSG","4326"]]
> > 
> > Origin = (-180.000000000000000,90.000000000000000)
> > Pixel Size = (0.900000000000000,-0.900000000000000)
> > 
> > Metadata:
> >   AREA_OR_POINT=Area
> > 
> > Image Structure Metadata:
> >   INTERLEAVE=BAND
> > 
> > Corner Coordinates:
> > Upper Left  (-180.0000000,  90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"N)
> > Lower Left  (-180.0000000, -90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"S)
> > Upper Right ( 180.0000000,  90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"N)
> > Lower Right ( 180.0000000, -90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"S)
> > Center      (   0.0000000,   0.0000000) (  0d 0' 0.01"E,  0d 0' 0.01"N)
> > 
> > $ python create_upsidedown_tif.py
> > ../autotest/gdrivers/data/small_world.tif reverse.tif
> > 
> > $ gdal_translate reverse.tif  out.xyz -of xyz
> > 
> > $ head out.xyz
> > -179.550000000000011 -89.5499999999999972 215
> > [...]
> > 
> > > 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
> > 
> > --
> > Geospatial professional services
> > http://even.rouault.free.fr/services.html

-- 
Geospatial professional services
http://even.rouault.free.fr/services.html

More information about the gdal-dev mailing list