[gdal-dev] writing XYZ file

[denis cohen]

Hi,

May be I wasn't clear. The x,y coordinates start from the lower left but
not the elevation (z coord).

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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20131109/a9f585f6/attachment.html>