Re: [gdal-dev] convert NETCDF to Geotiff upside down
The following works for the proper orientation but there's a shift. I would consult the documentation as they usually specify conversion params. gdal_translate -of GTiff -a_srs EPSG:4326 -a_ullr -179.9 -72.7 180 72.7 -co COMPRESS=LZW NETCDF:GLOBCOMPLIR_nc.2023080819:data ofn.tif From: gdal-dev on behalf of Dori Sent: Wednesday, August 9, 2023 2:29 PM To: gdal-dev@lists.osgeo.org Subject: [gdal-dev] convert NETCDF to Geotiff upside down Hi, I tried to convert the NETCDF stored in https://noaa-gmgsi-pds.s3.amazonaws.com/index.html#GMGSI_LW/202308/08/19/GLOBCOMPLIR_nc.2023080819 to Geotiff format using following command: gdal_translate 0ol GTiff -a_srs EPSG:4326 -aullr -180 -72.7 179.9 -72.7 -co “COMPRESS=LZW” NETCDF:GLOBCOMPLIR_nc.2023080819:data ofn.tif ofn.tif is upside down. North latitude points are in South and viceversa. How can I solve this upside down output? What causes this? I did notice that the points in the NETCDF file are not equidistant what could be my issue. Thank you for any hint on how to go about getting the correct Geotiff output ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Speeding up gdalwarp process
We've been using Geotrellis (https://geotrellis.io/) for Spark jobs. -marius From: gdal-dev on behalf of Simon Sent: Thursday, January 23, 2020 6:30 PM To: gdal-dev@lists.osgeo.org Subject: [gdal-dev] Speeding up gdalwarp process Hi gdal-devs, I have a question, if there is some way to use gdalwarp in a clustering system (e.g., Sparks) or with GPU to speed up the process? My task involves re-projection and re-sampling of hundreds of high-resolution images. Any ideas to make use of Sparks or GPU is welcomed. Thank you. Simon ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] gdal_calc.py: produce median raster
That's because numpy.median expects one array and keyword arguments, not three arrays. It's unclear to me what gdal_calc.py does under the hood (I admit I did not read it) but the correct way to get the median of the three bands would be: --calc median([A,B,C], 0) The strange part is that--calc median([A,B,C])works as well and generates a different output. When axis is not specified median should return a single value. At this point I suggest you take three 2D arrays, calculate their median as a 2D array, convert the three arrays to rasters, run gdal_calc.py and compare results. -marius Date: Tue, 25 Nov 2014 17:56:46 +0900 From: simlan...@gmail.com To: chris.yes...@ioz.ac.uk CC: gdal-dev@lists.osgeo.org Subject: Re: [gdal-dev] gdal_calc.py: produce median raster Dear Dr Chris Yesson Thanks for your quick response. I got the following error: ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()` Perhaps the expression for calc seems to be revised... On Tue, Nov 25, 2014 at 4:59 PM, Chris Yesson chris.yesson@gmail.com wrote: Dear Simen, You can refer to each band explicitly using the --A_band option.Try something like this gdal_calc.py -A infile --A_band 1 -B infile --B_band 2 -C infile --C_band 3 --outfile outfile --calc median(A,B,C) - ChrisDr Chris Yesson Institute of Zoology, Zoological Society of LondonTel: 020 7449 6267 email: chris.yes...@ioz.ac.uk web: http://www.zsl.org/chrisyesson/ Note: I work Mon-Thurs and do not check email on Fri On 25 November 2014 at 07:44, Simen Langseth simlan...@gmail.com wrote: Dear Respected Authors: I would be grateful if you could share any hints on it. On Tue, Nov 25, 2014 at 4:33 PM, Simen Langseth simlan...@gmail.com wrote: Dear Marius Jigmond: Thanks for your attempt. Unfortunately, it is producing 3 band outfile, rather than 1 band median image. I could not figure out what this code computed, all the resulted 3 bands images are also different. I hope someone can help me. On Tue, Nov 25, 2014 at 12:09 PM, Marius Jigmond mariusjigm...@hotmail.com wrote: Does the following work: gdal_calc.py -A infile --allBands=A --outfile=outfile --calc=median(A) Date: Mon, 24 Nov 2014 19:14:06 +0900 From: simlan...@gmail.com To: gdal-dev@lists.osgeo.org Subject: [gdal-dev] gdal_calc.py: produce median raster I have one GeoTiff file with 3 bands. I want to produce a raster file computing pixel wise median value for the three raster bands. How can I do that? I tried as follows: gdal_calc.py -A infile --allBands --outfile outfile --calc median(A) But could not get any out file. ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev This message has been scanned for viruses by MailControl, a service from BlackSpider Technologies. Click here to report this email as spam. The Zoological Society of London is incorporated by Royal Charter Principal Office England. Company Number RC000749 Registered address: Regent's Park, London, England NW1 4RY Registered Charity in England and Wales no. 208728 _ This e-mail has been sent in confidence to the named addressee(s). If you are not the intended recipient, you must not disclose or distribute it in any form, and you are asked to contact the sender immediately. Views or opinions expressed in this communication may not be those of The Zoological Society of London and, therefore, The Zoological Society of London does not accept legal responsibility for the contents of this message. The recipient(s) must be aware that e-mail is not a secure communication medium and that the contents of this mail may have been altered by a third party in transit. If you have any issues regarding this mail please contact: administra...@zsl.org. ___ This message has been scanned for viruses by MailControl, a service from BlackSpider Technologies. ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] gdal_calc.py: produce median raster
Does the following work: gdal_calc.py -A infile --allBands=A --outfile=outfile --calc=median(A) Date: Mon, 24 Nov 2014 19:14:06 +0900 From: simlan...@gmail.com To: gdal-dev@lists.osgeo.org Subject: [gdal-dev] gdal_calc.py: produce median raster I have one GeoTiff file with 3 bands. I want to produce a raster file computing pixel wise median value for the three raster bands. How can I do that? I tried as follows: gdal_calc.py -A infile --allBands --outfile outfile --calc median(A) But could not get any out file. ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] gdal_translate produces tif with different colors as original
I noticed you're using ArcMap for visualizing the results. Did you check what histogram/gamma stretch gets applied by default when you load the layer? I run into this issue often as results from GDAL operations appear much worse when loaded into ArcMap compared to originals. When I turn off stretching, both histogram/gamma, results are nearly identical to original. -marius To: gdal-dev@lists.osgeo.org From: andre+jo...@nurfuerspam.de Date: Thu, 7 Nov 2013 16:54:32 +0100 Subject: Re: [gdal-dev] gdal_translate produces tif with different colors as original Am 06.11.2013 21:32, schrieb Jukka Rahkonen: You're advised to preprocess your rasters with other tools, such as pct2rgb.py or gdal_translate -expand RGB to operate gdalbuildvrt on RGB rasters instead. I proposed that, but have a look at the color palette of the VRT: http://gis.stackexchange.com/questions/76662/gdal-translate-produces-tif-with-different-colors-as-original Is it a fault of the source or the driver? Greetings, André Joost ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] GDAL 1.10 build fails on FreeBSD
Frank, You'll need a newer PostGIS (2.x+) if you want to compile GDAL 1.10 with PG support. http://trac.osgeo.org/gdal/wiki/frmts_wtkraster.html -marius Date: Tue, 3 Sep 2013 17:02:44 +0200 From: sc...@sarvision.nl To: gdal-dev@lists.osgeo.org Subject: Re: [gdal-dev] GDAL 1.10 build fails on FreeBSD On Tue, 03 Sep 2013 14:55:20 +0200 Frank Broniewski b...@metrico.lu wrote: Hi all, I tried updating Gdal to the latest version on my FreeBSD system, but the build fails: snip c postgisrasterdataset.cpp -fPIC -DPIC -o ../o/.libs/postgisrasterdataset.o postgisrasterdataset.cpp: In function 'void GDALRegister_PostGISRaster()': postgisrasterdataset.cpp:2096: error: 'GDAL_DMD_SUBDATASETS' was not declared in this scope gmake[2]: *** [../o/postgisrasterdataset.lo] Error 1 gmake[2]: Leaving directory `/usr/ports/graphics/gdal/work/gdal-1.10.1/frmts/postgisraster' gmake[1]: *** [postgisraster-install-obj] Error 2 gmake[1]: Leaving directory `/usr/ports/graphics/gdal/work/gdal-1.10.1/frmts' gmake: *** [frmts-target] Error 2 *** [do-build] Error code 1 /snip PostGIS is POSTGIS=1.5.3 GEOS=3.3.8-CAPI-1.7.8 PROJ=Rel. 4.8.0, 6 March 2012 LIBXML=2.7.8 USE_STATS (procs from 1.5 r5385 need upgrade) so apparently no raster support ... Is this a GDAL issue or shall I report this to the FreeBSD crew? Or is my PostGIS too old ??? Frank, just fyi: I recently upgraded gdal on my freebsd system to 1.10, and all went fine. My postgis is significantly newer than yours, though (2.1 I think, or at least 2.0), so that might be the culprit indeed. Best, Vincent. Frank ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Bathymetry data for US reservoirs.
For Texas reservoirs you may contact Jason Kemp of the Texas Water Development Board (see http://www.twdb.texas.gov/surfacewater/surveys/index.asp for contact info). He may have some data or may be able to point you in the right direction. -marius Date: Mon, 6 May 2013 20:58:27 -0700 From: manoj.bi...@gmail.com To: gdal-dev@lists.osgeo.org Subject: [gdal-dev] Bathymetry data for US reservoirs. Hi, I have been using gdal since past few months (without any issue). I am running into a non gdal related issue pertaining to getting depth/bathy data for reservoirs that I am interested in. I am looking for bathymetry data for some USACE reservoirs. I am told this data used to be available in the USACE GIS data clearing house website. But this website has been down for quite some time now. Does anyone know if I can find the same data elsewhere? The following is the link to USACE GIS clearing house but it has been down for quite a while now. http://corpsgeo1.usace.army.mil/ I have wasted quite a bit of time searching for this (incl. contacting USACE unsuccessfully both via phone and email). I would really appreciate any pointers regarding this. Thanks, --MB. ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Question about using crop_to_cutline
At the time I created the ticket I had created a workaround that uses gdal_translate to clip a rectangular window and gdal_rasterize to burn an inverse image of the vector. Interested folks can email me for the script (does the list take attachments?). The script makes some assumptions about your data type (e.g. 0 is assigned as nodata) but you can change those assumptions to suit your needs. I plan to make the feature request I mentioned in #3947. -marius To: gdal-dev@lists.osgeo.org From: jukka.rahko...@mmmtike.fi Date: Fri, 26 Oct 2012 09:05:06 + Subject: Re: [gdal-dev] Question about using crop_to_cutline Even Rouault even.rouault at mines-paris.org writes: This is the very same topic that is discussed in http://trac.osgeo.org/gdal/ticket/3947 . There's no solution to your problem, but some background discussion that explains the current behaviour. This is close to a problem I had once with creating mosaics from individually warped orthophotos. Accurately calculated extents make pixels to slide a bit and individually warped images do not share any common canvas for their pixels. I sketched a plan for forcing the warped image to use a common canvas and found a python guy to make a program. I have been satisfied with the result. I believe that the same solution will work for you. You must widen the -te parameters that is calculated by crop_to_cutline to each direction so that they match exactly with some pixel row and line of the original image. See figures 3, 4, and 5 in http://www.scangis.org/scangis2007/papers/r3_rahkonen.pdf. The python code is there too. Without understanding anything about programming I suppose that the job is done with this: minmax = get_minmax(tm32_coords) minmax_wider = [ (int(math.floor(minmax[0][0])), int(math.floor(minmax[0][1]))), (int(math.ceil (minmax[1][0])), int(math.ceil (minmax[1][1]))), ] log(Extents (min,max): + str(minmax), outfilename) log(Widened extents (min,max): + str(minmax_wider), outfilename) -Jukka Rahkonen- ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Upgraded to 1.9 - Python Bindings Error
On 05/04/2012 08:16 PM, Jay L. wrote:
List,
Just upgraded to 1.9.0 on my windows 7 machine and I can no longer
access gdal via the python bindings. Command line (ogr2ogr, gdalinfo,
etc. are functioning as expected). Specifically, running 'from osgeo
import gdal', returns:
File C:\Python26\ArcGIS10.0\lib\site-packages\osgeo\__init.py__,
line 17, in
swig_import_helper
_mod = imp.load_module('_gdal', fp, pathname, descritpion)
ImportError: DLL load failed: The specified procedure could not be found.
GDAL 1.8 was working wonderfully with ArcGIS. No environmental
variables were altered between installations. Path still contains
C:\Program Files (x86)\GDAL and GDAL_DATA is stills set to C:\Program
Files (x86)\GDAL\gdal-data.
I did try adding the PYTHONPATH variable and setting that to
C:\Python26\ArcGIS10.0, but no change in the error.
GDAL1.9 and the python bindings were downloaded from
http://www.gisinternals.com/sdk/ (the MSVC2010-32bit stable version).
Python 2.6 that comes with ArcGIS 10 was compiled with MSVC2008. That
could be the reason why the bindings aren't working. Try uninstalling
and then installing the matching version.
-marius
Oddly, the gdal python utilities are functioning. That is
gdal2tiles.py for example.
Any suggestions appreciated. I assume a PATH error that I am missing,
but am stumped. The mailing list archives were sparse on this error...
Thanks,
Jay
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] gdal/ogr version (was NotImplementedError: Wrong number...)
On 09/01/2011 08:23 PM, Ole Nielsen wrote: Many thanks for the reply, but I need the Python bindings as we are using ogr deeply embedded in an application and creating of new layers using the SetField command (amongst others) is critical. Hence my worry when we get the error NotImplementedError: Wrong number of arguments for overloaded function 'Feature_SetField'. I am sure it is something simple that changed in more recent versions - such as an extra required argument perhaps - and I'd be grateful if someone who knows the codebase could tell me what this error means and what I must do to make the call succeed. Many thanks Ole Nielsen On Fri, Sep 2, 2011 at 8:15 AM, Eli Adam ea...@co.lincoln.or.us mailto:ea...@co.lincoln.or.us wrote: Ole, Is there a way of telling which one I am using (both when running commandline org2org and when importing the ogr module in python)? ogr2ogr --version will work on the command line for at least the major version number. I'm not sure about in python. This should give you the bindings version picked up by Python: from osgeo import gdal gdal.__version__ You can also check the PYTHONPATH with import sys sys.path to see which location gets picked up first when you import osgeo modules likely locations you'll find the modules: /usr/(local)/lib/python2.x/dist-packages -marius HTH, Eli Cheers and thanks Ole ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] gdal/ogr version (was NotImplementedError: Wrong number...)
From personal experience (and others') you get some odd error messages when you start mixing different versions of libgdal. I would expect the error to go away if you stick to one version of libgdal. Ubuntu-GIS should have all the packages you need without having to mix GDAL versions. -marius On 09/01/2011 09:08 PM, Ole Nielsen wrote: good idea, thanks. However, it appears that qgis 1.7 is the one that depends on gdal 1.6. I am confused. (riab_env)nielso@shakti:~/dev$ sudo dpkg --purge libgdal1-1.6.0 dpkg: dependency problems prevent removal of libgdal1-1.6.0: libqgis1.7.0 depends on libgdal1-1.6.0. qgis-providers depends on libgdal1-1.6.0. qgis depends on libgdal1-1.6.0. python-qgis depends on libgdal1-1.6.0. In any case, we should be able to work with either, so I am hoping someone will shed some light on what the error message really means. Many thanks Ole On Fri, Sep 2, 2011 at 8:28 AM, Ariel Nunez ingenieroar...@gmail.com mailto:ingenieroar...@gmail.com wrote: Ole, why don't you uninstall libgdal1-1.6.0 and the grass one and run ldconfig to make sure everything is using 1.8? Ariel. On Thu, Sep 1, 2011 at 9:23 PM, Ole Nielsen ole.moller.niel...@gmail.com mailto:ole.moller.niel...@gmail.com wrote: Many thanks for the reply, but I need the Python bindings as we are using ogr deeply embedded in an application and creating of new layers using the SetField command (amongst others) is critical. Hence my worry when we get the error NotImplementedError: Wrong number of arguments for overloaded function 'Feature_SetField'. I am sure it is something simple that changed in more recent versions - such as an extra required argument perhaps - and I'd be grateful if someone who knows the codebase could tell me what this error means and what I must do to make the call succeed. Many thanks Ole Nielsen On Fri, Sep 2, 2011 at 8:15 AM, Eli Adam ea...@co.lincoln.or.us mailto:ea...@co.lincoln.or.us wrote: Ole, Is there a way of telling which one I am using (both when running commandline org2org and when importing the ogr module in python)? ogr2ogr --version will work on the command line for at least the major version number. I'm not sure about in python. HTH, Eli Cheers and thanks Ole ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] build error from svn
David, Revision 22725 builds just fine for me. I normally do a make clean before every build. Maybe something went wrong in your build system. -marius On Fri, 2011-07-15 at 15:15 -0400, David Burken wrote: Hi Frank, all, Just fyi, getting build error from svn. Looks like cpl_serv.h is gone now. See below for minor patch. Take care, Dave $ svn update At revision 22725. gt_wkt_srs.cpp:33:22: fatal error: cpl_serv.h: No such file or directory compilation terminated. make[1]: *** [../o/gt_wkt_srs.lo] Error 1 make[1]: Leaving directory `/home/work/osgeo/gdal/frmts/gtiff' make: *** [gtiff-install-obj] Error 2 Compilation exited abnormally with code 2 at Fri Jul 15 15:06:55 Patch: svn diff frmts/gtiff/gt_wkt_srs.cpp Index: frmts/gtiff/gt_wkt_srs.cpp === --- frmts/gtiff/gt_wkt_srs.cpp(revision 22725) +++ frmts/gtiff/gt_wkt_srs.cpp(working copy) @@ -30,7 +30,9 @@ * DEALINGS IN THE SOFTWARE. / -#include cpl_serv.h +// #include cpl_serv.h + #include cpl_csv.h + #include geo_tiffp.h #define CPL_ERROR_H_INCLUDED ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
RE: [gdal-dev] OGR Geometry methods
Indeed, SetAttributeFilter returns the right feature. Thanks Luke.
FYI the ExecuteSQL bug submission is http://trac.osgeo.org/gdal/ticket/4156
-marius
Date: Tue, 12 Jul 2011 12:31:03 -0400
Subject: Re: [gdal-dev] OGR Geometry methods
From: luke.peter...@gmail.com
To: mariusjigm...@hotmail.com
CC: gdal-dev@lists.osgeo.org
On Tue, Jul 12, 2011 at 11:56 AM, Marius Jigmond mariusjigm...@hotmail.com
wrote:
Strangely enough, ExecuteSQL works (it returns an osgeo.ogr.layer object with
one feature) with
'select FID from judeteEPSG3844 where DENJUD = CLUJ'
but later a GetField('FID') request returns Illegal field requested in
GetField()
-marius
Another way to do this is to use the AttributeFilter. Run this by itself and
see if you get the expected result:
from osgeo import ogrfileName = '/data/romania/judeteEPSG3844.shp'shapeDS =
ogr.Open(fileName)shapeLayer =
shapeDS.GetLayerByIndex(0)shapeLayer.SetAttributeFilter('DENJUD=CLUJ')
targetFeature = shapeLayer.GetNextFeature()print
targetFeature.GetField('DENJUD') # should print the DENJUD of the feature
fitting the attribute filter.
Date: Tue, 12 Jul 2011 09:13:00 -0400
Subject: Re: [gdal-dev] OGR Geometry methods
From: luke.peter...@gmail.com
To: mariusjigm...@hotmail.com
-
Luke Peterson
- sent via mobile device -
On Jul 11, 2011 11:00 PM, Marius Jigmond mariusjigm...@hotmail.com wrote:
After some more investigation that is likely NOT the issue. I have an
ExecuteSQL statement which selects a certain polygon based on an attribute
value. Unfortunately it seems to return the wrong feature. The feature I
query for is unique so a duplicate is out of the question. Here's the code:
#!/usr/bin/python
from osgeo import ogr
import sys, math, time, os
aquifer = '/data/romania/judeteEPSG3844.shp'
aqDS = ogr.Open(aquifer)
sql = 'select * from judeteEPSG3844 where DENJUD = CLUJ'
aqLayer = aqDS.ExecuteSQL(sql)
feat = aqLayer.GetFeature(0)
print aqLayer.GetFeatureCount()
print feat.GetField('DENJUD')
aqDS.ReleaseResultSet(aqLayer)
The result:
marius@mobi:~/temp$ run_sql.py
1
TELEORMAN
marius@mobi:~/temp$
This explains why none of my centroids were remotely close to the polygon.
Is there something wrong with my query or should I file a bug?
Can you select the feature by its FID instead? (apologies, I'm banging this out
on my phone ... probably won't run but you should get the idea):
aqLayer = aqDS.GetLayerByIndex(0) #gets the whole shapefile
aqFID = aqLayer.GetFIDColumn() # I forget the method name here... hope this is
right but may have to pull in the layer def
sql = 'select %s from judeteEPSG3844 where DENJUD = CLUJ' % aqFID #
specifically asking for just the FID field in the resultset
sqlResults = aqDS.ExecuteSQL(sql) # sql results same as before
featFID = sqlResults.GetFeature(0).GetFieldAsInt(aqFID) #this should be the
actual FID of the results (problems would occur if your sql returned more than
one hit)
feat = aqLayer.GetFeature(featFID) # using the fid to grab the feature from the
0th layer in the shapefile
This may be a long way around ...
-marius
On Mon, 2011-07-11 at 20:46 -0500, Marius Jigmond wrote:
I suppose a piece of code speaks louder :):
xsect = False
for i in range(gridLayer.GetFeatureCount()):
feat = gridLayer.GetFeature(i)
geom = feat.GetGeometryRef()
point = geom.Centroid()
for j in range(aqLayer.GetFeatureCount()):
aqfeat = aqLayer.GetFeature(j)
aqgeom = aqfeat.GetGeometryRef()
if point.Intersect(aqgeom):
xsect = True
print 'ibound = 1'
break
if xsect:
feat.SetField('IBOUND', 1)
gridLayer.SetFeature(feat)
xsect = False
-marius
On Mon, 2011-07-11 at 20:34 -0500, Marius Jigmond wrote:
Hi everyone,
I am trying to test whether centroids of polygons lie/intersect within
another polygon. I have tried Intersect, Within, and Contains but they
always return false. Should these methods work for my intended purpose or
do I need to implement a point in polygon function? Thanks.
-marius
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] OGR Geometry methods
Hi everyone, I am trying to test whether centroids of polygons lie/intersect within another polygon. I have tried Intersect, Within, and Contains but they always return false. Should these methods work for my intended purpose or do I need to implement a point in polygon function? Thanks. -marius ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] OGR Geometry methods
I suppose a piece of code speaks louder :):
xsect = False
for i in range(gridLayer.GetFeatureCount()):
feat = gridLayer.GetFeature(i)
geom = feat.GetGeometryRef()
point = geom.Centroid()
for j in range(aqLayer.GetFeatureCount()):
aqfeat = aqLayer.GetFeature(j)
aqgeom = aqfeat.GetGeometryRef()
if point.Intersect(aqgeom):
xsect = True
print 'ibound = 1'
break
if xsect:
feat.SetField('IBOUND', 1)
gridLayer.SetFeature(feat)
xsect = False
-marius
On Mon, 2011-07-11 at 20:34 -0500, Marius Jigmond wrote:
Hi everyone,
I am trying to test whether centroids of polygons lie/intersect within
another polygon. I have tried Intersect, Within, and Contains but they
always return false. Should these methods work for my intended purpose
or do I need to implement a point in polygon function? Thanks.
-marius
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] OGR Geometry methods
After some more investigation that is likely NOT the issue. I have an
ExecuteSQL statement which selects a certain polygon based on an
attribute value. Unfortunately it seems to return the wrong feature. The
feature I query for is unique so a duplicate is out of the question.
Here's the code:
#!/usr/bin/python
from osgeo import ogr
import sys, math, time, os
aquifer = '/data/romania/judeteEPSG3844.shp'
aqDS = ogr.Open(aquifer)
sql = 'select * from judeteEPSG3844 where DENJUD = CLUJ'
aqLayer = aqDS.ExecuteSQL(sql)
feat = aqLayer.GetFeature(0)
print aqLayer.GetFeatureCount()
print feat.GetField('DENJUD')
aqDS.ReleaseResultSet(aqLayer)
The result:
marius@mobi:~/temp$ run_sql.py
1
TELEORMAN
marius@mobi:~/temp$
This explains why none of my centroids were remotely close to the
polygon.
Is there something wrong with my query or should I file a bug?
-marius
On Mon, 2011-07-11 at 20:46 -0500, Marius Jigmond wrote:
I suppose a piece of code speaks louder :):
xsect = False
for i in range(gridLayer.GetFeatureCount()):
feat = gridLayer.GetFeature(i)
geom = feat.GetGeometryRef()
point = geom.Centroid()
for j in range(aqLayer.GetFeatureCount()):
aqfeat = aqLayer.GetFeature(j)
aqgeom = aqfeat.GetGeometryRef()
if point.Intersect(aqgeom):
xsect = True
print 'ibound = 1'
break
if xsect:
feat.SetField('IBOUND', 1)
gridLayer.SetFeature(feat)
xsect = False
-marius
On Mon, 2011-07-11 at 20:34 -0500, Marius Jigmond wrote:
Hi everyone,
I am trying to test whether centroids of polygons lie/intersect
within another polygon. I have tried Intersect, Within, and Contains
but they always return false. Should these methods work for my
intended purpose or do I need to implement a point in polygon
function? Thanks.
-marius
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Unwanted partial transparency when clipping
1.8 compiled from source -marius On Fri, 2011-07-08 at 18:02 -0700, Michael Corey wrote: That certainly could be. I'm running GDAL 1.8 from the kyngchaos site. Which version are you running? Thanks, Michael Corey On 7/8/11 5:56 PM, Marius Jigmond wrote: On Fri, 2011-07-08 at 08:14 -0700, Michael Corey wrote: OK, I did a little more work on this, and I've narrowed down what's going on, but I could still use some help in figuring out how to solve it. Here's my original image: http://mikejcorey.com/spatial/diablo-orig-5pct.tif The original appears to have an RGBA setup (RGB channels and an alpha channel). When I run this: gdalwarp -crop_to_cutline -cutline cutout.shp sourceimage.tif diablo-cutline.tif Here's what I get: http://mikejcorey.com/spatial/diablo-cutline-5pct.tif This comes out as an RGB image with no alpha channel, with each channel being semi-transparent. Not in my case. Could there be something wrong with your GDAL setup? I drew a shapefile mask and ran the above command but my result is RGBA, the alpha band is not lost. See mask and result here: http://ubuntuone.com/p/13WR/ the output of gdalinfo: marius@mobi:~/Downloads$ gdalinfo diablo_cutl.tif Driver: GTiff/GeoTIFF Files: diablo_cutl.tif Size is 266, 224 Coordinate System is: PROJCS[NAD83 / UTM zone 10N, GEOGCS[NAD83, DATUM[North_American_Datum_1983, SPHEROID[GRS 1980,6378137,298.2572221010002, AUTHORITY[EPSG,7019]], AUTHORITY[EPSG,6269]], PRIMEM[Greenwich,0], UNIT[degree,0.0174532925199433], AUTHORITY[EPSG,4269]], PROJECTION[Transverse_Mercator], PARAMETER[latitude_of_origin,0], PARAMETER[central_meridian,-123], PARAMETER[scale_factor,0.9996], PARAMETER[false_easting,50], PARAMETER[false_northing,0], UNIT[metre,1, AUTHORITY[EPSG,9001]], AUTHORITY[EPSG,26910]] Origin = (693798.721929854014888,3902305.751849883235991) Pixel Size = (20.022668877344305,-20.040546259961307) Metadata: AREA_OR_POINT=Area Image Structure Metadata: INTERLEAVE=PIXEL Corner Coordinates: Upper Left ( 693798.722, 3902305.752) (120d52'12.04W, 35d14'42.43N) Lower Left ( 693798.722, 3897816.669) (120d52'15.84W, 35d12'16.81N) Upper Right ( 699124.752, 3902305.752) (120d48'41.44W, 35d14'38.67N) Lower Right ( 699124.752, 3897816.669) (120d48'45.35W, 35d12'13.05N) Center ( 696461.737, 3900061.211) (120d50'28.67W, 35d13'27.75N) Band 1 Block=266x7 Type=Byte, ColorInterp=Red Mask Flags: PER_DATASET ALPHA Band 2 Block=266x7 Type=Byte, ColorInterp=Green Mask Flags: PER_DATASET ALPHA Band 3 Block=266x7 Type=Byte, ColorInterp=Blue Mask Flags: PER_DATASET ALPHA Band 4 Block=266x7 Type=Byte, ColorInterp=Alpha -marius However, if I do this: gdalwarp -dstalpha -crop_to_cutline -cutline cutout.shp sourceimage.tif diablo-dstalpha-cutline.tif I get this: http://mikejcorey.com/spatial/diablo-dstalpha-cutline-5pct.tif This one is strange, because it appears to be a grayscale image. But when I open it in Photoshop, I see that it actually has 4 alpha channels. I suspect that those are just getting set incorrectly as alpha and are in fact the RGB channels, but can someone explain that behavior or how to fix it? What I want to end up with is a clipped RGB image (or RGBA) image where nodata is transparent and the RGB isn't translucent. Thanks again! Michael Corey On 7/7/11 7:13 AM, Eli Adam wrote: Michael, On 7/6/2011 at 5:35 PM, in message4e14ff42.50...@cironline.org, Michael Coreymco...@cironline.org wrote: Sure, I've uploaded samples here. http://www.mikejcorey.com/spatial/diablo-box-sample.tif http://www.mikejcorey.com/spatial/diablo-cutout-sample.tif I don't notice the semi-transparency in these scaled down images. Perhaps it is the way your viewer reads the mask? These are the same as the images created by the process I described (but scaled down). To your point about specifying size in the first step -- will that make the process run faster, or does it do the scaling down after it builds the full-resolution image? Also, I notice that my filesize always gets significantly bigger when I do the cutout step, which seems counter-intuitive to me since in theory shouldn't there be less information present once the cutout is done? -cutline does not 'discard' any data. The extent of the data remains the same unless you reset those extents. You can do that with -crop_to_cutline. Here are some details from the gdalwarp page, http://gdal.org/gdalwarp.html : -crop_to_cutline: (GDAL= 1.8.0) Crop the extent
Re: [gdal-dev] Unwanted partial transparency when clipping
I use QGIS which knows how to interpret the bands. I have no idea how Photoshop handles the alpha band. You can try specifying -dstnodata: gdalwarp -crop_to_cutline -cutline mask.shp -dstnodata 0 0 0 source.tif dest.tif this should output a RGBA tif with NoData=0 (black) and maybe Photoshop can interpret the NoData tag. You might also want to try nearblack: http://www.gdal.org/nearblack.html I suppose it would be important to know what the final purpose of these images is because as far as GIS software goes just about all can read the cropped tif just fine with proper transparency. -marius On Fri, 2011-07-08 at 18:24 -0700, Michael Corey wrote: Actually when I opened up your sample file in Photoshop I got the same results as I had earlier. What viewer are you using to look at your output files? Also, I did find a workaround that at least improves the situation for me. First I translate the source images: gdal_translate source1.tif source1-noalpha.tif -b 1 -b 2 -b 3 -mask 4 -co COMPRESS=JPEG -co PHOTOMETRIC=YCBCR --config GDAL_TIFF_INTERNAL_MASK YES Then use warp to do my merging and clipping: gdalwarp -crop_to_cutline -cutline ~/Documents/GIS/usa/California/doq/diablo-fullzoom-cutout.shp source1-noalpha.tif source2-noalpha.tif fullzoom-clipped.tif I still have to use Photoshop to get rid of the nodata section, but at least my transparency isn't affected. Any ideas how to get rid of the black and have transparent pixels instead, anyone? Thanks again, Michael Corey On 7/8/11 6:12 PM, Marius Jigmond wrote: 1.8 compiled from source -marius On Fri, 2011-07-08 at 18:02 -0700, Michael Corey wrote: That certainly could be. I'm running GDAL 1.8 from the kyngchaos site. Which version are you running? Thanks, Michael Corey On 7/8/11 5:56 PM, Marius Jigmond wrote: On Fri, 2011-07-08 at 08:14 -0700, Michael Corey wrote: OK, I did a little more work on this, and I've narrowed down what's going on, but I could still use some help in figuring out how to solve it. Here's my original image: http://mikejcorey.com/spatial/diablo-orig-5pct.tif The original appears to have an RGBA setup (RGB channels and an alpha channel). When I run this: gdalwarp -crop_to_cutline -cutline cutout.shp sourceimage.tif diablo-cutline.tif Here's what I get: http://mikejcorey.com/spatial/diablo-cutline-5pct.tif This comes out as an RGB image with no alpha channel, with each channel being semi-transparent. Not in my case. Could there be something wrong with your GDAL setup? I drew a shapefile mask and ran the above command but my result is RGBA, the alpha band is not lost. See mask and result here: http://ubuntuone.com/p/13WR/ the output of gdalinfo: marius@mobi:~/Downloads$ gdalinfo diablo_cutl.tif Driver: GTiff/GeoTIFF Files: diablo_cutl.tif Size is 266, 224 Coordinate System is: PROJCS[NAD83 / UTM zone 10N, GEOGCS[NAD83, DATUM[North_American_Datum_1983, SPHEROID[GRS 1980,6378137,298.2572221010002, AUTHORITY[EPSG,7019]], AUTHORITY[EPSG,6269]], PRIMEM[Greenwich,0], UNIT[degree,0.0174532925199433], AUTHORITY[EPSG,4269]], PROJECTION[Transverse_Mercator], PARAMETER[latitude_of_origin,0], PARAMETER[central_meridian,-123], PARAMETER[scale_factor,0.9996], PARAMETER[false_easting,50], PARAMETER[false_northing,0], UNIT[metre,1, AUTHORITY[EPSG,9001]], AUTHORITY[EPSG,26910]] Origin = (693798.721929854014888,3902305.751849883235991) Pixel Size = (20.022668877344305,-20.040546259961307) Metadata: AREA_OR_POINT=Area Image Structure Metadata: INTERLEAVE=PIXEL Corner Coordinates: Upper Left ( 693798.722, 3902305.752) (120d52'12.04W, 35d14'42.43N) Lower Left ( 693798.722, 3897816.669) (120d52'15.84W, 35d12'16.81N) Upper Right ( 699124.752, 3902305.752) (120d48'41.44W, 35d14'38.67N) Lower Right ( 699124.752, 3897816.669) (120d48'45.35W, 35d12'13.05N) Center ( 696461.737, 3900061.211) (120d50'28.67W, 35d13'27.75N) Band 1 Block=266x7 Type=Byte, ColorInterp=Red Mask Flags: PER_DATASET ALPHA Band 2 Block=266x7 Type=Byte, ColorInterp=Green Mask Flags: PER_DATASET ALPHA Band 3 Block=266x7 Type=Byte, ColorInterp=Blue Mask Flags: PER_DATASET ALPHA Band 4 Block=266x7 Type=Byte, ColorInterp=Alpha -marius However, if I do this: gdalwarp -dstalpha -crop_to_cutline -cutline cutout.shp sourceimage.tif diablo-dstalpha-cutline.tif I get this: http://mikejcorey.com/spatial/diablo-dstalpha-cutline-5pct.tif
Re: [gdal-dev] Mask and GEOTiff
In QGIS you can use Raster Calculator to multiply the rasters. -marius On Thu, 2011-04-28 at 00:30 -0700, canduc17 wrote: Hi everyone. Is there a way to crop a GEOTiff image using a Mask of the same dimension? For example, a mask with a filled circle: every pixel inside the circle is 1 and the others are 0. I would like to obtain a new GEOTiff with only the pixel inside the circle visible and the other set to 0. Is that possible with tools like gdal_translate or other ones? Thanks in advance! -- View this message in context: http://osgeo-org.1803224.n2.nabble.com/Mask-and-GEOTiff-tp6312284p6312284.html Sent from the GDAL - Dev mailing list archive at Nabble.com. ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] FGDB Opening Sample File
They are links between tables. Usually, between the attribute table of a feature class and a simple (non-spatial, as Mike said) table. But you can relate any kind of attribute tables. -marius On Thu, 2011-04-28 at 14:30 -0700, Paul Ramsey wrote: Are they tables or links between tables? I can see this is going to be a hard project to do without ArcGIS handy :) P. On Thu, Apr 28, 2011 at 10:26 AM, Smith, Michael ERDC-CRREL-NH michael.sm...@usace.army.mil wrote: Relations are joins to non-spatial tables containing domain values etc. They are relational tables. Mike -- Michael Smith Remote Sensing/GIS Center US Army Corps of Engineers On 4/28/11 1:19 PM, Even Rouault even.roua...@mines-paris.org wrote: Le jeudi 28 avril 2011 02:11:30, Paul Ramsey a écrit : I assume the failures are relationships we don't handle yet (or ever and will need to be silenced). Not sure what those relations things are. There is a test file in the FGDB API itself, and it fails to open entirely, which seems odd, [pramsey@localhost data]$ ogrinfo ./TestData.gdb ERROR 1: GDB Error: Failed to open Geodatabase long:-2147467259 ERROR 1: GDB Error: Failed to open Geodatabase long:-2147467259 FAILURE: Unable to open datasource `./TestData.gdb' with the following drivers. Do other folks have trouble opening the ESRI sample file? The error code doesn't make much sense, it is supposed to mean If the path is seriously in error, say pointing to the wrong drive, a -2147467259 (E_FAIL) error is returned and relates to the CreateGeodatabase method, not the OpenGeodatabase method. Yes I did notice that issue too. I assume either the example is corrupted, either it is a bug in the FileGDB API, either... Perhaps reporting to ESRI might be usefull. Or perhaps Ragi has some clues. Even ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Extract Extent or Bounding Box for a subset of shapes for a shapefile
Luca,
You can also loop through the features with relatively simple Python
code to get the extents. See sample below (adapt at will):
def findPoints(geometry, results):
for i in range(geometry.GetPointCount()):
x,y,z = geometry.GetPoint(i)
if results['north'] == None or results['north'][1] y:
results['north'] = (x,y)
if results['south'] == None or results['south'][1] y:
results['south'] = (x,y)
if results['east'] == None or results['east'][0] x:
results['east'] = (x,y)
if results['west'] == None or results['west'][0] x:
results['west'] = (x,y)
for i in range(geometry.GetGeometryCount()):
findPoints(geometry.GetGeometryRef(i), results)
inSHP = 'test.shp'
shapefile = osgeo.ogr.Open(inSHP)
layer = shapefile.GetLayer(0)
for feat in range(layer.GetFeatureCount()):
feature = layer.GetFeature(feat)
geometry = feature.GetGeometryRef()
results = {'north' : None, 'south' : None, 'east' : None, 'west' :
None}
findPoints(geometry, results)
-marius
On Fri, 2011-04-22 at 17:39 +0200, Luca Sigfrido Percich wrote:
Thank you Chaitanya.
I think that this makes great sense when dealing with a summary of the
whole layer, especially because a lot of formats store the extent in the
layer metadata/header.
But if I need to dump all the features, or to get a dump or summary of a
filtered set of features (in which case I presume ogr will scan all the
features to apply the filter), checking for boundary limits shouldn't
add much cost to the already expensive operation.
Probably this goes beyond the intended scope of ogrinfo. And Matthew's
workaround works fine.
Thanks again
Sig
Il giorno ven, 22/04/2011 alle 20.20 +0530, Chaitanya kumar CH ha
scritto:
Sig,
ogrinfo doesn't force the recomputation of extent. It is considered as
a costly operation. ogrinfo just reports the stored extent.
Whereas in ogr2ogr, creating a new file automatically computes the
extents.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Calculate footprints of shapefiles
Armin, Not sure if this would be straightforward in GDAL but you might want to consider a combination of GDAL + Shapely (http://gispython.org/shapely/docs/1.2/manual.html#cascading-unions). What you're trying to do is spatial analysis not directly implemented in GDAL. Shapely requires that you're working with Cartesian coordinates. The upside is that it has WKT/WKB support for direct loading into PostGIS. -marius On Wed, 2011-03-16 at 01:08 +0100, Armin Burger wrote: Hi all I would like to catalogue shapefiles scattered over lots of directories of the file system and store retrievable information of the shapefiles in a PostGIS layer. Extracting parameters like extent, projection, fields, etc works very fine with GDAL's Python bindings. But I would also like to store a sort of footprint of the whole shapefile as a geometry object since the extent is a bit coarse geographic representation of the shapefile. So far I have no better idea than eg. for polygon shapefiles looping through all features, applying a Union function on them. And at the end trying to use the Simplify method on the resulting polygon that will be used as the footprint. This is for sure not very efficient for larger shapefiles with lots of records. And for line and point shapefiles I still don't have a clue how their records could be represented by an enclosing polygon (maybe the Boundary functions does something like this...). Any ideas how this footprint generation could be achieved in a feasible way using GDAL/OGR Python? Cheers, Armin ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] crop a non-square subwindow?
On Wed, 2011-03-09 at 12:35 -0500, Frank Warmerdam wrote: On 11-03-09 11:51 AM, Wendell Turner wrote: Do any of the gdal utilities crop an image on non-square boundaries? I use gdal_translate with -projwin ulx uly lrx lry to crop an image that is square (left and right edges are vertical, and the top and bottom edges are horizontal). Is there a command that will crop other shapes, for instance a diamond shape, where the edges are not horizontal/vertical? If it must be done by programming (pref. Python), is there a sample that I could use to start with? Wendell, I think gdalwarp with -crop_to_cutline and a cutline specified with the shape will accomplish what you want. http://www.gdal.org/gdalwarp.html The cutline comes from an OGR supported file but this can be as simple as a .csv file with a WKT geometry in it. The -crop_to_cutline may be new to 1.8. Best regards, Wendell, Please be aware that this method may introduce raster shifting and possibly resampling (depending on the original cell size). It might or might not matter for your purposes. A better approach is to use gdal_translate with -projwin to extract to shape envelope and then gdal_rasterize with -i to burn nodata where shape is missing. -marius ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Mosaicking Orthos
On Mon, 2011-02-21 at 13:31 +0530, Chaitanya kumar CH wrote: Zoltan, The creation of 0394w.tif failed. Does the process have write permission on directory? I note that you did not specify the -of option for nearblack. In it's absence, nearblack creates the output file in Erdas Imagine format. True :), regardless of the extension you specify. On Sun, Feb 20, 2011 at 10:50 PM, Zoltan Szecsei zolt...@geograph.co.za wrote: Hi Chaitanya ( all), I tried your nearblack comment, and got: nearblack -white -o 0394w.tif 0394.tif Warning 1: TIFFFetchNormalTag:ASCII value for tag ImageDescription does not end in null byte Warning 1: TIFFFetchNormalTag:ASCII value for tag ImageDescription does not end in null byte Warning 1: TIFFFetchNormalTag:ASCII value for tag Software contains null byte in value; value incorrectly truncated during reading due to implementation limitati ERROR 4: Creation of file 0394w.tif failed. I think I am having trouble understanding the setting/using of nodata values. I have now also downloaded and compiled gdal 1.8.0 (but nearblack remained on v1.7.3) on ubuntu 10.04 LTS, and recreated my vrt, to no improvement. To recap: I have 179 3*8bit RGB band aerial images that have been orthorectified. The resulting skewness in the images has created white slivers in the corners, and I am presuming that these white slivers are recognised as nodata. (My first mistake?) I use the following to create a vrt: gdalbuildvrt -addalpha -hidenodata -srcnodata 255 255 255 -vrtnodata 0 0 0 -input_file_list files 177_untiled_orthos_nodata_180.vrt Thinking that the above will tell gdalbuildvrt that, for the input orthos, white should be recognised as 'nodata', and in the VRT, those white pixels should be converted to black, as nodata. Well, as you've guessed, it didn't happen this way, and I still have the white bits after mosaicking my VRT thus: gdal_translate -projwin -16000 -3520650 -13000 -3525650 177_untiled_orthos_nodata.vrt try1.tif I tried re-running gdalbuildvrt specifying -vrtnodata 255 255 255 (to match the input ortho files) and the outside border of my vrt turned white, but the slivers between the mosaicked orthos were still present. Please can someone show me the correct way to address this problem. I have included gdalinfo results for both an individual ortho, and for the vrt itself. Thanks regards, Zoltan ** GDALINFO on indicidual ORTHO geograph@gs0:~/vrt_177$ gdalinfo --version GDAL 1.8.0, released 2011/01/12 geograph@gs0:~/vrt_177$ geograph@gs0:~/vrt_177$ geograph@gs0:~/vrt_177$ geograph@gs0:~/vrt_177$ gdalinfo /mnt/geo_lvm0/ngi_jobs/orthos_8bit/NGI_177_untiled_orthos/0003.tif Warning 1: TIFFFetchNormalTag:ASCII value for tag ImageDescription does not end in null byte Warning 1: TIFFFetchNormalTag:ASCII value for tag Software contains null byte in value; value incorrectly truncated during reading due to implementation limitations Mask Flags: PER_DATASET ALPHA Band 4 Block=128x128 Type=Byte, ColorInterp=Alpha geograph@gs0:~/vrt_177$ *** On 2011-02-16 19:39, Chaitanya kumar CH wrote: Zoltan, Have you tried the nearblack GDAL utility on the original 176 images? http://www.gdal.org/nearblack.html On Wed, Feb 16, 2011 at 5:25 PM, Zoltan Szecsei zolt...@geograph.co.za wrote: On 2011-02-16 12:20, Zoltan Szecsei wrote: Hi All, I've built a vrt with 176 images. I now need to mosaic this whole area into roughly 5 x 5 ortho images. How can I force the white edges into background of the overlapping orthos? I'm a little cautious of simply (finding a way to) making white transparent as it is only the white on the edges that need to be
Re: [gdal-dev] Coulour balancing and feathering, in a VRT
With gdal_translate you can do per band re-scaling. I've used it on Landsat, it's not magic but it works pretty well. You can group your images, find the right formula for each group and you can quickly verify the results with a VRT. -marius On Mon, 2011-02-21 at 10:39 +0200, Zoltan Szecsei wrote: Hi, I've been googling around a bit, but no luck yet. Now that I have my fancy vrt with a 170 orthos in it, is there a (opensource) product I can colour balance and feather with, before I cut out the required orthophoto regions? Thanks in advance, Zoltan ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Mosaicking Orthos
Hi Zoltan, Are you sure that the white you see is a pure 255/255/255? Albeit, the nearblack with its near default of 15 would have likely taken care of that. If you load two othophotos in QGIS and specify the transparent pixels, does the collar between the orthophotos disappear? Regarding the warnings see: http://www.mail-archive.com/gdal-dev@lists.osgeo.org/msg06960.html -marius ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] gdal_rasterize 1.8.0 options
Hi Even, I'm not trying to re-ignite the stereo debate but I think there's more to the issue with the warning. I've realized that the particular EPSG code 3844 (not sure about others) is reported totally different by ogrinfo and gdalinfo. It's not the actual reporting that's wrong but the way -a_srs is applied by ogr2ogr versus gdal_translate. Below is an example: marius@mobi:~/data/romania/test$ gdalinfo RGB321-K-34-005.tif Driver: GTiff/GeoTIFF Files: RGB321-K-34-005.tif Size is 2985, 2775 Coordinate System is: PROJCS[Pulkovo 1942(58) / Stereo70, GEOGCS[Pulkovo 1942(58), DATUM[Pulkovo_1942_58, SPHEROID[Krassowsky 1940,6378245,298.29998, AUTHORITY[EPSG,7024]], AUTHORITY[EPSG,6179]], PRIMEM[Greenwich,0], UNIT[degree,0.0174532925199433], AUTHORITY[EPSG,4179]], PROJECTION[Oblique_Stereographic], PARAMETER[latitude_of_origin,46], PARAMETER[central_meridian,25], PARAMETER[scale_factor,0.99975], PARAMETER[false_easting,50], PARAMETER[false_northing,50], UNIT[metre,1, AUTHORITY[EPSG,9001]], AUTHORITY[EPSG,3844]] Origin = (96676.434833699619048,290247.117718913941644) Pixel Size = (14.256214353408733,-14.256214353408733) ... marius@mobi:~/data/romania/test$ ogrinfo -where 'FID=0' -geom=summary -fields=no ../caroiaje/caroiaje_100k_EPSG3844.shp caroiaje_100k_EPSG3844 INFO: Open of `../caroiaje/caroiaje_100k_EPSG3844.shp' using driver `ESRI Shapefile' successful. Layer name: caroiaje_100k_EPSG3844 Geometry: Polygon Feature Count: 1 Extent: (94505.736838, 203681.624680) - (905494.263162, 771251.195001) Layer SRS WKT: PROJCS[Stereographic, GEOGCS[GCS_Krassovsky, 1942, DATUM[unknown, SPHEROID[krass,6378245,298.3]], PRIMEM[Greenwich,0], UNIT[Degree,0.017453292519943295]], PROJECTION[Stereographic], PARAMETER[latitude_of_origin,46], PARAMETER[central_meridian,25], PARAMETER[scale_factor,0.99975], PARAMETER[false_easting,50], PARAMETER[false_northing,50], UNIT[Meter,1]] INDICATIV: String (25.0) OGRFeature(caroiaje_100k_EPSG3844):0 POLYGON : 65 points Is the disconnect simply a reporting issue and reprojections are handled properly behind the scenes when using this EPSG code? It seems that what gdalinfo reports is in accordance to http://spatialreference.org/ref/epsg/3844/html/ . -marius On Sat, 2011-02-05 at 00:06 +0100, Even Rouault wrote: Le vendredi 04 février 2011 20:25:43, Marius Jigmond a écrit : Hi Everyone, Some of the options (I didn't test all) introduced with version 1.8.0 seem to behave strangely. Here's what I'm doing: 1. gdal_rasterize -burn 0 -i -l mask mask.shp sample.tif works fine but I get the following warning: Warning : the output raster dataset and the input vector layer do not have the same SRS. Results might be incorrect (no on-the-fly reprojection of input data). even though both inputs use the same projection, EPSG:3844. Unrelated with GDAL 1.8.0. The issue was that the transformation between OGC WKT and ESRI WKT didn't handle well the Oblique Stereographic projection. Fixed in trunk in http://trac.osgeo.org/gdal/changeset/21627 2. gdal_rasterize -burn 0 -i -a_srs EPSG:3844 -l mask mask.shp sample.tif fails with: '-tr xres yes' or '-ts xsize ysize' is required. 3. gdal_rasterize -burn 0 -i -a_nodata 0 -l mask mask.shp sample.tif fails with: '-tr xres yes' or '-ts xsize ysize' is required. From http://www.gdal.org/gdal_rasterize.html it doesn't seem that -a_srs or -a_nodata would require -tr or -ts. Am I doing something wrong? When those 2 new options are used, they imply creating a new dataset. So they actually need to be used with -tr or -ts option. -marius ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] gdal_rasterize 1.8.0 options
Should it then match http://spatialreference.org/ref/epsg/3844/esriwkt/ because it does not seem to? I apologize if my questions seem silly. -marius On Sun, 2011-02-20 at 00:38 +0100, Even Rouault wrote: Le dimanche 20 février 2011 00:24:16, Marius Jigmond a écrit : Is the disconnect simply a reporting issue and reprojections are handled properly behind the scenes when using this EPSG code? It seems that what gdalinfo reports is in accordance to http://spatialreference.org/ref/epsg/3844/html/ . Marius, The .prj file of shapefiles stores WKT in ESRI WKT, not OGC WKT. http://trac.osgeo.org/gdal/changeset/21627 (which is a fix not found in 1.8.0) translates the ESRI Stereographic into the OGC Oblique_Stereographic, to be consistant with what is done in the reverse transformation -marius On Sat, 2011-02-05 at 00:06 +0100, Even Rouault wrote: Le vendredi 04 février 2011 20:25:43, Marius Jigmond a écrit : Hi Everyone, Some of the options (I didn't test all) introduced with version 1.8.0 seem to behave strangely. Here's what I'm doing: 1. gdal_rasterize -burn 0 -i -l mask mask.shp sample.tif works fine but I get the following warning: Warning : the output raster dataset and the input vector layer do not have the same SRS. Results might be incorrect (no on-the-fly reprojection of input data). even though both inputs use the same projection, EPSG:3844. Unrelated with GDAL 1.8.0. The issue was that the transformation between OGC WKT and ESRI WKT didn't handle well the Oblique Stereographic projection. Fixed in trunk in http://trac.osgeo.org/gdal/changeset/21627 2. gdal_rasterize -burn 0 -i -a_srs EPSG:3844 -l mask mask.shp sample.tif fails with: '-tr xres yes' or '-ts xsize ysize' is required. 3. gdal_rasterize -burn 0 -i -a_nodata 0 -l mask mask.shp sample.tif fails with: '-tr xres yes' or '-ts xsize ysize' is required. From http://www.gdal.org/gdal_rasterize.html it doesn't seem that -a_srs or -a_nodata would require -tr or -ts. Am I doing something wrong? When those 2 new options are used, they imply creating a new dataset. So they actually need to be used with -tr or -ts option. -marius ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] error adding ecw on Ubuntu 10.4 (lucid)
If that doesn't work, going back to what Even said earlier: 1. create a gdal_plugins.conf file in /etc/ld.so.conf.d/ with one line: /usr/lib/gdal17plugins 2. sudo ldconfig -marius On Sun, 2011-02-06 at 19:08 +0100, Even Rouault wrote: Le dimanche 06 février 2011 19:00:03, Matt Wilkie a écrit : Check that the gdal_ECW.so (or whatever is name is) is installed in the GDAL plugin directory and then do a ldd on this .so to see if the ECW SDK .so are reachable. You might need to adjust the /etc/ld.so.conf or LD_LIBRARY_PATH to add them. Is this searchable: ? Looks ok. Perhaps try setting the environment variable GDAL_DRIVER_PATH to /usr/lib/gdal17plugins ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Re: gdal_rasterize 1.8.0 options
Even, Hermann might be right. From reading several sources: http://gge.unb.ca/Pubs/TR46.pdf http://www.manifold.net/doc/double_stereographic.htm it seems that double stereographic is oblique stereographic as long as the origin is non-equatorial or non-polar. The actual denomination Double Stereographic seems to be a ESRI by-product that refers to Oblique Stereographic. However, ESRI's own descriptions of these projections are, at best, vague. -marius On Sat, 2011-02-05 at 13:12 +0100, Even Rouault wrote: Le samedi 05 février 2011 12:49:15, Hermann Peifer a écrit : On 05/02/2011 00:06, Even Rouault wrote: Unrelated with GDAL 1.8.0. The issue was that the transformation between OGC WKT and ESRI WKT didn't handle well the Oblique Stereographic projection. Fixed in trunk in http://trac.osgeo.org/gdal/changeset/21627 Even, Once upon a time I opened http://trac.osgeo.org/gdal/ticket/2869, thinking that the reported issue could be fixed with a simple remapping between Double_Stereographic and Oblique_Stereographic projections. Could you perhaps have a look at this 2-year old ticket? Thanks in advance, Hermann Hermann, I've just had a look, but not being neither Dutch nor a specialist of (stereographic) projections, I don't feel competent enough to do any action on it. Those tickets plus the reading of http://udig.refractions.net/files/docs/api- geotools/org/geotools/referencing/operation/projection/Stereographic.html make me deeply confused and wondering if r21627 was really appropriate. ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] error adding ecw on Ubuntu 10.4 (lucid)
Did you run sudo ldconfig? http://trac.osgeo.org/ubuntugis/wiki/UserTutorials -marius On Sat, 2011-02-05 at 20:24 +0100, Even Rouault wrote: Le samedi 05 février 2011 20:18:16, Matt Wilkie a écrit : So you should try sudo gdal-ecw-build /usr/local because you ./configure the ECW SDK without any installation prefix. Ahh, thanks Even. That let me run gdal-ecw-build without error, though there were lots of deprecation warnings. gdalinfo --formats |grep -i ecw still comes back empty however. http://lists.osgeo.org/pipermail/ubuntu/2009-June/54.html indicates this should be enough, though apparently not for me. Any other suggestions? thanks, Check that the gdal_ECW.so (or whatever is name is) is installed in the GDAL plugin directory and then do a ldd on this .so to see if the ECW SDK .so are reachable. You might need to adjust the /etc/ld.so.conf or LD_LIBRARY_PATH to add them. ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] gdalbuildvrt
Aníbal, Have you tried opening the .vrt file with a text editor? For each source .tif there should be three entries in the .vrt file (one per band). What does gdalinfo source.tif say? -marius On Fri, 2011-02-04 at 11:41 -0200, Aníbal Pacheco wrote: On 02/03/2011 04:59 PM, Frank Warmerdam wrote: On 11-02-03 01:30 PM, Aníbal Pacheco wrote: Hi, I'm using gdalbuildvrt to produce a virtual raster layer which can then used in Mapnik and its ogcserver(WMS), everything works except that the result image is black white and with a very poor resolution. The 2 source tif files used in this test were both colored with a good resolution. Any ideas? Aníbal, I would suggest examining the resulting vrt to see if there is anything odd about it. Start with gdalinfo to see if it's resolution seems reasonable and if in fact it is 3 bands. You might want to run gdalinfo -stats xx.vrt to get band statistics just to verify that the bands are radiometrically distinct. Also, check that the band pixel type is reasonable - the same as the source files. If that all looks ok, then perhaps visually examine the VRT with something like QGIS. After that it comes down to talking to Mapnik people about why there are problems. Best regards, I think there is something wrong with the vrt because QGis does not find the 3rd band, here[0] is a screenshot to help. thanks [0] http://imgur.com/OHnFg ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] gdal_rasterize 1.8.0 options
Hi Everyone, Some of the options (I didn't test all) introduced with version 1.8.0 seem to behave strangely. Here's what I'm doing: 1. gdal_rasterize -burn 0 -i -l mask mask.shp sample.tif works fine but I get the following warning: Warning : the output raster dataset and the input vector layer do not have the same SRS. Results might be incorrect (no on-the-fly reprojection of input data). even though both inputs use the same projection, EPSG:3844. 2. gdal_rasterize -burn 0 -i -a_srs EPSG:3844 -l mask mask.shp sample.tif fails with: '-tr xres yes' or '-ts xsize ysize' is required. 3. gdal_rasterize -burn 0 -i -a_nodata 0 -l mask mask.shp sample.tif fails with: '-tr xres yes' or '-ts xsize ysize' is required. From http://www.gdal.org/gdal_rasterize.html it doesn't seem that -a_srs or -a_nodata would require -tr or -ts. Am I doing something wrong? -marius ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] gdal_rasterize 1.8.0 options
Eli, Both input datasets (mask and raster) exist, which is why 1. works. Option 1 doesn't use any of version 1.8.0 flags, which is why I thought something else was introduced behind the scenes with the new flags. -marius On Fri, 2011-02-04 at 12:08 -0800, Eli Adam wrote: Marius, I don't know about 1. but possibilities for 2. and 3 below. Hi Everyone, Some of the options (I didn't test all) introduced with version 1.8.0 seem to behave strangely. Here's what I'm doing: 1. gdal_rasterize -burn 0 -i -l mask mask.shp sample.tif works fine but I get the following warning: Warning : the output raster dataset and the input vector layer do not have the same SRS. Results might be incorrect (no on-the-fly reprojection of input data). even though both inputs use the same projection, EPSG:3844. 2. gdal_rasterize -burn 0 -i -a_srs EPSG:3844 -l mask mask.shp sample.tif fails with: '-tr xres yes' or '-ts xsize ysize' is required. 3. gdal_rasterize -burn 0 -i -a_nodata 0 -l mask mask.shp sample.tif fails with: '-tr xres yes' or '-ts xsize ysize' is required. From http://www.gdal.org/gdal_rasterize.html it doesn't seem that -a_srs or -a_nodata would require -tr or -ts. Am I doing something wrong? Does sample.tif exist? If not, the docs say you need to use -tr or -ts, Since GDAL 1.8.0, the target GDAL file can be created by gdal_rasterize. One of -tr or -ts option must be used in that case. I don't know why that error would not also come from 1. Bests, Eli -marius ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] gdal_rasterize 1.8.0 options
On Sat, 2011-02-05 at 00:06 +0100, Even Rouault wrote: Le vendredi 04 février 2011 20:25:43, Marius Jigmond a écrit : Hi Everyone, Some of the options (I didn't test all) introduced with version 1.8.0 seem to behave strangely. Here's what I'm doing: 1. gdal_rasterize -burn 0 -i -l mask mask.shp sample.tif works fine but I get the following warning: Warning : the output raster dataset and the input vector layer do not have the same SRS. Results might be incorrect (no on-the-fly reprojection of input data). even though both inputs use the same projection, EPSG:3844. Unrelated with GDAL 1.8.0. The issue was that the transformation between OGC WKT and ESRI WKT didn't handle well the Oblique Stereographic projection. Fixed in trunk in http://trac.osgeo.org/gdal/changeset/21627 Thanks. 2. gdal_rasterize -burn 0 -i -a_srs EPSG:3844 -l mask mask.shp sample.tif fails with: '-tr xres yes' or '-ts xsize ysize' is required. 3. gdal_rasterize -burn 0 -i -a_nodata 0 -l mask mask.shp sample.tif fails with: '-tr xres yes' or '-ts xsize ysize' is required. From http://www.gdal.org/gdal_rasterize.html it doesn't seem that -a_srs or -a_nodata would require -tr or -ts. Am I doing something wrong? When those 2 new options are used, they imply creating a new dataset. So they actually need to be used with -tr or -ts option. Thanks for the clarification as it's not apparent from the web page. Would have loved the -a_nodata to be available for existing datasets. -marius ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
RE: [gdal-dev] Gdalwarp -crop_to_cutline results in shifted raster
Submitted as ticket# 3947 but I had to link the sample data because of size. -marius Date: Tue, 1 Feb 2011 19:00:13 +0530 Subject: Re: [gdal-dev] Gdalwarp -crop_to_cutline results in shifted raster From: chaitanya...@gmail.com To: mariusjigm...@hotmail.com CC: gdal-dev@lists.osgeo.org Marius, Can you provide some sample data to investigate this? You can raise a ticket at http://trac.osgeo.org/gdal/newticket and attach the data. On Tue, Feb 1, 2011 at 8:12 AM, Marius Jigmond mariusjigm...@hotmail.com wrote: Hi Everyone, I am using the -crop_to_cutline option of gdalwarp but the resulting raster is slightly shifted. Both cutline shapefile and raster are using the same projection. On the sample raster dataset (pixel=28.5m) the shift is about 5m to the east and 3.5m to the north. Pixel values remain the same (so resampling doesn't take place, as expected). I should add that not adding -crop_to_cutline results in no shift. Commands I ran: (1) gdalwarp -cutline mask.shp -crop_to_cutline sample.tif cropped.tif result is shifted (2) gdalwarp -cutline mask.shp sample.tif out.tif result is not shifted A workaround would be using gdal_rasterize on after option (2) above. -marius ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev -- Best regards, Chaitanya kumar CH. /tʃaɪθənjə/ /kʊmɑr/ +91-9494447584 17.2416N 80.1426E ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] JPEG compressed GeoTIFF ignores Nodata
Thanks for all the suggestions. The following website: http://sites.google.com/site/bpederse/caliwms gave me an idea of how to treat the JPEG artifacts. Going back to our previous discussion, the OpenJPEG implementation of JPEG2000 creates very clean (at the boundary of data/nodata) compressed images when compared to JPEG. But, because it's a new format, I'll have to wait until it's implemented in the jai-imageio-ext package. -marius Le samedi 29 janvier 2011 22:40:20, Marius Jigmond a écrit : The link should work for another day or so (automatic server purge). I can re-upload it if you don't get the chance to download it. Warning: 1.2GB. https://www.twdb.state.tx.us/filetxfr/emaillinkdownload.aspx?id=FileTransf erID=289562464 ok, I've reproduced the bug and fixed it in http://trac.osgeo.org/gdal/ticket/3938 This was one of the most interesting bug since months I had to track down ! Thanks for reporting. Even ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] JPEG compressed GeoTIFF ignores Nodata
Hi Everyone, I am having some trouble figuring out why adding JPEG compression to a RGB GeoTIFF results in Nodata pixel triplets becoming data pixel triplets. Is this expected behavior due to the noisy nature of JPEG compression or to the fact that these are border Nodata triplets? I know certain algorithms process the boundary between data and Nodata differently. I'm running GDAL 1.8.0 with internal libtiff. The command I issued was: gdal_translate -co COMPRESS=JPEG -co TILED=YES -co JPEG_QUALITY=90 -co PHOTOMETRIC=YCBCR s70rgb321.tif test.tif In QGIS using the identify feature on the same triplets results in: s70rgb321.tif: Band 1: null(nodata) Band 2: null(nodata) Band 3: null(nodata) test.tif: Band 1: 1 Band 2: 1 Band 3: 3 gdalinfo s70rgb321.tif yields: Driver: GTiff/GeoTIFF Files: s70rgb321.tif s70rgb321.tif.ovr Size is 32261, 30632 Coordinate System is: PROJCS[Pulkovo_1942_58_Stereo70, GEOGCS[GCS_Pulkovo 1942(58), DATUM[Pulkovo_1942_58, SPHEROID[Krassowsky_1940,6378245,298.3, AUTHORITY[EPSG,7024]], AUTHORITY[EPSG,6179]], PRIMEM[Greenwich,0], UNIT[degree,0.0174532925199433]], PROJECTION[Oblique_Stereographic], PARAMETER[latitude_of_origin,46], PARAMETER[central_meridian,25], PARAMETER[scale_factor,0.99975], PARAMETER[false_easting,50], PARAMETER[false_northing,50], UNIT[metre,1, AUTHORITY[EPSG,9001]]] Origin = (27432.591760581413837,941642.761553813470528) Pixel Size = (28.512429802339831,-28.512429802339831) Metadata: AREA_OR_POINT=Area Image Structure Metadata: INTERLEAVE=PIXEL Corner Coordinates: Upper Left ( 27432.592, 941642.762) ( 18d26'17.19E, 49d47'26.40N) Lower Left ( 27432.592, 68250.012) ( 19d18' 9.44E, 41d58' 4.57N) Upper Right ( 947272.090, 941642.762) ( 31d12'45.10E, 49d48'33.95N) Lower Right ( 947272.090, 68250.012) ( 30d23'36.76E, 41d58'59.57N) Center ( 487352.341, 504946.387) ( 24d50'11.61E, 46d 2'39.82N) Band 1 Block=32261x1 Type=Byte, ColorInterp=Red NoData Value=0 Overviews: 16131x15316, 8066x7658, 4033x3829, 2017x1915, 1009x958, 505x479, 253x240 Metadata: LAYER_TYPE=athematic Band 2 Block=32261x1 Type=Byte, ColorInterp=Green NoData Value=0 Overviews: 16131x15316, 8066x7658, 4033x3829, 2017x1915, 1009x958, 505x479, 253x240 Metadata: LAYER_TYPE=athematic Band 3 Block=32261x1 Type=Byte, ColorInterp=Blue NoData Value=0 Overviews: 16131x15316, 8066x7658, 4033x3829, 2017x1915, 1009x958, 505x479, 253x240 Metadata: LAYER_TYPE=athematic gdalinfo test.tif yields: Driver: GTiff/GeoTIFF Files: test.tif Size is 32261, 30632 Coordinate System is: PROJCS[Pulkovo_1942_58_Stereo70, GEOGCS[GCS_Pulkovo 1942(58), DATUM[Pulkovo_1942_58, SPHEROID[Krassowsky 1940,6378245,298.29998, AUTHORITY[EPSG,7024]], AUTHORITY[EPSG,6179]], PRIMEM[Greenwich,0], UNIT[degree,0.0174532925199433]], PROJECTION[Oblique_Stereographic], PARAMETER[latitude_of_origin,46], PARAMETER[central_meridian,25], PARAMETER[scale_factor,0.99975], PARAMETER[false_easting,50], PARAMETER[false_northing,50], UNIT[metre,1, AUTHORITY[EPSG,9001]]] Origin = (27432.591760581413837,941642.761553813470528) Pixel Size = (28.512429802339831,-28.512429802339831) Metadata: AREA_OR_POINT=Area Image Structure Metadata: SOURCE_COLOR_SPACE=YCbCr COMPRESSION=YCbCr JPEG INTERLEAVE=PIXEL Corner Coordinates: Upper Left ( 27432.592, 941642.762) ( 18d26'17.19E, 49d47'26.40N) Lower Left ( 27432.592, 68250.012) ( 19d18' 9.44E, 41d58' 4.57N) Upper Right ( 947272.090, 941642.762) ( 31d12'45.10E, 49d48'33.95N) Lower Right ( 947272.090, 68250.012) ( 30d23'36.76E, 41d58'59.57N) Center ( 487352.341, 504946.387) ( 24d50'11.61E, 46d 2'39.82N) Band 1 Block=256x256 Type=Byte, ColorInterp=Red NoData Value=0 Metadata: LAYER_TYPE=athematic Band 2 Block=256x256 Type=Byte, ColorInterp=Green NoData Value=0 Metadata: LAYER_TYPE=athematic Band 3 Block=256x256 Type=Byte, ColorInterp=Blue NoData Value=0 Metadata: LAYER_TYPE=athematic Thanks for any insights you might have. -marius ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] JPEG compressed GeoTIFF ignores Nodata
Just realized that earlier I hit the wrong reply button (apologies). The COPY_SRC_OVERVIEWS yielded a wild result. Here's what I did: 1. gdaladdo -ro s70rgb321.tif 2 4 8 16 32 64 128 (note: input GTiff is 3 band with nodata=0 and no mask) 2. gdal_translate -co COMPRESS=JPEG -co TILED=YES -co JPEG_QUALITY=90 -co PHOTOMETRIC=YCBCR -co COPY_SRC_OVERVIEWS=YES s70rgb321.tif test1.tif The result: http://tinypic.com/r/110gig4/7 The top (~)quarter of the original image is replicated downwards up to the original extent of the image. I would have expected the artifacts but not this :). According to the http://www.gdal.org/frmt_gtiff.html COPY_SRC_OVERVIEWS is applied as long as no general options are passed. -marius -- Sent from Ubuntu On Sat, 2011-01-29 at 21:04 +0100, Even Rouault wrote: Le samedi 29 janvier 2011 20:28:02, Marius Jigmond a écrit : Thanks Even. I ran some more tests but that didn't quite work for me. I still get those noisy pixels after gdal_translate (in QGIS I set value 0 for transparency and those pixels stand out). However, I noticed the following: * nearblack -setmask shrinks the significant pixels area (expected?) * the near black noise pixels seem to follow a stair step pattern similar to the edge of the significant area but at a different scale (larger steps, if what I just said makes sense) and extend past the significant data areas in both masked and unmasked datasets. Ah yes, nearblack has some default parameters that can account for this. See the -near and -nb parameters of http://gdal.org/nearblack.html You can set them to 0 to have strict nearblack behaviour. I was aware of the fact that nodata values are band specific but I had a little brain glitch from working with gdaladdo (which honors nodata only as a triplet nodata). That's how I actually ran into the problem. I am trying to compress a RGB GeoTIFF using JPEG and then generate compressed overviews. We plan to serve some imagery through Geoserver and would like to stick to open-source solutions. The compressed datasets are significantly smaller and with overviews they are extremely responsive, rendering very fast). OpenJPEG v2 looks great but I haven't found support for it in Geoserver, yet. To limit the issues with artifacts, I would compute the overviews on the uncompressed dataset and then translate the whole full resolution + overviews into a JPEG compressed TIFF with the new COPY_SRC_OVERVIEWS option of the GTiff driver. I'm not sure that JPEG2000 would solve the issue. This would need to be checked if alpha channel is strictly preserved. For the record, MapServer is able to use nodata mask band as transparency channel. -marius On Sat, 2011-01-29 at 16:13 +0100, Even Rouault wrote: Le samedi 29 janvier 2011 15:47:43, Marius Jigmond a écrit : Hi Everyone, I am having some trouble figuring out why adding JPEG compression to a RGB GeoTIFF results in Nodata pixel triplets becoming data pixel triplets. Is this expected behavior due to the noisy nature of JPEG compression or to the fact that these are border Nodata triplets? I know certain algorithms process the boundary between data and Nodata differently. Marius, Yes you did find the cause : the pixels maching nodata values will be JPEG compressed, so you have no guarantee that the 0 value is preserved through compression/decompression. I'd also note that the nodata concept in GDAL is a per-band value (so it should not be interpreted as a nodata triplet) A possibility to workaround this would be to create a nodata mask band (see http://trac.osgeo.org/gdal/wiki/rfc15_nodatabitmask), but I'm pretty sure that QGIS ignores those mask bands. Anyway, if you want to try and create one, you can : 1) Use nearblack to create a nodata mask band based on (0,0,0) triplet nearblack -setmask -of GTiff -o out_with_mask.tif in.tif (the -setmask option is new in GDAL 1.8.0) 2) Then compress the dataset with gdal_translate -co COMPRESS=JPEG -co TILED=YES -co JPEG_QUALITY=90 -co PHOTOMETRIC=YCBCR out_with_mask.tif out_compressed.tif Only the RGB channels will be JPEG compressed : the mask band will use lossless compression so you won't have artifacts at the borders between significant areas and nodata areas. (You can add --config GDAL_TIFF_INTERNAL_MASK YES to make the nodata mask internal to the GTiff file, instead of having a .tif.msk external file, which will not make any difference for software using GDAL API but can be more convenient to have a sinle file) FYI, you can later reconvert the dataset into an uncompressed RGBA dataset : gdal_translate out_compressed.tif out_rgba.tif -b 1 -b 2 -b 3 -b mask Best regards, Even ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http
Re: [gdal-dev] JPEG compressed GeoTIFF ignores Nodata
The link should work for another day or so (automatic server purge). I can re-upload it if you don't get the chance to download it. Warning: 1.2GB. https://www.twdb.state.tx.us/filetxfr/emaillinkdownload.aspx?id=FileTransferID=289562464 -marius -- Sent from Ubuntu On Sat, 2011-01-29 at 22:25 +0100, Even Rouault wrote: Le samedi 29 janvier 2011 22:17:45, Marius Jigmond a écrit : Just realized that earlier I hit the wrong reply button (apologies). The COPY_SRC_OVERVIEWS yielded a wild result. Here's what I did: 1. gdaladdo -ro s70rgb321.tif 2 4 8 16 32 64 128 (note: input GTiff is 3 band with nodata=0 and no mask) 2. gdal_translate -co COMPRESS=JPEG -co TILED=YES -co JPEG_QUALITY=90 -co PHOTOMETRIC=YCBCR -co COPY_SRC_OVERVIEWS=YES s70rgb321.tif test1.tif The result: http://tinypic.com/r/110gig4/7 The top (~)quarter of the original image is replicated downwards up to the original extent of the image. I would have expected the artifacts but not this :). According to the http://www.gdal.org/frmt_gtiff.html COPY_SRC_OVERVIEWS is applied as long as no general options are passed. This is weird. Do you have a link so I can download s70rgb321.tif and test ? -marius Le samedi 29 janvier 2011 20:28:02, Marius Jigmond a écrit : Thanks Even. I ran some more tests but that didn't quite work for me. I still get those noisy pixels after gdal_translate (in QGIS I set value 0 for transparency and those pixels stand out). However, I noticed the following: * nearblack -setmask shrinks the significant pixels area (expected?) * the near black noise pixels seem to follow a stair step pattern similar to the edge of the significant area but at a different scale (larger steps, if what I just said makes sense) and extend past the significant data areas in both masked and unmasked datasets. Ah yes, nearblack has some default parameters that can account for this. See the -near and -nb parameters of http://gdal.org/nearblack.html You can set them to 0 to have strict nearblack behaviour. I was aware of the fact that nodata values are band specific but I had a little brain glitch from working with gdaladdo (which honors nodata only as a triplet nodata). That's how I actually ran into the problem. I am trying to compress a RGB GeoTIFF using JPEG and then generate compressed overviews. We plan to serve some imagery through Geoserver and would like to stick to open-source solutions. The compressed datasets are significantly smaller and with overviews they are extremely responsive, rendering very fast). OpenJPEG v2 looks great but I haven't found support for it in Geoserver, yet. To limit the issues with artifacts, I would compute the overviews on the uncompressed dataset and then translate the whole full resolution + overviews into a JPEG compressed TIFF with the new COPY_SRC_OVERVIEWS option of the GTiff driver. I'm not sure that JPEG2000 would solve the issue. This would need to be checked if alpha channel is strictly preserved. For the record, MapServer is able to use nodata mask band as transparency channel. -marius On Sat, 2011-01-29 at 16:13 +0100, Even Rouault wrote: Le samedi 29 janvier 2011 15:47:43, Marius Jigmond a écrit : Hi Everyone, I am having some trouble figuring out why adding JPEG compression to a RGB GeoTIFF results in Nodata pixel triplets becoming data pixel triplets. Is this expected behavior due to the noisy nature of JPEG compression or to the fact that these are border Nodata triplets? I know certain algorithms process the boundary between data and Nodata differently. Marius, Yes you did find the cause : the pixels maching nodata values will be JPEG compressed, so you have no guarantee that the 0 value is preserved through compression/decompression. I'd also note that the nodata concept in GDAL is a per-band value (so it should not be interpreted as a nodata triplet) A possibility to workaround this would be to create a nodata mask band (see http://trac.osgeo.org/gdal/wiki/rfc15_nodatabitmask), but I'm pretty sure that QGIS ignores those mask bands. Anyway, if you want to try and create one, you can : 1) Use nearblack to create a nodata mask band based on (0,0,0) triplet nearblack -setmask -of GTiff -o out_with_mask.tif in.tif (the -setmask option is new in GDAL 1.8.0) 2) Then compress the dataset with gdal_translate -co COMPRESS=JPEG -co TILED=YES -co JPEG_QUALITY=90 -co PHOTOMETRIC=YCBCR out_with_mask.tif out_compressed.tif Only the RGB channels will be JPEG compressed : the mask band will use lossless compression so you won't have artifacts at the borders between significant areas and nodata areas. (You
[gdal-dev] Gdalwarp -crop_to_cutline results in 0 filled raster
Hi Everyone, Has anyone experienced the problem in the subject line? I ran: $gdalwarp -cutline shapefile -crop_to_cutline -s_srs EPSG:31700 -t_srs EPSG:3844 source.tif dest.tif and I ended up with a zero filled raster with the correct extent however. Using 1.8dev. I have a rudimentary script that does the extract by mask similar to -crop_to_cutline feature, and I was looking forward to this much simpler approach. -marius -- Sent from Ubuntu ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
