<div dir="ltr"><div><div>Morning,<br></div>Sorry about this. You are right. It works. My mistake.<br>Thanks a lot for your help.<br></div>Denis<br><br></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Sat, Nov 9, 2013 at 11:01 PM, Even Rouault <span dir="ltr"><<a href="mailto:even.rouault@mines-paris.org" target="_blank">even.rouault@mines-paris.org</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Le samedi 09 novembre 2013 21:17:07, denis cohen a écrit :<br>
> Hi,<br>
><br>
<div class="im">> May be I wasn't clear. The x,y coordinates start from the lower left but<br>
> not the elevation (z coord).<br>
<br>
</div>Well, it seems to work as expected with my test dataset. Are you sure about<br>
that ? Perhaps there's something particular with your dataset that I can't see<br>
at my end, although I don't see what could go wrong for you and not for me.<br>
<div class="HOEnZb"><div class="h5"><br>
><br>
> Denis<br>
><br>
><br>
><br>
> On Sat, Nov 9, 2013 at 9:12 PM, Even Rouault<br>
><br>
> <<a href="mailto:even.rouault@mines-paris.org">even.rouault@mines-paris.org</a>>wrote:<br>
> > Le samedi 09 novembre 2013 20:59:11, denis cohen a écrit :<br>
> > > Salut,<br>
> > ><br>
> > > Almost there, it looks like the data is not inverted, it still starts<br>
> ><br>
> > from<br>
> ><br>
> > > the upper left corner instead of lower left.<br>
> ><br>
> > Really ? Not in my testing.<br>
> ><br>
> > $ gdalinfo ../autotest/gdrivers/data/small_world.tif<br>
> > Driver: GTiff/GeoTIFF<br>
> > Files: ../autotest/gdrivers/data/small_world.tif<br>
> > Size is 400, 200<br>
> > Coordinate System is:<br>
> > GEOGCS["WGS 84",<br>
> ><br>
> > DATUM["WGS_1984",<br>
> ><br>
> > SPHEROID["WGS 84",6378137,298.257223563,<br>
> ><br>
> > AUTHORITY["EPSG","7030"]],<br>
> ><br>
> > AUTHORITY["EPSG","6326"]],<br>
> ><br>
> > PRIMEM["Greenwich",0],<br>
> > UNIT["degree",0.0174532925199433],<br>
> > AUTHORITY["EPSG","4326"]]<br>
> ><br>
> > Origin = (-180.000000000000000,90.000000000000000)<br>
> > Pixel Size = (0.900000000000000,-0.900000000000000)<br>
> ><br>
> > Metadata:<br>
> > AREA_OR_POINT=Area<br>
> ><br>
> > Image Structure Metadata:<br>
> > INTERLEAVE=BAND<br>
> ><br>
> > Corner Coordinates:<br>
> > Upper Left (-180.0000000, 90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"N)<br>
> > Lower Left (-180.0000000, -90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"S)<br>
> > Upper Right ( 180.0000000, 90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"N)<br>
> > Lower Right ( 180.0000000, -90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"S)<br>
> > Center ( 0.0000000, 0.0000000) ( 0d 0' 0.01"E, 0d 0' 0.01"N)<br>
> ><br>
> > $ python create_upsidedown_tif.py<br>
> > ../autotest/gdrivers/data/small_world.tif reverse.tif<br>
> ><br>
> > $ gdal_translate reverse.tif out.xyz -of xyz<br>
> ><br>
> > $ head out.xyz<br>
> > -179.550000000000011 -89.5499999999999972 215<br>
> > [...]<br>
> ><br>
> > > Thanks<br>
> > ><br>
> > > Denis<br>
> > ><br>
> > ><br>
> > ><br>
> > > On Sat, Nov 9, 2013 at 8:46 PM, Even Rouault<br>
> > ><br>
> > > <<a href="mailto:even.rouault@mines-paris.org">even.rouault@mines-paris.org</a>>wrote:<br>
> > > > Le samedi 09 novembre 2013 20:39:50, denis cohen a écrit :<br>
> > > > > Hi,<br>
> > > > ><br>
> > > > > Thanks for the script. It works fine except for the actual data<br>
> > > ><br>
> > > > (elevation<br>
> > > ><br>
> > > > > in my case, this is a DEM).<br>
> > > > > It does not write the correct values.<br>
> > > ><br>
> > > > Ah right, this version only worked for Byte data type. The below<br>
> ><br>
> > version<br>
> ><br>
> > > > should work for any GDAL data type<br>
> > > ><br>
> > > > ----- Begin of script invert.py ------<br>
> > > ><br>
> > > > # -*- coding: utf-8 -*-<br>
> > > > import sys<br>
> > > > from osgeo import gdal<br>
> > > ><br>
> > > > src_ds = gdal.Open(sys.argv[1])<br>
> > > > xsize = src_ds.RasterXSize<br>
> > > > ysize = src_ds.RasterYSize<br>
> > > > bandcount = src_ds.RasterCount<br>
> > > > datatype = src_ds.GetRasterBand(1).DataType<br>
> > > > out_drv = gdal.GetDriverByName('GTiff')<br>
> > > > out_ds = out_drv.Create(sys.argv[2], xsize, ysize, bandcount,<br>
> > > > datatype) out_ds.SetProjection(src_ds.GetProjectionRef())<br>
> > > > src_gt = src_ds.GetGeoTransform()<br>
> > > > out_gt = [ src_gt[i] for i in range(6) ]<br>
> > > > out_gt[3] = src_gt[3] + ysize * src_gt[5]<br>
> > > > out_gt[5] = -src_gt[5]<br>
> > > > out_ds.SetGeoTransform(out_gt)<br>
> > > ><br>
> > > > for j in range(src_ds.RasterYSize):<br>
> > > > data = src_ds.ReadRaster(0, ysize - 1 - j, xsize, 1, buf_type =<br>
> > > ><br>
> > > > datatype)<br>
> > > ><br>
> > > > out_ds.WriteRaster(0, j, xsize, 1, data, buf_type = datatype)<br>
> > > ><br>
> > > > out_ds = None<br>
> > > > src_ds = None<br>
> > > ><br>
> > > > ----- End of script ------<br>
> > > ><br>
> > > > > Merci<br>
> > > > ><br>
> > > > > Denis<br>
> > > > ><br>
> > > > ><br>
> > > > ><br>
> > > > > On Sat, Nov 9, 2013 at 8:02 PM, Even Rouault<br>
> > > > ><br>
> > > > > <<a href="mailto:even.rouault@mines-paris.org">even.rouault@mines-paris.org</a>>wrote:<br>
> > > > > > Le samedi 09 novembre 2013 10:46:20, denis cohen a écrit :<br>
> > > > > > > Hello,<br>
> > > > > > ><br>
> > > > > > > I am a new user of gdal.<br>
> > > > > > > I used gdal_translate to convert a netcdf file to xyz:<br>
> > > > > > > gdal_translate -of XYZ <a href="http://file.nc" target="_blank">file.nc</a> file.xyz<br>
> > > > > > ><br>
> > > > > > > However, the xyz file starts with the upperleft coordinate.<br>
> > > > > > > Is it possible to have it write the xyz file starting with the<br>
> > > > > > > lower left coordinate?<br>
> > > > > ><br>
> > > > > > Not directly although it wouldn't be complicated to add an option<br>
> ><br>
> > to<br>
> ><br>
> > > > > > do that.<br>
> > > > > ><br>
> > > > > > You can try the following GDAL Python script that will create an<br>
> > > > > > intermediary<br>
> > > > > > TIFF image that has both inverted georeferencing and imagery in a<br>
> > > > > > consistant<br>
> > > > > > way.<br>
> > > > > ><br>
> > > > > > python invert.py <a href="http://file.nc" target="_blank">file.nc</a> file_reserved.tif<br>
> > > > > > gdal_translate -of XYZ file_reversed.tif file.xyz<br>
> > > > > ><br>
> > > > > ><br>
> > > > > > ----- Begin of script invert.py ------<br>
> > > > > ><br>
> > > > > > # -*- coding: utf-8 -*-<br>
> > > > > > import sys<br>
> > > > > > from osgeo import gdal<br>
> > > > > ><br>
> > > > > > src_ds = gdal.Open(sys.argv[1])<br>
> > > > > > xsize = src_ds.RasterXSize<br>
> > > > > > ysize = src_ds.RasterYSize<br>
> > > > > > bandcount = src_ds.RasterCount<br>
> > > > > > out_drv = gdal.GetDriverByName('GTiff')<br>
> > > > > > out_ds = out_drv.Create(sys.argv[2], xsize, ysize, bandcount)<br>
> > > > > > out_ds.SetProjection(src_ds.GetProjectionRef())<br>
> > > > > > src_gt = src_ds.GetGeoTransform()<br>
> > > > > > out_gt = [ src_gt[i] for i in range(6) ]<br>
> > > > > > out_gt[3] = src_gt[3] + ysize * src_gt[5]<br>
> > > > > > out_gt[5] = -src_gt[5]<br>
> > > > > > out_ds.SetGeoTransform(out_gt)<br>
> > > > > ><br>
> > > > > > for j in range(src_ds.RasterYSize):<br>
> > > > > > data = src_ds.ReadRaster(0, ysize - 1 - j, xsize, 1)<br>
> > > > > > out_ds.WriteRaster(0, j, xsize, 1, data)<br>
> > > > > ><br>
> > > > > > out_ds = None<br>
> > > > > > src_ds = None<br>
> > > > > ><br>
> > > > > > ----- End of script ------<br>
> > > > > ><br>
> > > > > > Even<br>
> > > > > ><br>
> > > > > > --<br>
> > > > > > Geospatial professional services<br>
> > > > > > <a href="http://even.rouault.free.fr/services.html" target="_blank">http://even.rouault.free.fr/services.html</a><br>
> > > ><br>
> > > > --<br>
> > > > Geospatial professional services<br>
> > > > <a href="http://even.rouault.free.fr/services.html" target="_blank">http://even.rouault.free.fr/services.html</a><br>
> ><br>
> > --<br>
> > Geospatial professional services<br>
> > <a href="http://even.rouault.free.fr/services.html" target="_blank">http://even.rouault.free.fr/services.html</a><br>
<br>
--<br>
Geospatial professional services<br>
<a href="http://even.rouault.free.fr/services.html" target="_blank">http://even.rouault.free.fr/services.html</a><br>
</div></div></blockquote></div><br></div>