[gdal-dev] writing XYZ file

denis cohen dcohen at iastate.edu
Sat Nov 9 20:27:40 PST 2013 

Morning,
Sorry about this. You are right. It works. My mistake.
Thanks a lot for your help.
Denis



On Sat, Nov 9, 2013 at 11:01 PM, Even Rouault
<even.rouault at mines-paris.org>wrote:

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

More information about the gdal-dev mailing list