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 > > <[email protected]>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 > > > > > > <[email protected]>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 > > > > > > > > > > <[email protected]>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 _______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
