[gdal-dev] writing XYZ file

[Even Rouault]

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