Re: [gdal-dev] GDAL precision for Geometry object
On 2015-10-16 18:20, Even Rouault wrote: Le vendredi 16 octobre 2015 17:57:10, Hermann Peifer a écrit : On 2015-10-12 11:28, Even Rouault wrote: I've changed to '%.15g' formatting for WKT in https://trac.osgeo.org/gdal/ticket/6145 Hi Even, Why did you choose '%.15g', rather than '%.17g'? The latter format should guarantee a lossless binary->decimal->binary round-trip, whereas the former one might also do -- or not. I don't think that any precision loss between '%.15g', and '%.17g' would result into "a real life problem", I am asking just out of curiosity. Hi Hermann, re'adding the list in CC. Your point is valid. My idea was to avoid making WKT returned by the modified version potentially longer than it is currently I am pretty sure you know: '%.17g' could also result into shorter string values, compared to the earlier '%.15f', as in the below (admittedly unfair) example. Hermann $ printf "%.15f\n" 50.375 50.375 $ $ printf "%.17g\n" 50.375 50.375 ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] GDAL precision for Geometry object
On 2015-10-10 13:45, Dimitrianos Savva wrote: For example the geometry POINT (4799826.09861662145704 2773995.445373429451138) cannot be represented by two double types. But still the ExportToWKT function produces the correct result with full precision for the values. While thinking twice about the issue (I might be wrong, though): Isn't it more of a generic problem that by far most reals (of which there are infinite many) cannot be *exactly* represented by half, single, double, quadruple-precision, or by whatever other floating-point format that could ever be invented in the future? In these cases, the "most accurate representation" must simply be good enough, as one can't do much about it? Your sample value above is 4799826.09861662145704 (decimal): The most accurate representation as double is: 4.79982609861662145704030990601E6 (decimal), 0x41524F54864FBC17 (hex) and : 0101 01010010 0100 01010100 1110 0100 1000 00010111 (binary). That's just how it is, or is there another problem, of which I wouldn't be aware? I took the "most accurate" values from http://www.binaryconvert.com/result_double.html @Even: I also tried to imagine a reasonable real-life coordinate value that would generate scientific notation when printed with "%.15g" or: "%.17g" (as hamish mentioned in this old ticket). Can you point me to one? Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] GDAL precision for Geometry object
On 2015-10-11 19:20, Dimitrianos Savva wrote: I just wondering if GDAL make use of any kind of custom-user-defined type that provides higher precision from a double. My question wasn’t so clear. The problem started comparing GDAL’s and GeoTools’ (Java) WKT representations of the same geometries that has been read from the same Shapefile. I noticed that GDAL provides higher accuracy than GeoTools. So I asked this question here, wondering whether GDAL is using something else than just doubles in memory. But I also concluded that in fact JAVA and C doubles have the same representation in memory by running some minimal tests. However in Java, when converting a double into string, it even with (“.15f”) format, it does not give the same string value as doing so in C. Hmm. I'm not a programmer of any kind and only know doubles and their siblings as defined by IEEE 754. I doubt that any serious programming language could have an interest in ignoring this standard. > GDAL provides higher accuracy.. Hmm. My observation is rather that when exporting into text formats (AAIGrid, WKT, etc.), GDAL seems to print more "numbers" (actually: longer string values) than meaningful in the given context. A double is a double is a double, and can only be as accurate as defined by IEEE 754. I am quoting hamish once again: "I'd point out that for IEEE double precision floating point you don't gain anything by going past %.17g." Just to quote another source (Wikipedia): If an IEEE 754 double precision is converted to a decimal string with at least 17 significant digits and then converted back to double, then the final number must match the original.[1] I assume it is clear now that the conclusion "I see many numbers", so this format "must have higher accuracy" can indeed be wrong. Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] List of reasonable default parameters for each projection
On 2015-10-11 21:09, Ray Gardener wrote: I was wondering if there is a list, for each projection supported by OGR, reasonable default values for each projection parameter. For example, if the user is switching projections in an app, a small table listing the projection parameters appears, but they need to be populated with reasonable default values. e.g. for the New Zealand Mining Grid, one could do spatialRef.SetNZMG(centerLat = 41, centerLong = 173, falseEasting=251, falseNorthing=6023150); [peifer@moby:gdal]> pwd /path/to/local/share/gdal [peifer@moby:gdal]> grep 6023150 *.csv pcs.csv:27200,"NZGD49 / New Zealand Map Grid",9001,4272,19917,9811,1,0,4400,8801,-41,9102,8802,173,9102,8806,251,9001,8807,6023150,9001, projop_wparm.csv:19917,New Zealand Map Grid,9811,8801,-41,9102,8802,173,9102,8806,251,9001,8807,6023150,9001, Seems easy, Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] GDAL precision for Geometry object
On 2015-10-10 17:56, Hermann Peifer wrote:
On 2015-10-10 14:06, Even Rouault wrote:
The implementation of ExportToWKT() in GDAL uses 15 decimal figures
("%.15f"
printf formatting) whatever the source precision was, and as you point,
whatever the intrinsic precision of the double value is.
Does "%.15f" have any added value over "%.15g" (in the usual real-life
GIS situations, I mean). I also vaguely remember some posting/ticket
about "%.15f" versus "%.15g", including a request to partially revert a
changeset. This was years ago, and I can't find at the moment.
The earlier issue was related to some change in the AAIGrid writer,
http://trac.osgeo.org/gdal/ticket/3732. As hamish (from the GRASS
community) wrote at the end of the ticket: "I'd point out that for IEEE
double precision floating point you don't gain anything by going past
%.17g."
Hermann
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] GDAL precision for Geometry object
On 2015-10-10 14:06, Even Rouault wrote:
The implementation of ExportToWKT() in GDAL uses 15 decimal figures ("%.15f"
printf formatting) whatever the source precision was, and as you point,
whatever the intrinsic precision of the double value is.
Does "%.15f" have any added value over "%.15g" (in the usual real-life
GIS situations, I mean). I also vaguely remember some posting/ticket
about "%.15f" versus "%.15g", including a request to partially revert a
changeset. This was years ago, and I can't find at the moment.
Hermann
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] assigning vertical datum to image files
About your initial problem: but gdalbuildvrt complained, informing me that they were not in the same projection. What about using the below gdalbuildvrt option? Hermann -allow_projection_difference: (starting with GDAL 1.7.0) When this option is specified, the utility will accept to make a VRT even if the input datasets have not the same projection. Note: this does not mean that they will be reprojected. Their projection will just be ignored. Source: http://www.gdal.org/gdalbuildvrt.html On 2015-06-08 21:01, Newcomb, Doug wrote: Hi Folks, I have a directory of 783 .img format Lidar - based DEMs in the same projection, North Carolina State Plane Feet (2011) , NAD 83 , NVD88. I was going to use gdalbuildvrt to create a single virtual image for the area, but gdalbuildvrt complained, informing me that they were not in the same projection. A couple of quick bash scripts/commands for x in *.img; do gdalinfo $x img_info.txt;done and grep PROJCS img_info.txt|sort|uniq -c gives me the following: 437 PROJCS[NAD_1983_2011_StatePlane_North_Carolina_FIPS_3200_Ft_US, 346 PROJCS[NAD_1983_StatePlane_North_Carolina_FIPS_3200_Feet_2011, gdalinfo gives the following for each type of file ( origin section clipped out) : Driver: HFA/Erdas Imagine Images (.img) Files: D05_37_20878102_20141231.img D05_37_20878102_20141231.img.aux.xml Size is 1000, 1000 Coordinate System is: PROJCS[NAD_1983_2011_StatePlane_North_Carolina_FIPS_3200_Ft_US, GEOGCS[GCS_NAD_1983_2011, DATUM[NAD_1983_2011, SPHEROID[GRS_1980,6378137.0,298.257222101]], PRIMEM[Greenwich,0.0], UNIT[Degree,0.017453292519943295]], PROJECTION[Lambert_Conformal_Conic_2SP], PARAMETER[False_Easting,200.0], PARAMETER[False_Northing,0.0], PARAMETER[Central_Meridian,-79.0], PARAMETER[Standard_Parallel_1,34.3], PARAMETER[Standard_Parallel_2,36.16667], PARAMETER[Latitude_Of_Origin,33.75], UNIT[Foot_US,0.30480060960121924], VERTCS[NAVD_1988_Foot_US, VDATUM[North_American_Vertical_Datum_1988], PARAMETER[Vertical_Shift,0.0], PARAMETER[Direction,1.0], UNIT[Foot_US,0.3048006096012192]]] Driver: HFA/Erdas Imagine Images (.img) Files: D05_37_20957301_20141231.img D05_37_20957301_20141231.img.aux.xml Size is 1000, 1000 Coordinate System is: PROJCS[NAD_1983_StatePlane_North_Carolina_FIPS_3200_Feet_2011, GEOGCS[GCS_NAD_1983_2011, DATUM[NAD_1983_2011, SPHEROID[GRS_1980,6378137.0,298.257222101]], PRIMEM[Greenwich,0.0], UNIT[Degree,0.0174532925199433]], PROJECTION[Lambert_Conformal_Conic_2SP], PARAMETER[False_Easting,200.00261], PARAMETER[False_Northing,0.0], PARAMETER[Central_Meridian,-79.0], PARAMETER[Standard_Parallel_1,34.34], PARAMETER[Standard_Parallel_2,36.16], PARAMETER[Latitude_Of_Origin,33.75], UNIT[Foot_US,0.30480060960121924], VERTCS[NAVD_1988_Foot_US, VDATUM[North_American_Vertical_Datum_1988], PARAMETER[Vertical_Shift,0.0], PARAMETER[Direction,1.0], UNIT[Foot_US,0.3048006096012192]]] In theory, these should all be EPSG:6543, so just assigning the correct horizontal projection/datum with the epsg code should make things usable. (i.e, gdal_translate -a_srs epsg:6543 --config GDAL_CACHEMAX 512 -of GTiff -co COMPRESS=DEFLATE -co PREDICTOR=3 in.img out.tif ) ( Thank you for the EPSG:6543 projection support in GDAL 2.0!) However, this only assigns the horizontal infomation, how does one assign a vertical datum with a horizontal EPSG code? I see the NVD88 code in vertcs.csv , but how would I implement it in the above command? Using gdal 2.0 Beta2 Doug -- Doug Newcomb USFWS Raleigh, NC 919-856-4520 ext. 14 doug_newc...@fws.gov mailto:doug_newc...@fws.gov - The opinions I express are my own and are not representative of the official policy of the U.S.Fish and Wildlife Service or Dept. of the Interior. Life is too short for undocumented, proprietary data formats. ___ 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] Can anyone else replicate this bug?
On 2014-09-03 1:24, Ravi Desai wrote: Hello everyone, I unziped this file: http://www12.statcan.gc.ca/census-recensement/2011/geo/bound-limit/files-fichiers/gct_000a11a_e.zip , and then fed the directory into the following command: ogr2ogr -t_srs epsg:4326 -f CSV output.csv canada_census_boundaries/ -lco GEOMETRY=AS_WKT Line 4 ends with ...,-76.52966604442 44.24964776838)),5210015.00,0015.00,521,Kingston,B,35521,35,Ontario but line 5 ends with ...-81.06617526813 42.,5550100.01,0100.01,555,London,B,3,35,Ontario As you can see, the geometry is truncated and the ending quotes are also not applied to it, creating invalid CSV. See the post OGR v1.11 CSV GEOM=AS_WKT truncated..., at http://lists.osgeo.org/pipermail/gdal-dev/2014-June/thread.html The issue has been fixed already: r27434 | rouault | 2014-06-04 21:38:05 +0200 (Wed, 04 Jun 2014) CSV: fix to avoid truncation of WKT geometries to 8000 characters I also mentioned a work-around in my reply back then. Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] OGR v1.11 CSV GEOM=AS_WKT truncated to 8000 characters (works normally in v1.10.1)
On 2014-06-03 17:45, lucvanlinden wrote: Hi All With upgrading to ogr 1.11 it seems that OGR2OGR -lco geom=AS_WKT, when used in writing to a CSV output format, the WKT notation is truncated to a length of 8000 characters. This is/was not the case with ogr 1.10.1. Can this be confirmed as being a bug? Same here with GDAL 1.12dev from trunk, on Mac OS This works fine and prints lines with 1.5+M characters, with my sample data: ogr2ogr -f csv out.csv source_data -sql select *,OGR_GEOM_WKT from layer Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] what version of epsg-db in gdal 1.11 ?
On 2014-05-06 15:04, Andrea Peri wrote: Hi, The IGM (Italian Geography Military) has update the definition of the Reference System in the EPSG DB. So now there some new codes usable for datasets in Italy: epsg::6007, 6008, and so on. I just checked at http://epsg-registry.org/ retrieve by code. There don't seem to be any CRS codes in this range. Their online database version is 8.4 Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] what version of epsg-db in gdal 1.11 ?
Aha. These codes are not included in GDAL's most recent gcs.csv/pcs.csv, see http://trac.osgeo.org/gdal/browser/trunk/gdal/data Hermann On 2014-05-06 21:30, Andrea Peri wrote: Oops, you have right I wrong to report the codes. the new codes are from 6704 to 6711 so the right for test was 6707 not 6007 Thx, Andrea. 2014-05-06 21:19 GMT+02:00 Hermann Peifer pei...@gmx.eu mailto:pei...@gmx.eu: On 2014-05-06 15:04, Andrea Peri wrote: Hi, The IGM (Italian Geography Military) has update the definition of the Reference System in the EPSG DB. So now there some new codes usable for datasets in Italy: epsg::6007, 6008, and so on. I just checked at http://epsg-registry.org/ retrieve by code. There don't seem to be any CRS codes in this range. Their online database version is 8.4 Hermann -- - Andrea Peri . . . . . . . . . qwerty àèìòù - ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] what version of epsg-db in gdal 1.11 ?
On 2014-05-06 15:04, Andrea Peri wrote: Hi, The IGM (Italian Geography Military) has update the definition of the Reference System in the EPSG DB. So now there some new codes usable for datasets in Italy: epsg::6007, 6008, and so on. I just checked at http://epsg-registry.org/ retrieve by code. There don't seem to be any CRS codes in this range. Their online database version is 8.4 Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] GML enhancement proposal: Getting all XML attributes as OGR fields
Dear All, I am wondering if someone else would be interested in a config option like: GML_ATTRIBUTES_TO_OGR_FIELDS that would enable auto-detection of all XML attributes and subsequently write them into the .gfs file, as described at http://trac.osgeo.org/gdal/ticket/5418 If there is no broader interest in such a feature, then I should close the ticket as invalid, I would guess. Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Extra Decimal places added to ASCII gdal_translate mosaic
On 2014-01-24 17:22, Norman Vine wrote: hmm this seems to be related to this fairly recent change http://trac.osgeo.org/gdal/ticket/3732 http://osgeo-org.1560.x6.nabble.com/gdal-dev-Change-of-DECIMAL-PRECISION-in-AAIGrid-td5075524.html I am just repeating what I wrote in the above-mentioned mail from 3 Sep 2013: --- snip --- DECIMAL_PRECISION actually became a misnomer now and should perhaps be called SIGNIFICANT_DIGITS -- or the format should be changed back to %precisionf. Maybe some of the GDAL devs has an opinion about the issue. --- snip --- Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] gdalinfo coordinates problem?
As far as I can see: lon_0 is the same in both cases and there is also no actual difference between +datum=NAD83 (gdal) and +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 (qgis). I am not sure where +x_0=60 (qgis) comes from, but gdal's +x_0 value is simply 6458320.41666 * 0.3048006096012192 = 1968500 *meters*, so in fact only +units=us-ft is wrong there. Hermann On 2014-01-18 20:09, Etienne Tourigny wrote: I noticed that qgis-2.0 also has a problem with this file, but places it around 54 degrees west longitude (if using on-the-fly reprojection to WGS84). The CRS is different from that identified by gdal (the lon_0, units and datum are different) qgis: +proj=lcc +lat_1=40.97 +lat_2=39.93 +lat_0=39.34 +lon_0=-77.75 +x_0=60 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs gdal: PROJ.4 : '+proj=lcc +lat_1=40.97 +lat_2=39.93 +lat_0=39.34 +lon_0=-77.75 +x_0=1968500 +y_0=0 +datum=NAD83 +units=us-ft +no_defs ' OGC WKT : PROJCS[NAD_1983_StatePlane_Pennsylvania_South_FIPS_3702_Feet, 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[Lambert_Conformal_Conic_2SP], PARAMETER[standard_parallel_1,40.97], PARAMETER[standard_parallel_2,39.93], PARAMETER[latitude_of_origin,39.34], PARAMETER[central_meridian,-77.75], PARAMETER[false_easting,6458320.41666], PARAMETER[false_northing,0], UNIT[us_survey_feet,0.3048006096012192], AUTHORITY[EPSG,32129]] On Fri, Jan 17, 2014 at 8:37 PM, Even Rouault even.roua...@mines-paris.org mailto:even.roua...@mines-paris.org wrote: Le vendredi 17 janvier 2014 21:28:35, Reynolds, Scott a écrit : Hi, I've downloaded several DEMs (e.g 27002570PAS_dem.zip) from ftp://pamap.pasda.psu.edu/pamap_lidar/cycle1/DEM/South/2008/2000/. Gdalinfo is reporting the corner coordinates to be in the vicinity of 91 degrees 31 minutes west longitude which should be more like 75 degrees 36 minutes west longitude. Can someone explain what is happening here? GeoTIFF support for non-meter linear units has been a major source of headaches in GDAL / libgeotiff. People will appreciate looking at frmts/gtiff/gt_wkt_srs.cpp I'm not sure if that GeoTIFF file is correct or not regarding the values of the false easting/northing with respect to the GeoTIFF spec (many in the wild are broken), but when looking at what happens, I noticed that a transformation from meters to feet happened twice. The attached patch fixes that (and the corner coordinates are around 75°), but I'm not sure if it is completely correct in all cases. You should likely create a GDAL Trac ticket with your report and the patch, but I'm not sure who will dare hurting his head against this wall... Even Thanks, Scott -- Geospatial professional services http://even.rouault.free.fr/services.html ___ gdal-dev mailing list gdal-dev@lists.osgeo.org mailto: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] Understanding gdal_rasterize
On 2013-10-21 14:46, Paul Meems wrote: I'm trying to understand how gdal_rasterize works. As a test I'm trying to convert a shapefile with all counties of the USA to a GeoTiff. I've added a field to the shapefile called myID which is equal to the shape number. The goal is to create a tiff with each county its own value. This is my command: gdal_rasterize -a myID -l counties -of GTiff -a_nodata -999 -ts 800 800 counties.shp GDALRasterizeTest.tif You could drop -of GTiff as GTiff is anyway the default format. -a_nodata -999 makes most sense in combination with -init -999. Anything else looks fine to me. I'm having trouble understanding the -ts or -tr parameters. Using the above command I get an 800*800 tiff which shows the USA, but it looks like all pixels are the same color. -ts obviously sets the output file size, to 800x800 pixels, as you noted yourself. Do all pixels also have the same *value* ? Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] gdal_edit.py -a_nodata none
On 2013-08-23 15:23, Even Rouault wrote: Selon peifer pei...@gmx.eu: Hi, http://www.gdal.org/gdal_edit.html states about -a_nodata Assign a specified nodata value to output bands. Can be set to none to remove a nodata value if one exists for the dataset. I assume that none means the literal string `none', but what happens is given below. Is this just a plain bug or am I doing something terribly wrong? Hi Hermann, This is a documentation bug. This is not supported by the code. And there's no (standard) way in the GDAL API to remove a nodata value once it is set. Could you file a ticket about that ? Done, see http://trac.osgeo.org/gdal/ticket/5213 Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] gdalwarp produces all black output
On 2013-08-08 9:55, ludwig.hilger wrote: Hello everybody, I am also a newbie and seem to have a similar problem, but cannot get it solved. I am trying to reproject the following file: http://www.altmuehlnet.de/~hilger/Test_Smooth_0208_4258.tif I am using the following line: gdalwarp -et 0 -r bilinear -s_srs EPSG:4258 -t_srs EPSG:25832 Test_Smooth_0208_4258.tif Test_Smooth_0208_4258_25832.tif I downloaded your test file and noted that it has 1x1 pixel and a real world extent of around 7x7 metres. So 1 pixel in your GeoTIFF represents less than 1x1 mm. I can only blindly assume that some rounding effects related to the small pixel size lead to the problem you observed. I noted the same problem in my testing (I am using the MacPort of GDAL 1.10, on my MacBook). I re-defined the extent of your file to 1x1 degree and was able to warp the resulting file without any problem: $ gdal_translate Test_Smooth_0208_4258.tif out_1x1.tif -a_ullr 10 47 11 46 $ gdalwarp out_1x1.tif out_25832.tif -t_srs epsg:25832 -et 0 -rb Maybe someone from the GDAL devs has an explanation for this behaviour. Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Bad result when writing into KML without -t_srs
On 2013-06-07 6:33, Even Rouault wrote: Selon Jukka Rahkonen jukka.rahko...@mmmtike.fi: Herman Peifer wrote but forgot to send to the list: Jukka, Could you please try with -a_srs epsg:2393 instead of -s_srs epsg:2393, see http://trac.osgeo.org/gdal/ticket/4496; This is funny, but using -a_srs:2393 in the command really leads to correct result. Following command, even being all wrong with the a_srs, gives good KML and it also overrides the erroneous prj file. ogr2ogr -f KML -a_srs epsg:2393 test6.kml test_region.shp This is also explainable by ogr2ogr and the KML driver works. -a_srs attaches a SRS to the feature that is passed to the output driver (overriding the SRS of the source feature if existing). The KML driver has a very unique behaviour in which it takes into account the SRS of the features that are passed to it to do reprojection to EPSG:4326 if needed. Said differently the KML driver has some sort of -s_srs -t_srs mechanism integrated into it. Said even differently, and from my perspective: the KML driver's behaviour is counter-intuitive. However, this seems to be a feature rather than bug, as far as I can gather from Even's explanations. Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] GoogleEarth wont read KML files created by ogr2ogr: Unknown element Schema
FYI, some 4 years ago, I reported an XML schema validation error related to the Schema element: http://trac.osgeo.org/gdal/ticket/2858. IIRC, the explanation at the time was that there are some limitations with the kml driver and the real solution is to switch to libkml... Hermann On 2013-05-01 19:40, Frank Warmerdam wrote: Nikos, Could you file a ticket, and attach a working and nonworking kml example? Any idea if this behavior varies across versions of google earth? I wonder if libkml can read the file(s) ok. Best regards, Frank On Wed, May 1, 2013 at 3:30 PM, Nikos Alexandris n...@nikosalexandris.net mailto:n...@nikosalexandris.net wrote: Hi! I work with GDAL 1.9.2, compiled from source with support for expat and KML under Kubuntu 12.10. Any attempt to convert a(ny) shapefile (be it of type point or polygon) using ogr2ogr (or exporting from GRASS GIS or/and QGIS -- all based on the same GDAL librairies) results in a .kml file that Google- Earth (version 7. ) wont read with the excuse: Unknown element Schema. Even the example KML file over at https://developers.google.com/kml/articles/vector is rejected by GE with the same error message. Removing the Schema entries altogether, makes the file readable by GE. Anyhow, GE reads without problems a KML example (that includes Schema elements) to be found over at https://developers.google.com/kml/documentation/extendeddata (look for the Example under Adding Typed Data to a Feature). What are the odds? Thank you, Nikos ___ gdal-dev mailing list gdal-dev@lists.osgeo.org mailto:gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev -- ---+-- I set the clouds in motion - turn up | Frank Warmerdam, warmer...@pobox.com mailto:warmer...@pobox.com light and sound - activate the windows | http://pobox.com/~warmerdam and watch the world go round - Rush| Geospatial Software Developer ___ 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] creating bounding box map of a contour map
On 2013-01-24 11:44, Ahmet Temiz wrote: hello I have a contour elevation (line) map. I create bounding box map ( as a rectangular polygon ) of this contour map using extend. how can I do that ? If I understood correctly, all you need to do is: ogrtindex output_dataset src_dataset http://www.gdal.org/ogrtindex.html Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] OSGB36 to WGS84 transformation
The lower left corner coordinates of the SD square are: 300 000 metres East 400 000 metres North So your test coordinate is at (Easting Northing): 380150 402109 echo 380150 402109 | gdaltransform -s_srs EPSG:27700 -t_srs EPSG:4326 -2.30083027968938 53.5152719980087 48.8662291513756 Please note that behind the scenes, these datum shift parameters are applied: +towgs84=375,-111,431,0,0,0,0 (which might be what you want .. or not). Hermann On 2013-01-08 21:17, Johnny wrote: Hi, I am trying to convert OSGB grid references to WGS84, but do not succeed and don't understand why? I use : echo SD8015002109 |gdaltransform -s_srs EPSG:27700 -t_srs EPSG:4326 which gives : -7.55729138544866 49.7667307394399 50.480780794 If I use [1], the result is : 53.515279 -2.300819 which can be verified as correct on maps. What error am I making in assuming this should work? I use GDAL 1.9.1, released 2012/05/15 Fedora 17 Thanks! Footnotes: [1] http://www.nearby.org.uk/coord.cgi?p=SD8015002109f=conv ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] OSGB36 to WGS84 transformation
On 2013-01-09 16:13, Johnny wrote: Are you aware of any way to batch process this in Linux to enable swift conversion of multiple references? You'd need a script that maps the 2-letter grid cell identifiers to the respective offsets, then calculates the Easting and Northing values. I can send you an example off-list. Please note that behind the scenes, these datum shift parameters are applied: +towgs84=375,-111,431,0,0,0,0 (which might be what you want .. or not). Sorry, I don't know much about cartography and datums, but guess these are the defaults by gdaltransform? Anyhow, it achieves what I want! :) Indeed, these are the default values for a 3-parameter transformation (I'm using GDAL 1.9.2). On page 33 of [1], OSGB gives the values for a 7-parameter transformation which should be more accurate in most cases. You have to revert the signs of the given values and use them like this: +towgs84=446.448,-125.157,542.06,0.1502,0.247,0.8421,-20.4894 And there is OSGB's grid shift file at [2] which leads to even more accurate results. Hermann [1] http://www.ordnancesurvey.co.uk/oswebsite/gps/docs/A_Guide_to_Coordinate_Systems_in_Great_Britain.pdf [2] http://www.ordnancesurvey.co.uk/oswebsite/gps/osnetfreeservices/furtherinfo/ostn02_ntv2.html Regards, Johnny Footnotes: [1] http://digimap.edina.ac.uk/webhelp/os/gazetteer_plus/grid_ref_conversion.htm ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Updating shapefile fields name with ogr2ogr
On 2012-06-20 19:50, Even Rouault wrote: I know this is going to sound a bit counter-intutive at first... ogrinfo opens datasources in read/write mode by default... Just out of curiosity: why does an ...info tool use read/write mode by default? One would expect that it only needs to *read* some info in the data source? Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: configure error
On 19/03/2012 09:55, Jean-Claude Repetto wrote: Hello , Since this morning, ./configure of the trunk version of GDAL returns an error : configure: error: FileGDBAPI not found. Tested on Linux 32 bits and Linux 64 bits. http://osgeo-org.1560.n6.nabble.com/gdal-dev-compiling-with-FileGeodb-API-1-2-td4623180.html And see the re-opened ticket: http://trac.osgeo.org/gdal/ticket/4570 Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: How to check version of ogr?
On 16/03/2012 10:20, Jean-Claude Repetto wrote: Le 16/03/2012 09:55, Chrishelring a écrit : Hi all, very simple question (atleast I hope so ;). I need to check the version of the ogr library included in my mapserver configuration, but cant seem to find a way to do it. If I just enter ogr2ogr ind the shell no versionnumber is presented. ogr2ogr --version This is not always reliable information. In my case, it says: GDAL 1.9.0, released 2011/12/29 But what I am actually using is GDAL/OGR from trunk, as of today, revision 24121. Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: ogr2ogr GML creation encoding issue
On 24/02/2012 20:05, Even Rouault wrote: CC'ing the list in case someone has a better idea. I don't. Should work... Selon Travis Kirstinetraviskirst...@gmail.com: My Command ogr2ogr -f GML -t_srs EPSG:4326 -s_srs EPSG:900913 -nln administrative planet_osm_polygon_boundary_administrative.gml PG:'dbname=osmdb user=osm' -sql SELECT admin_level, name, way FROM planet_osm_polygon WHERE boundary = 'administrative') Are you using ogr2ogr frim trunk? Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: org2ogr, KML, and srs definition
On 18/02/2012 20:11, Brent Fraser wrote: I have some shapefiles I thought I would look at in Google Earth, so I used ogr2ogr (v1.8.1) to convert them to KML: ogr2ogr -f KML Lines_OGR.kml Lines.shp but they appeared about 100 meters Northeast of my expectation, leading me to believe there was a Datum problem. OGR interprets your .prj file like this: +proj=utm +zone=35 +ellps=intl +units=m +no_defs ...and the EPSG code like this: +proj=utm +zone=35 +ellps=intl +towgs84=-87,-98,-121,0,0,0,0 +units=m +no_defs As there is no +towgs84 parameter in the 1st case, you end up with the 100m shift. In trunk, there is some code that should enable some sort of EPSG code auto-detection, but this is still experimental... Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: Getting SRID from ESRI shapefile
On 06/02/2012 14:52, Etienne Tourigny wrote: I have been working on experimental EPSG lookup, see ticket #4345 http://trac.osgeo.org/gdal/ticket/4345 It requires gdalsrsinfo from gdal-1.9 and some extra data files installed in the gdal data directory (see the ticket) It could perhaps be an idea to make it more clear in the ticket which extra data files are needed (as there is quite a list of attachments). I also tried to run gen_epsg_wkt.sh and ended up with a lot of: ERROR 1: ERROR - wkt_old output not supported I guess this should read: wkt_simple Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: Problems with transform raster image - try again !
Another option would be to gdal_translate to a .vrt by using -a_srs srs_def EPSG:31466 -a_nodata value # in case you happen to know -a_ullr ulx uly lrx lry # coordinates at pixel edges, not their centre points Think about the .vrt as being the big brother of a world file. You might want to spend some time reading the documentation at [1], [2] and [3]. [1] http://www.gdal.org/gdal_translate.html [2] http://www.gdal.org/gdal_vrttut.html [3] http://www.gdal.org/frmt_various.html#WLD Hermann On 25/11/2011 11:49, Jan Tappenbeck wrote: HI ! i try my translation in use of an other tif-source files! here even in the tif-header are informations about the reference - worst coordintes: C:\Program Files (x86)\FWTools2.4.7gdalinfo D:\City\raster\gas\500\g8746-1_ cit.tif Driver: GTiff/GeoTIFF Files: D:\City\raster\gas\500\g8746-1_cit.tif Size is 16215, 8904 Coordinate System is `' Origin = (-948545665.88952053000,-1021132473.1104795) Pixel Size = (3209.779041131308100,-3209.779041131308100) Image Structure Metadata: COMPRESSION=LZW INTERLEAVE=BAND Corner Coordinates: Upper Left (-948545665.890,-1021132473.110) Lower Left (-948545665.890,-1049712345.693) Upper Right (-896499098.738,-1021132473.110) Lower Right (-896499098.738,-1049712345.693) Center (-922522382.314,-1035422409.402) Band 1 Block=16215x1 Type=Byte, ColorInterp=Palette Color Table (RGB with 256 entries) 0: 255,255,255,255 i create following twf-file manuell: 0.032098 0.00 0.00 -0.032098 2586989.395890 5746263.495720 but when i now try to translate by C:\Program Files (x86)\FWTools2.4.7gdalwarp -s_srs EPSG:31466 -t_srs EPSG:25832 D:\City\raster\gas\500\g8746-1_cit.tif D:\City\raster\gas\500\g8746-1_cit_epsg_utm.tif i get following message: Copying color table from D:\City\raster\gas\500\g8746-1_cit.tif to new file. Creating output file that is 0P x 0L. ERROR 1: Attempt to create 0x0 dataset is illegal,sizes must be larger than zero . did anyone / you have a idea how to translate successfully? alternative to edit the coordinates manuell in the tif-file? it is very important for me. reagrads Jan :-) ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] A gdal-users mailing list?
Hi all, Could it (perhaps) make sense to have a gdal-users mailing list, in analogy to grass-users, qgis-users, etc. ? I assume that, according to some logic, GDAL/OGR is considered to be a library from/for developers, not a GIS application for end-users. On the other hand, there are quite some documented utilities around [1], utilised by gdal-users like me. Actually, if I compare [1] with [2], I see that my make install installs many more utilities/scripts than mentioned on the GDAL Utilities page. Users might have (..hmm..) simple questions that real developers might consider to be (..hmm..) not so interesting. gdal-users could be a list where experienced users reply to newcomers, with the side-effect of reducing some traffic on gdal-dev. Hermann [1] http://www.gdal.org/gdal_utilities.html [2] $ make install 1 libtool: install: (...) gdalinfo 2 libtool: install: (...) gdal_translate 3 libtool: install: (...) gdaladdo 4 libtool: install: (...) gdalwarp 5 libtool: install: (...) nearblack 6 libtool: install: (...) gdalmanage 7 libtool: install: (...) gdalenhance 8 libtool: install: (...) gdaltransform 9 libtool: install: (...) gdaldem 10 libtool: install: (...) gdallocationinfo 11 libtool: install: (...) gdalsrsinfo 12 libtool: install: (...) gdal_contour 13 libtool: install: (...) gdaltindex 14 libtool: install: (...) gdal_rasterize 15 libtool: install: (...) gdal_grid 16 libtool: install: (...) ogrinfo 17 libtool: install: (...) ogr2ogr 18 libtool: install: (...) ogrtindex 19 libtool: install: (...) testepsg 20 libtool: install: (...) gdalbuildvrt 21 libtool: install: (...) scripts/epsg_tr.py 22 libtool: install: (...) scripts/esri2wkt.py 23 libtool: install: (...) scripts/gcps2vec.py 24 libtool: install: (...) scripts/gcps2wld.py 25 libtool: install: (...) scripts/gdal2tiles.py 26 libtool: install: (...) scripts/gdal2xyz.py 27 libtool: install: (...) scripts/gdal_calc.py 28 libtool: install: (...) scripts/gdalchksum.py 29 libtool: install: (...) scripts/gdal_fillnodata.py 30 libtool: install: (...) scripts/gdalident.py 31 libtool: install: (...) scripts/gdalimport.py 32 libtool: install: (...) scripts/gdal_merge.py 33 libtool: install: (...) scripts/gdal_polygonize.py 34 libtool: install: (...) scripts/gdal_proximity.py 35 libtool: install: (...) scripts/gdal_retile.py 36 libtool: install: (...) scripts/gdal_sieve.py 37 libtool: install: (...) scripts/mkgraticule.py 38 libtool: install: (...) scripts/pct2rgb.py 39 libtool: install: (...) scripts/rgb2pct.py ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: Strange behavior of gdalwarp with -overwrite option
On 08/11/2011 20:28, Even Rouault wrote: gdalwarp something.tif something.tif -t_srs EPSG: makes no sense... Indeed: typing 'something.tif' twice should not be necessary. But the following could theoretically make sense (in analogy to GNU sed's -i/--in-place switch (which creates a temporary file and later *overwrites* the original file): gdalwarp -i something.tif -t_srs EPSG: However, as one hardly get things right in the first run: there is a good chance that the result is not what you want -- and the original file is gone :-(. At least this happens to me whenever I use sed's -i switch. So I wouldn't see a need for introducing this behaviour in GDAL command line utilitiies. Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: Using a gdal generated tiff to generate DEM tiles - losing registration
On 30/10/2011 20:16, Iain wrote: But the tiles I generate appear to be out about 1/4 degree too far north (look OK east west) over the whole area. (the mapnik generated tiles, using background from hill shade + colour-relief and data from openstreetmap in a postgis database are all fine and match exactly the position of tiles from the openstreetmap server. - OSM tiles come in a *Spherical* Mercator projection - what you observe seems to be a latitude-only shift This type of datum shift often occurs when projecting from a sphere to an ellpsoid or vice-versa, although 1/4 degree (~28 km) sounds like quite a lot. Here 2 links to recent gdal-dev postings on this issue: http://osgeo-org.1803224.n2.nabble.com/gdal-dev-Ogr2ogr-projection-problem-features-are-shifted-12-miles-to-the-north-How-do-I-fix-it-td6779672.html http://osgeo-org.1803224.n2.nabble.com/gdal-dev-Ogr2ogr-datum-shift-td6795580.html ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: 254 into 255
On 27/10/2011 14:48, Chaitanya kumar CH wrote:
However it [gdalbuildvrt] doesn't support complex sources..
Chaitanya, could you perhaps elaborate on the above statement?
joolek,
About the batch-processing of 3000 files: in my
Linux/Bash/GDAL-from-trunk environment, I would do something like
indicated below.
Hope this helps, Hermann
for f in path/to/mytiffs/*.tif; do
gdal_translate -of vrt $f /dev/stdout |
awk '
BEGIN { lut = LUT0:0,253:253,254:255/LUT }
{ gsub(/SimpleSource/, ComplexSource) }
$1 == /ComplexSource { $0 = lut ORS $0 }
{ print }' tmp.vrt
gdal_translate tmp.vrt ${f%.tif}_modified.tif
done
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: 254 into 255
On 27/10/2011 10:58, Chaitanya kumar CH wrote: A vrt file can implement your requirements. It can create a 'lookup table' to translate your pixel values... Chaitanya, I just tried to translate an Int32 GeoTIFF with pixel values 111..523 into a byte-type GeoTIFF with values 1..44, using the LUT as indicated below. I noted that in the resulting GeoTIFF, the previous NODATA areas have now a pixel value of 1. Is this expected behaviour or rather some sort of over-extrapolation of my lookup table (as I would assume) ? Thanks in advance, Hermann VRTDataset rasterXSize=8786 rasterYSize=6443 ... VRTRasterBand dataType=Int32 band=1 NoDataValue0.00E+00/NoDataValue ColorInterpGray/ColorInterp ComplexSource SourceFilename relativeToVRT=1ch00_form.tif/SourceFilename SourceBand1/SourceBand SourceProperties RasterXSize=8786 RasterYSize=6443 DataType=Int32 BlockXSize=256 BlockYSize=256 / SrcRect xOff=0 yOff=0 xSize=8786 ySize=6443 / DstRect xOff=0 yOff=0 xSize=8786 ySize=6443 / LUT111:1,112:2,121:3,122:4,123:5,124:6,131:7,132:8,.../LUT /ComplexSource /VRTRasterBand /VRTDataset ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: 254 into 255
On 28/10/2011 20:02, Chaitanya kumar CH wrote: What was the nodata value in the original? NoDataValue0.00E+00/NoDataValue Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: 254 into 255
On 28/10/2011 20:08, Chaitanya kumar CH wrote: I meant the nodata value of ch00_form.tif This *is* the nodata value of ch00_form.tif. The vrt which I posted has been generated through a simple: $ gdal_translate -of vrt ch00_form.tif ch00_form.vrt The only modifications I made was to change SimpleSource into ComplexSource, and to add the LUT. Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Another PGeo driver issue, on 64-bit Debian Linux
Hi, I have been trying to get the PGeo driver working, following the hints at [0]. I have gotten to a partial success in so far that the driver connects to the mdb and finds the layers and the features [1]. However, the geometry is not understood, which seems to be related to the importFromWKT() failed error message, which is repeated for each layer [2]. Querying individual layers gives an additional GetFeatureCount() failed error, see under [3]. I searched on Internet and found a few (older) postings mentioning similar errors where the WKT string is always truncated to 12 characters, but nothing could really help me fixing the issue. Does anyone know how to fix the issue? (I am using GDAL from trunk, on 64-bit Debian Linux, unixodbc 2.2.14p2-4, mdbtools 0.5.99.0.6pre1.) Thanks in advance, Hermann [0] http://gdal.org/ogr/drv_pgeo.html [1] $ ogrinfo --debug on pgeo:mt PGeo: MDB Tools driver: /usr/lib/libmdbodbc.so.0 PGeo: MDB Tools driver installed successfully! PGeo: EstablishSession(mt) ODBC: SQLConnect(mt) INFO: Open of `pgeo:mt' using driver `PGeo' successful. OGR: GetLayerCount() = 26 1: CoastL 2: DamL 3: WatrcrsL ... [2] ERROR 1: importFromWKT() failed on SRS 'GEOGCS[GCS_'. [3] $ ogrinfo --debug on PGeo:mt CoastL PGeo: MDB Tools driver: /usr/lib/libmdbodbc.so.0 PGeo: MDB Tools driver installed successfully! PGeo: EstablishSession(mt) ODBC: SQLConnect(mt) OGR: OGROpen(PGeo:mt/0x8804b0) succeeded as PGeo. INFO: Open of `PGeo:mt' using driver `PGeo' successful. OGR: GetLayerCount() = 26 Layer name: CoastL Geometry: Unknown (any) Error at Line : syntax error near ( ERROR 1: GetFeatureCount() failed on query SELECT COUNT(*) FROM CoastL. Feature Count: 11 Extent: (14.183600, 35.806450) - (14.576150, 36.082430) Layer SRS WKT: (unknown) FID Column = OBJECTID Geometry Column = Shape OBJECTID: Integer (0.0) FCsubtype: Integer (0.0) gfid: String (0.0) F_CODE: String (0.0) ICC: String (0.0) SN: Integer (0.0) Shape_Length: Real (0.0) PGeo: 11 features read on layer 'CoastL'. ODBC: SQLDisconnect() ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Another PGeo driver issue, on 64-bit Debian Linux
On 26/10/2011 18:11, Even Rouault wrote: IMHO, you can't. I've spent some time investigating this issues a few months ago and my conclusion was that mdbtools was too buggy (and unmaintained), and in particular not 64bit ready. I tried a few fixes, but finally I gave up. I've developped the MDB driver that relies on the Jackess Java library to workaround mdbtools problems and it works quite well (although not particularly fast due to the invokation of Java code from native code). Thanks for the quick response. If you didn't get it to work, then I should better not have tried in the first place. Perhaps some warnings for Linux users should be added to http://gdal.org/ogr/drv_pgeo.html (?) I'll have a look at the MBD driver. Thanks again, Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: OGR FileGDB driver: Failed at creating table ... (General function failure.)
On 16/09/2011 16:49, Paul Ramsey wrote: ...perhaps in the FGDB driver we can try and avoid using the WKT at all when we have a WKID available. I just disabled the generation of the WKT element in my local copy of FGdbLayer.cpp [0]: /* CPLCreateXMLElementAndValue(srs_xml,WKT, wkt); */ If I now run: ogr2ogr -f filegdb out.gdb in.gdb -a_srs epsg:3035, the generated XML looks like [1], i.e. the only meaningful information about the ProjectedCoordinateSystem is: WKID3035/WKID. out.gdb is created without any error message and ogrinfo tells me that its SRS is EPSG:3035, see [2]. However, ArcCatalog 10 considers the same out.gdb to have an unknown spatial reference :-( Hermann [0] http://trac.osgeo.org/gdal/browser/trunk/gdal/ogr/ogrsf_frmts/filegdb/FGdbLayer.cpp#L371 [1] $ ogr2ogr -f filegdb out.gdb in.gdb -a_srs epsg:3035 ... SpatialReference xsi:type=esri:ProjectedCoordinateSystem XOrigin-2147483647/XOrigin YOrigin-2147483647/YOrigin XYScale10/XYScale ZOrigin-2147483647/ZOrigin ZScale10/ZScale XYTolerance0.0001/XYTolerance ZTolerance0.0001/ZTolerance HighPrecisiontrue/HighPrecision WKID3035/WKID /SpatialReference ... [2] $ ogrinfo -al -so out.gdb OGR: OGROpen(out.gdb/0x10d7140) succeeded as FileGDB. INFO: Open of `out.gdb' using driver `FileGDB' successful. Layer name: eea_1Kgrid Geometry: Multi Polygon Feature Count: 9 Extent: (-17.00, 28.00) - (-17.00, 28.00) Layer SRS WKT: PROJCS[ETRS 1989 / LAEA, GEOGCS[ETRS89, DATUM[European_Terrestrial_Reference_System_1989, SPHEROID[GRS 1980,6378137,298.257222101, AUTHORITY[EPSG,7019]], TOWGS84[0,0,0,0,0,0,0], AUTHORITY[EPSG,6258]], PRIMEM[Greenwich,0, AUTHORITY[EPSG,8901]], UNIT[degree,0.0174532925199433, AUTHORITY[EPSG,9122]], AUTHORITY[EPSG,4258]], PROJECTION[Lambert_Azimuthal_Equal_Area], PARAMETER[latitude_of_center,52], PARAMETER[longitude_of_center,10], PARAMETER[false_easting,4321000], PARAMETER[false_northing,321], UNIT[metre,1, AUTHORITY[EPSG,9001]], AXIS[Y,NORTH], AXIS[X,EAST], AUTHORITY[EPSG,3035]] FID Column = OBJECTID Geometry Column = SHAPE CELLCODE: String (0.0) Shape_Length: Real (0.0) Shape_Area: Real (0.0) ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Re: OGR FileGDB driver: Failed at creating table ... (General function failure.)
On 19/09/2011 18:41, Paul Ramsey wrote: Wow, we just can't win can we? So the only way forward is to get morphToEsri to produce the extract literal string representation expected by the API? What fun. It is really a pity: over some 500 lines of code, OGR's morphToESRI() function tries its best to generate an ESRI-friendly WKT representation of EPSG:3035 and comes pretty close.. but eventually, the API aborts layer creation because the generated WKT string starts with ETRS89_LAEA_Europe, whereas ETRS_1989_LAEA is expected. OTOH: I can imagine that amending the code so that an exact string match will be achieved for most SRSs would really be a nightmare. Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: OGR FileGDB driver: Failed at creating table ... (General function failure.)
On 15/09/2011 16:46, Hermann Peifer wrote: My quick-fix: I edited ../local/share/gdal/pcs.csv and changed ETRS89 / LAEA Europe into ETRS 1989 / LAEA, ..et voilà: the error is gone. ETRS 1989 / LAEA is ESRI'ified into ETRS_1989_LAEA, which is exactly the expected string for this SRS. I am not sure if it is worth filing a ticket for this specific case, as it is part of a wider ESRI-WKT-versus-GDAL/OGR-WKT dilemma, which in return is known already, see e.g. Uwe's post from some days ago. I am not sure if people are still interested but for the sake of the exercise I just created a new FileGDB by assigning a spatial reference which only contains ESRI's string identifier of EPSG:3035, see [1]. This works fine and does not generate any error. However, any minor change in the string ETRS_1989_LAEA leads to the FileGDB API's General function failure. Frank, Even, Paul, et al.: Could it perhaps make sense if the morphToESRI() function paid more attention to these string identifiers in order to produce even more ESRI-friendly WKT? Hermann [1] $ ogr2ogr -f filegdb out.gdb in.gdb -a_srs nonsense.wkt $ cat nonsense.wkt PROJCS[ETRS_1989_LAEA] $ ogrinfo -al -so out.gdb OGR: OGROpen(out.gdb/0x7a2240) succeeded as FileGDB. INFO: Open of `out.gdb' using driver `FileGDB' successful. Layer name: eea_1Kgrid Geometry: Multi Polygon Feature Count: 9 Extent: (-17.00, 28.00) - (-17.00, 28.00) Layer SRS WKT: PROJCS[ETRS_1989_LAEA, UNIT[Meter,1]] FID Column = OBJECTID Geometry Column = SHAPE CELLCODE: String (0.0) Shape_Length: Real (0.0) Shape_Area: Real (0.0) ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Re: OGR FileGDB driver: Failed at creating table ... (General function failure.)
On 16/09/2011 16:49, Paul Ramsey wrote: On Fri, Sep 16, 2011 at 3:16 AM, Hermann Peiferpei...@gmx.eu wrote: Frank, Even, Paul, et al.: Could it perhaps make sense if the morphToESRI() function paid more attention to these string identifiers in order to produce even more ESRI-friendly WKT? Yes, probably we will have to. I wonder how we would go about figuring out the relevant differences...? In the meanwhile, since you've found that using the WKID works perhaps in the FGDB driver we can try and avoid using the WKT at all when we have a WKID available. P Isn't the WKID a numeric identifier, and not the string identifier I was mentioning in my previous mail (?) My little exercise only showed that with an appropriate string in the WKT (here: ETRS_1989_LAEA), the General function failure can be avoided when creating out.gdb through the FileGDB driver. It doesn't mean that the string value ETRS_1989_LAEA is enough to define a valid spatial reference. When opening the generated out.gdb in ArcCatalog, I get a General function failure, Cannot read SRS... failure. AFAICT: in order to create a meaningful spatial reference, the complete WKT has to be passed over: morphed string identifiers, morphed projection and parameter names. etc. Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: OGR FileGDB driver: Failed at creating table ... (General function failure.)
AFAICT, the issue is not so much that the FileGDB only supports certain projections, but that it only accepts a very specific WKT string for a given SRS. A single character of difference can cause the FileGDB API to not recognise the SRS, then fall back to a General function failure. GDAL/OGR has its own WKT definition of a given SRS, which is slightly different to ESRI's WKT of the same SRS. I noted that when creating a FileGDB, GDAL/OGR tries to ESRI'ify its WKT through the morphToESRI() function. For EPSG:3035, the result is obviously not ESRI-friendly enough. In fact, it is not what I first thought: an issue with the parameter name (Central_Meridian versus longitude_of_center) but simply the SRS's string identifier which causes the problem. My quick-fix: I edited ../local/share/gdal/pcs.csv and changed ETRS89 / LAEA Europe into ETRS 1989 / LAEA, ..et voilà: the error is gone. ETRS 1989 / LAEA is ESRI'ified into ETRS_1989_LAEA, which is exactly the expected string for this SRS. I am not sure if it is worth filing a ticket for this specific case, as it is part of a wider ESRI-WKT-versus-GDAL/OGR-WKT dilemma, which in return is known already, see e.g. Uwe's post from some days ago. Hermann On 15/09/2011 13:28, Smith, Michael ERDC-CRREL-NH wrote: I've run into this. It seems to be that there are only certain projections that are supported in FileGDB. I tried to convert data in to the Spherical Mercator projection and none of the projections I used (either epsg, esri epsg, or esri wkt) would work. I believe that Even R. determined that it was a limitation of the FileGDB API currently. A ticket should be filed so that this can be documented. Mike ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] OGR FileGDB driver: ERROR 1: 'OBJECTID' not recognised as an available field.
Hi All, I am not quite sure if this is a feature or rather a bug, but OGR's -where switch seems to be aware of the OBJECTID field in a FileGDB, whereas the -sql option reports: 'OBJECTID' not recognised as an available field, see below. Would it be worth filing a ticket? Hermann $ ogrinfo -al -so out.gdb INFO: Open of `out.gdb' using driver `FileGDB' successful. Layer name: eea_1Kgrid Geometry: Multi Polygon Feature Count: 9 Extent: (-17.00, 28.00) - (-17.00, 28.00) Layer SRS WKT: GEOGCS[WGS 84, DATUM[WGS_1984, SPHEROID[WGS 84,6378137,298.257223563, AUTHORITY[EPSG,7030]], AUTHORITY[EPSG,6326]], PRIMEM[Greenwich,0, AUTHORITY[EPSG,8901]], UNIT[degree,0.0174532925199433, AUTHORITY[EPSG,9122]], AUTHORITY[EPSG,4326]] FID Column = OBJECTID Geometry Column = SHAPE CELLCODE: String (0.0) Shape_Length: Real (0.0) Shape_Area: Real (0.0) $ $ ogrinfo -q -geom=no out.gdb eea_1Kgrid -where OBJECTID=1 Layer name: eea_1Kgrid OGRFeature(eea_1Kgrid):1 CELLCODE (String) = E1647N1038 Shape_Length (Real) = 0.0383874724895616 Shape_Area (Real) = 9.18055152110593e-05 $ ogrinfo -q -geom=no out.gdb -sql select * from eea_1Kgrid where OBJECTID=1 ERROR 1: 'OBJECTID' not recognised as an available field. $ ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] OGR FileGDB driver: Failed at creating table ... (General function failure.)
Hi again, Once I am at it, I could also report another issue with the FileGDB driver. The driver re-projects my test geodatabase from WGS84 to NAD83 without any problems, but when trying to re-project to EPSG:3035 (ETRS89/LAEA Europe), I end up with a General function failure, see more details below. I am not sure what the actual issue is, but a quick Internet search suggests that the FileGDB API's General function failure might be related to problems with understanding the SRS. Maybe OGR's SRS definition is not ESRI-friendly enough when trying to create the geodatabase(?). According to my experience, the classic problem with the LAEA projection is the parameter name: Central_Meridian(ESRI) versus longitude_of_center (GDAL/OGR). Would it be worth filing a ticket? Hermann $ ogrinfo -al -so out.gdb INFO: Open of `out.gdb' using driver `FileGDB' successful. Layer name: eea_1Kgrid Geometry: Multi Polygon Feature Count: 9 Extent: (-17.00, 28.00) - (-17.00, 28.00) Layer SRS WKT: GEOGCS[WGS 84, DATUM[WGS_1984, SPHEROID[WGS 84,6378137,298.257223563, AUTHORITY[EPSG,7030]], AUTHORITY[EPSG,6326]], PRIMEM[Greenwich,0, AUTHORITY[EPSG,8901]], UNIT[degree,0.0174532925199433, AUTHORITY[EPSG,9122]], AUTHORITY[EPSG,4326]] FID Column = OBJECTID Geometry Column = SHAPE CELLCODE: String (0.0) Shape_Length: Real (0.0) Shape_Area: Real (0.0) $ $ $ ogr2ogr -f filegdb out_4269.gdb out.gdb -t_srs epsg:4269 $ ogrinfo -al -so out_4269.gdb INFO: Open of `out_4269.gdb' using driver `FileGDB' successful. Layer name: eea_1Kgrid Geometry: Multi Polygon Feature Count: 9 Extent: (-17.00, 28.00) - (-17.00, 28.00) Layer SRS WKT: GEOGCS[NAD83, DATUM[North_American_Datum_1983, SPHEROID[GRS 1980,6378137,298.257222101, AUTHORITY[EPSG,7019]], TOWGS84[0,0,0,0,0,0,0], AUTHORITY[EPSG,6269]], PRIMEM[Greenwich,0, AUTHORITY[EPSG,8901]], UNIT[degree,0.0174532925199433, AUTHORITY[EPSG,9122]], AUTHORITY[EPSG,4269]] FID Column = OBJECTID Geometry Column = SHAPE CELLCODE: String (0.0) Shape_Length: Real (0.0) Shape_Area: Real (0.0) $ $ ogr2ogr -f filegdb out_3035.gdb out.gdb -t_srs epsg:3035 ERROR 1: Error: Failed at creating table for \eea_1Kgrid (General function failure.) ERROR 1: Terminating translation prematurely after failed translation of layer eea_1Kgrid (use -skipfailures to skip errors) ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: AW: ESRI products have problems reading gdal spatialreference entry
On 06/09/2011 17:39, Frank Warmerdam wrote: On 11-09-06 12:54 AM, Schmitz, Uwe wrote: So if I conclude: ESRI and GDAL have different WKT formats. It is not possible to write GeoTIFF files which are fully ESRI *and* GDAL compatible. Uwe, I do not agree with this conclusion! However, it is difficult to produce a GeoTIFF file that will exactly preserve and ESRI Project Engine string (ie. ESRI WKT) while at the same time being fully compatible with all other software. The issues are primarily around limits of the GeoTIFF format, and to some extend somewhat incomplete code implementations in GDAL / ArcGIS. I would really be interested in such an ESRI-friendly GeoTIFF file. Could it be that in theory it is just difficult to produce such a GeoTIFF (as Frank writes above) whereas in practice and as of today, it is actually impossible (as Uwe wrote)? ;-) Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: ESRI products have problems reading gdal spatialreference entry
On 05/09/2011 16:03, Schmitz, Uwe wrote: I wonder, if anyone else has experienced this or similar behaviour. And how can I achieve that ESRI- *and* gdal-Tools can identify the correct reference system. What I asked on this list some weeks ago: --- snip --- From experience, I know that ESRI users have trouble with some projection and parameter names generated by GDAL and vice versa. Is it possible to convince GDAL to store coordinate system information in a way that would be ESRI-friendly and GDAL-friendly at the same time? http://osgeo-org.1803224.n2.nabble.com/gdal-dev-gdal-translate-eats-projection-parameter-values-td6680528.html#d1313170021000-317 --- snip --- My understanding of Frank's answer is that this is basically not possible. Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: ogr2ogr - output WKT to console?
On 31/08/2011 23:28, Even Rouault wrote: ... But it was easy to also support /dev/stdout as an alias for /vsistdout/. Added in r23016. Works as expected. Thanks for your time. Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] OGR GML and other drivers producing text output: stdout, /dev/stdout and /vsistdout/ (was: ogr2ogr - output WKT to console?)
On 01/09/2011 00:28, Even Rouault wrote: Le jeudi 01 septembre 2011 00:19:52, Elijah Robison a écrit : ogr2ogr -f GML /vsistdout/ C:\xData\CountyUpdate\GeoFix\feb16_polygons.shp -sql select * from feb16_polygons where objectid= 2126 ogr2ogr -f KML /vsistdout/ C:\xData\CountyUpdate\GeoFix\feb16_polygons.shp -sql select * from feb16_polygons where objectid= 2126 ogr2ogr -f GeoJSON /vsistdout/ C:\xData\CountyUpdate\GeoFix\feb16_polygons.shp -sql select * from feb16_polygons where objectid= 2126 You can also add the BNA, PGDump, GeoRSS and GPX drivers to the list About the GML driver: I was just wondering where the .xsd file is written to when using: ogr2ogr -f gml /vsistdout/. According to the gml file's xsi:schemaLocation attribute, the schema file is `.xsd', see [1]. I can't find `.xsd' locally. When using /dev/stdout rather than /vsistdout/, I receive an error, which tells me that the schema file can't be created, see [2]. Another observation with the gml driver: stdout is defined as an alias for /vsistdout/ (rather than /dev/stdout). When using the alias, xsi:schemaLocation tells me that the schema file is `stdout.xsd', see [3]. I can't find `stdout.xsd' locally. About other drivers producing text output: There seem to be several options about how to alias (/dev/)stdout and /vsistdout/. The current practice is perhaps somewhat inconsistent, see [4]. IMO, Even's recently added /dev/stdout - /vsistdout/ alias looks most appropriate, I am however not so sure about non-Unix OSs. Finally, only drv_pgdump.html seems to mention /vsistdout/ as an option for writing to standard output, see [5]. Perhaps other relevant driver pages should also mention this option ? Would it be worth creating a ticket (or two) related to the above observations? Hermann [1] $ ogr2ogr -f gml /vsistdout/ test.shp ?xml version=1.0 encoding=utf-8 ? ogr:FeatureCollection xmlns:xsi=http://www.w3.org/2001/XMLSchema-instance; xsi:schemaLocation=http://ogr.maptools.org/ .xsd xmlns:ogr=http://ogr.maptools.org/; xmlns:gml=http://www.opengis.net/gml; gml:featureMember ogr:test fid=F0 (...) [2] ERROR 4: Failed to open file /dev/stdout.xsd for schema output. [3] $ ogr2ogr -f gml stdout test.shp ?xml version=1.0 encoding=utf-8 ? ogr:FeatureCollection xmlns:xsi=http://www.w3.org/2001/XMLSchema-instance; xsi:schemaLocation=http://ogr.maptools.org/ stdout.xsd xmlns:ogr=http://ogr.maptools.org/; xmlns:gml=http://www.opengis.net/gml; gml:featureMember ogr:test fid=F0 (...) [4] ./gpx/ogrgpxdatasource.cpp:if( EQUAL(pszFilename,stdout) || EQUAL(pszFilename,/vsistdout/)) ./gml/ogrgmldatasource.cpp:if( EQUAL(pszFilename,stdout) || EQUAL(pszFilename,/vsistdout/)) ./kml/ogrkmldatasource.cpp:if( EQUAL(pszName, stdout) || EQUAL(pszName, /vsistdout/) ) ./csv/ogrcsvdriver.cpp:if (strcmp(pszName, /dev/stdout) == 0) ./bna/ogrbnadatasource.cpp:if( EQUAL(pszFilename,stdout) ) ./georss/ogrgeorssdatasource.cpp:if( EQUAL(pszFilename,stdout) ) ./geojson/ogrgeojsondatasource.cpp:if( EQUAL( pszName, stdout ) ) [5] $ find ogr/ -name '*.html' | xargs grep /vsistdout/ ogr/ogrsf_frmts/pgdump/drv_pgdump.html:to specify /vsistdout/ as output file to output on the standard output.p ogr/ogrsf_frmts/pgdump/drv_pgdump.html:% ogr2ogr --config PG_USE_COPY YES -f PGDump /vsistdout/ abc.shp | psql -d my_dbname -f - ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: OGR GML and other drivers producing text output: stdout, /dev/stdout and /vsistdout/
On 01/09/2011 13:33, Even Rouault wrote: But yes, it lacks the special case to recognize /dev/stdout as a particular value where you must not try to write a file, which would avoid error [2]. Ticket created at http://trac.osgeo.org/gdal/ticket/4225 Yes we probably lack consistency in that area. I think the stdout case should be removed (should not affect many people as it was unadvertized feature), and only accept /dev/stdout and /vsistdout/ which are less hacky solutions. Ticket created: http://trac.osgeo.org/gdal/ticket/4226 Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: ogr2ogr - output WKT to console?
On 31/08/2011 22:52, Even Rouault wrote: Yes, as a output filename, you can use the special name /vsistdout/ (you need the trailing slash) . You likely need GDAL 1.8 or later for it to work. Even, Thanks for the hint. Is there a specific reason why the csv driver isn't able to use /dev/stdout, similar to the kml driver [1] ? Hermann [1] $ ogr2ogr -f kml /dev/stdout test.shp ?xml version=1.0 encoding=utf-8 ? kml xmlns=http://www.opengis.net/kml/2.2; DocumentFoldernametest/name Schema name=test id=test SimpleField name=Name type=string/SimpleField SimpleField name=Description type=string/SimpleField SimpleField name=Latitude type=string/SimpleField SimpleField name=Longitude type=string/SimpleField SimpleField name=Value type=string/SimpleField /Schema Placemark nameFirst point/name ExtendedDataSchemaData schemaUrl=#test SimpleData name=NameFirst point/SimpleData SimpleData name=Latitude50/SimpleData SimpleData name=Longitude50/SimpleData SimpleData name=Value1/SimpleData /SchemaData/ExtendedData Pointcoordinates50,50/coordinates/Point /Placemark (...) ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Re: gdal_rasterize output off by one pixel/cell
On 28/08/2011 14:25, Even Rouault wrote: Le samedi 27 août 2011 11:14:00, Hermann Peifer a écrit : Matthew, An absolute kludge-type of solution would be a) to make sure that the GeoTIFF is larger than the shapefile when rasterising (in order to not loose points at the edges) and b) to shift the GeoTIFF NW by 0.5 pixel after rasterising, as I already wrote in my previous mail. There has been no traffic on Eli's ticket http://trac.osgeo.org/gdal/ticket/3774 during the last 11 months, so my guess is that the bug will remain for a while. It's strange that there aren't more people stumbling across the issue. I've applied a fix in trunk in http://trac.osgeo.org/gdal/changeset/22998 . Please test and confirm if it works for you. Even, After some testing, I can confirm that both issues (half-pixel shift and lost point at edges) are gone. This is great. Thanks for spending your time fixing this issue, on a sunny Sunday afternoon. Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: gdal_rasterize output off by one pixel/cell
Matthew, An absolute kludge-type of solution would be a) to make sure that the GeoTIFF is larger than the shapefile when rasterising (in order to not loose points at the edges) and b) to shift the GeoTIFF NW by 0.5 pixel after rasterising, as I already wrote in my previous mail. There has been no traffic on Eli's ticket http://trac.osgeo.org/gdal/ticket/3774 during the last 11 months, so my guess is that the bug will remain for a while. It's strange that there aren't more people stumbling across the issue. Hermann On 26/08/2011 22:07, Burton-Kelly, Matthew wrote: Eli, I finally got around to looking at the example shapefile and running the same operation on my machine, and it does indeed seem like the same bug. Thanks for the good eye, this almost seems like two separate issues; first being the missing edge points and second being the weird offset. I only say this because the offset persists even when I set the extents of the output raster as larger than the input layer. Matt On 22 Aug, 2011, at 11:50 AM, Eli Adam wrote: Matt, You may want to look at this ticket and see if it is the same thing. If so, you can add yourself to the cc list and you will get emails when there are updates to the ticket. You can also add any additional relevant information to the ticket. http://trac.osgeo.org/gdal/ticket/3774 Regards, Eli On 8/20/2011 at 11:06 AM, in message aa52f8ac-89dc-4543-a210-0dd0ba9e5...@my.und.edu, Burton-Kelly, Matthew matthew.burton.ke...@my.und.edu wrote: Hello, I'm attempting to create a raster from a shapefile of point data, with grid cells 1 degree square. Although the area I have defined matches up with the origin coordinates I want, and the grid cells match up with a vector 1-degree grid I produced in QGIS, the squares supposedly containing the point data are shifted, it looks like down and to the right. Has anyone encountered this issue before? Here is the command I used: gdal_rasterize -acolumn name -ts 360 180 -te -180 -90 180 90 -lsource layer name in source file source file destination file I have uploaded a screenshot of the output I am describing here: http://www.flickr.com/photos/matthewbk/6059112633 Thanks for any insight, Matt ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] OGR csv driver: feature numbering
Hi, I am just wondering why the csv driver starts numbering at 1, rather than at 0, see e.g. the First point in the test.csv example at [1], which ends up as OGRFeature(test):1. However, if I take the test.vrt/test.csv example and convert it to shapefile format, I end up with features 0, 1 and 2, see [2]. Is this how it should be? Hermann [1] http://www.gdal.org/ogr/drv_csv.html [2] $ ogr2ogr test.shp test.vrt $ ogrinfo -al test.shp INFO: Open of `test.shp' using driver `ESRI Shapefile' successful. Layer name: test Geometry: Point Feature Count: 3 Extent: (0.25, 47.50) - (1.10, 49.20) Layer SRS WKT: GEOGCS[GCS_WGS_1984, DATUM[WGS_1984, SPHEROID[WGS_84,6378137,298.257223563]], PRIMEM[Greenwich,0], UNIT[Degree,0.017453292519943295]] Latitude: String (80.0) Longitude: String (80.0) Name: String (80.0) OGRFeature(test):0 Latitude (String) = 48.1 Longitude (String) = 0.25 Name (String) = First point POINT (0.25 48.1) OGRFeature(test):1 Latitude (String) = 49.2 Longitude (String) = 1.1 Name (String) = Second point POINT (1.1 49.2) OGRFeature(test):2 Latitude (String) = 47.5 Longitude (String) = 0.75 Name (String) = Third point POINT (0.75 47.5) ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] OGR csv driver: feature numbering
On 27/08/2011 12:25, Even Rouault wrote: Le samedi 27 août 2011 12:07:12, Hermann Peifer a écrit : Hi, I am just wondering why the csv driver starts numbering at 1, rather than at 0, see e.g. the First point in the test.csv example at [1], which ends up as OGRFeature(test):1. However, if I take the test.vrt/test.csv example and convert it to shapefile format, I end up with features 0, 1 and 2, see [2]. Is this how it should be? AFAIK, there is no convention imposed by the OGR model on how numbering should be done and the current situation depends much of the mood of the driver authors when they coded it. I want to draw your attention that, pedantically speaking, this number is not a sequence number, but a feature ID, with no guarantee that it is sequential when iterating over the features. For shapefiles or database drivers, you can have holes in the numbering due to deleted features. For database drivers when issuing a SQL request, the features can also even be returned in random order w.r.t to their feature ids. But I agree it would have been good to behave consistently in the case of no holes. Not sure if it is worth changing that now and potentially causing issues to people (wrongly) relying on how the numbering is currently done. Even, Thanks once again for your quick reply. I just came across when trying to get the first point from a csv file into a shapefile, using something like: $ ogr2ogr test.shp test.vrt -where fid=0 The resulting shapefile was empty, to my surprise. I see that the kml driver also starts numbering at 1. I do agree that it might not be worth changing the current, somewhat inconsistent behaviour. Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] ERROR 4: `test.asc' not recognised as a supported file format.
Hi, Occasionally I use a small .asc file, which I then gdal_translate into a bigger blank GeoTIFF. I tried to make the smallest possible test.asc file with 1 col x 1 row and 1 cell value only, but I ended up with ERROR 4: file not recognised as a supported file format. The error message seems to go away after adding (useless) values to the file, see below. I am wondering if this is a feature or rather a bug? Thanks for your time, Hermann -- $ gdalinfo --debug on test.asc ERROR 4: `test.asc' not recognised as a supported file format. gdalinfo failed - unable to open 'test.asc'. $ gdalinfo --debug on test2.asc GDAL: GDALOpen(test2.asc, this=0xab0610) succeeds as AAIGrid. Driver: AAIGrid/Arc/Info ASCII Grid Files: test2.asc Size is 1, 1 Coordinate System is `' Origin = (0.000,1.000) Pixel Size = (1.000,-1.000) Corner Coordinates: Upper Left ( 0.000, 1.000) Lower Left ( 0.000, 0.000) Upper Right ( 1.000, 1.000) Lower Right ( 1.000, 0.000) Center ( 0.500, 0.500) Band 1 Block=1x1 Type=Int32, ColorInterp=Undefined GDAL: GDALClose(test2.asc, this=0xab0610) $ gdalinfo --debug on test3.asc GDAL: GDALOpen(test3.asc, this=0x16c5610) succeeds as AAIGrid. Driver: AAIGrid/Arc/Info ASCII Grid Files: test3.asc Size is 1, 1 Coordinate System is `' Origin = (0.000,1.000) Pixel Size = (1.000,-1.000) Corner Coordinates: Upper Left ( 0.000, 1.000) Lower Left ( 0.000, 0.000) Upper Right ( 1.000, 1.000) Lower Right ( 1.000, 0.000) Center ( 0.500, 0.500) Band 1 Block=1x1 Type=Int32, ColorInterp=Undefined GDAL: GDALClose(test3.asc, this=0x16c5610) $ cat test.asc ncols 1 nrows 1 xllcorner 0 yllcorner 0 cellsize 1 0 $ cat test2.asc ncols 1 nrows 1 xllcorner 0 yllcorner 0 cellsize 1 0 0 0 0 0 0 0 0 0 0 $ cat test3.asc ncols 1 nrows 1 xllcorner 0 yllcorner 0 cellsize 1 123456789 987654321 $ ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] ERROR 4: `test.asc' not recognised as a supported file format.
On 26/08/2011 16:33, Even Rouault wrote: Selon Hermann Peiferpei...@gmx.eu: Hi, Occasionally I use a small .asc file, which I then gdal_translate into a bigger blank GeoTIFF. I tried to make the smallest possible test.asc file with 1 col x 1 row and 1 cell value only, but I ended up with ERROR 4: file not recognised as a supported file format. The error message seems to go away after adding (useless) values to the file, see below. I am wondering if this is a feature or rather a bug? Feature/bug/hazard... who knows ;-) The test that identifies if a file is AAIGRID first checks if there are at least 100 bytes in the file, which looks like a rough estimate of the size taken by the minimum header lines and should be fine for any sensible raster. I guess this test could be just killed since the Identify() also checks the presence of one of the keywords. But you have found the workaround, just add enough spaces at the end of your files and it should be OK. Anyway you can create a ticket if you wish Even, Thanks for the explanation. Ticket created at http://trac.osgeo.org/gdal/ticket/4209. For the time being, I just add more spaces into the header lines, like this example, which has just about enough bytes to get through the test: ncols 1 nrows 1 xllcorner 0 yllcorner 0 cellsize 1 0 Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: gdal_rasterize output off by one pixel/cell
Eli, Matt, I just rasterised some point data with GDAL 1.9dev from svn and can confirm both observations: a) points are lost at the edges b) points and pixels do not line up correctly in qgis I can align points and pixels by shifting the corner coordinates of the resulting GeoTIFF by 0.5 pixel, but this is hardly a fix: gdal_translate test.tif test_shifted.tif -a_ullr -180.5 90.5 179.5 -89.5 I was able to reproduce the same results with GDAL 1.6.3, so I am wondering for how long the issue is already lurking around -- or are we just doing something terribly wrong? Hopefully not. Hermann On 22/08/2011 18:50, Eli Adam wrote: Matt, You may want to look at this ticket and see if it is the same thing. If so, you can add yourself to the cc list and you will get emails when there are updates to the ticket. You can also add any additional relevant information to the ticket. http://trac.osgeo.org/gdal/ticket/3774 Regards, Eli On 8/20/2011 at 11:06 AM, in message aa52f8ac-89dc-4543-a210-0dd0ba9e5...@my.und.edu, Burton-Kelly, Matthew matthew.burton.ke...@my.und.edu wrote: Hello, I'm attempting to create a raster from a shapefile of point data, with grid cells 1 degree square. Although the area I have defined matches up with the origin coordinates I want, and the grid cells match up with a vector 1-degree grid I produced in QGIS, the squares supposedly containing the point data are shifted, it looks like down and to the right. Has anyone encountered this issue before? Here is the command I used: gdal_rasterize -acolumn name -ts 360 180 -te -180 -90 180 90 -lsource layer name in source file source file destination file I have uploaded a screenshot of the output I am describing here: http://www.flickr.com/photos/matthewbk/6059112633 Thanks for any insight, Matt ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: Dump Shapefile into DB using ogr2ogr
On 22/08/2011 13:13, Ingo Weinzierl wrote: Hi *, I already tried the approach below - without success. The result of my operation is: The command: ogr2ogr -f ESRI Shapefile -sql SELECT cast(999 as integer(3)) as river_id from SRCShape DESTShape.shp SRCShape.shp The command's result: - ERROR 1: SQL: Unrecognised field name 999 What does ogr2ogr --version say? On my side, it says: GDAL 1.9dev, released 2011/01/18 BTW: you could drop `-f ESRI Shapefile', as it is the default. Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: Dump Shapefile into DB using ogr2ogr
On 22/08/2011 14:08, Ingo Weinzierl wrote: It says GDAL 1.6.3, released 2009/11/19. This version comes along with Debian Squeeze. You have to upgrade to GDAL/OGR = 1.8.0 in order to make use of the new OGR SQL functionality, http://www.gdal.org/ogr/ogr_sql.html Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: Dump Shapefile into DB using ogr2ogr
On 19/08/2011 15:05, Ingo Weinzierl wrote: Thanks for the quick reply, but the -select option only solves one of the two parts of my problem. I still need to append a further column that does not exist in the shapefile and which should be defined manually. I am not sure if I fully understood the issue, but... what about: ... -sql select cast(999 as integer(3)) as field from ... (If you want an additional integer column with name `field', where all values are defined to be 999) Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] gdal_translate eats projection parameter values
Hi, PARAMETER[Central_Meridian,10.0] from the projection file comes out as: PARAMETER[longitude_of_center,0] PARAMETER[Latitude_Of_Origin,52.0] is changed to PARAMETER[latitude_of_center,0] See the details below. Is there something wrong with my projection file or rather with GDAL? Wondering, Hermann PS: I am using gdal from trunk, revision 22757 from 2011-07-19T18:57:35.878941Z PPS: I reported a somewhat similar (?) issue 2 years ago, see http://trac.osgeo.org/gdal/ticket/3016 --- $ cat dummy.asc ncols1 nrows1 xllcorner0.0 yllcorner0.0 cellsize 50.0 NODATA_value - - $ $ cat ETRS_1989_LAEA_L52_M10.prj PROJCS[ETRS_1989_LAEA_L52_M10,GEOGCS[GCS_ETRS_1989,DATUM[D_ETRS_1989,SPHEROID[GRS_1980,6378137.0,298.257222101]],PRIMEM[Greenwich,0.0],UNIT[Degree,0.0174532925199433]],PROJECTION[Lambert_Azimuthal_Equal_Area],PARAMETER[False_Easting,4321000.0],PARAMETER[False_Northing,321.0],PARAMETER[Central_Meridian,10.0],PARAMETER[Latitude_Of_Origin,52.0],UNIT[Meter,1.0]]$ $ $ gdal_translate dummy.asc dummy.tif -a_srs ETRS_1989_LAEA_L52_M10.prj Input file size is 1, 1 0...10...20...30...40...50...60...70...80...90...100 - done. $ $ gdalinfo dummy.tif Driver: GTiff/GeoTIFF Files: dummy.tif Size is 1, 1 Coordinate System is: PROJCS[ETRS_1989_LAEA_L52_M10, GEOGCS[GCS_ETRS_1989, DATUM[D_ETRS_1989, SPHEROID[GRS_1980,6378137,298.257222101]], PRIMEM[Greenwich,0], UNIT[degree,0.0174532925199433]], PROJECTION[Lambert_Azimuthal_Equal_Area], PARAMETER[latitude_of_center,0], PARAMETER[longitude_of_center,0], PARAMETER[false_easting,4321000], PARAMETER[false_northing,321], UNIT[metre,1, AUTHORITY[EPSG,9001]]] Origin = (0.000,50.000) Pixel Size = (50.000,-50.000) Metadata: AREA_OR_POINT=Area Image Structure Metadata: INTERLEAVE=BAND Corner Coordinates: Upper Left ( 0.000, 50.000) ( 43d40' 8.77W, 27d18'43.27S) Lower Left ( 0.000, 0.000) ( 43d40' 9.28W, 27d18'44.80S) Upper Right ( 50.000, 50.000) ( 43d40' 6.89W, 27d18'43.45S) Lower Right ( 50.000, 0.000) ( 43d40' 7.41W, 27d18'44.97S) Center ( 25.000, 25.000) ( 43d40' 8.09W, 27d18'44.12S) Band 1 Block=1x1 Type=Int32, ColorInterp=Gray NoData Value=- ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] gdal_translate eats projection parameter values
On 12/08/2011 17:37, Even Rouault wrote: Selon Hermann Peiferpei...@gmx.eu: In fact, it is not a bug, but just one of those annoying subtelties you encounter when playing with projections. (...) So : gdal_translate dummy.asc dummy.tif -a_srs ESRI::ETRS_1989_LAEA_L52_M10.prj Or you could just copy ETRS_1989_LAEA_L52_M10.prj as dummy.prj, because the AAIGRID driver will automatically consider that a .prj file is a ESRI one. Best regards, Even Frank, Even, Thanks for the quick response, I will do as you suggested. I wasn't aware of the ESRI:: prefix and wrongly assumed that GDAL would always auto-detect ESRI-isms in the projection files and morph them. The background of the issue is that we disseminate a series of GDAL-generated GeoTIFFs (all in EPSG:3035) via a public website, http://www.eea.europa.eu/data-and-maps/data#c17=raster+data From experience, I know that ESRI users have trouble with some projection and parameter names generated by GDAL and vice versa. Is it possible to convince GDAL to store coordinate system information in a way that would be ESRI-friendly and GDAL-friendly at the same time (and follow the OGC-WKT standard, as far as it exists ... I can't find it at http://www.opengeospatial.org/standards) ? Thanks again for your time, Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: gdalinfo: corner coords in other projection
On 26/05/2011 09:21, Jukka Rahkonen wrote: Lefman, Jonathan ERDC-TEC-VAJonathan.Lefmanat usace.army.mil writes: Hi all, Is there an executable utility that takes the corner coordinates from a gtiff and translates them into another projection system? For example, I have a gtiff in UTM and I want to translate the coordinates of the bounding box (or corners) into WGS84. Or is there a way to do this non-programmatically using a combination of executable utilities? The programmers should stop reading now, this is a farmer's son recipe for doing what you want with something that already exists. - Use gdaltindex for sending the bounding box of your image into a shapefile - Use ogr2ogr with s_srs and t_srs and reproject the shapefile into the system you are interested in. - Use ogrinfo for listing the extents of the reprojected shapefile. Just to add that the extents of the reprojected shapefile might well be what the OP wanted to know, but it is not what he asked for: the corner coordinates from a gtiff translated into another projection system Here a picture from osgeo.org/mapguide which shows the difference (although the picture was made to illustrate something different) http://trac.osgeo.org/mapguide/attachment/wiki/MapGuideRfc51/MapToLayerCS.png Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: gdalinfo: corner coords in other projection
On 24/05/2011 18:15, Lefman, Jonathan ERDC-TEC-VA wrote: Hi all, Is there an executable utility that takes the corner coordinates from a gtiff and translates them into another projection system? For example, I have a gtiff in UTM and I want to translate the coordinates of the bounding box (or corners) into WGS84. Or is there a way to do this non-programmatically using a combination of executable utilities? Couldn't you just simply use gdalinfo's output, see below? Hermann $ gdalinfo highway_number_3035.tif Driver: GTiff/GeoTIFF Files: highway_number_3035.tif Size is 10384, 7031 Coordinate System is: PROJCS[ETRS89 / ETRS-LAEA, GEOGCS[ETRS89, DATUM[European_Terrestrial_Reference_System_1989, SPHEROID[GRS 1980,6378137,298.2572221010002, AUTHORITY[EPSG,7019]], AUTHORITY[EPSG,6258]], PRIMEM[Greenwich,0], UNIT[degree,0.0174532925199433], AUTHORITY[EPSG,4258]], PROJECTION[Lambert_Azimuthal_Equal_Area], PARAMETER[latitude_of_center,52], PARAMETER[longitude_of_center,10], PARAMETER[false_easting,4321000], PARAMETER[false_northing,321], UNIT[metre,1, AUTHORITY[EPSG,9001]], AUTHORITY[EPSG,3035]] Origin = (270937.427360458299518,5929526.137714201584458) Pixel Size = (701.323508343892058,-701.323508343892058) Metadata: AREA_OR_POINT=Area Image Structure Metadata: INTERLEAVE=BAND Corner Coordinates: Upper Left ( 270937.427, 5929526.138) ( 67d52'39.91W, 53d13'23.14N) Lower Left ( 270937.427, 998520.551) ( 29d49'49.65W, 22d27'34.90N) Upper Right ( 7553480.738, 5929526.138) ( 81d 0' 6.31E, 59d42'39.88N) Lower Right ( 7553480.738, 998520.551) ( 42d23'31.15E, 25d49'31.87N) Center ( 3912209.083, 3464023.344) ( 3d44'23.41E, 54d 7'25.38N) ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: FileGDB OGR driver test
On 04/04/2011 05:05, Ragi Burhum wrote: Hello list, I am trying to test a new version of the FileGDB driver for OGR, but I lack enough FileGDBs to test :) Is the lack of FileGDBs still an issue? At work, we have several hundreds of them and I could check if I can make some available. If you have an *ArcGIS 10* FileGDB and would like to test the FileGDB driver (and report back), feel free to grab a gdal-trunk binary I made for Windows. A related question: will there ever be a FileGDB driver for a non-Windows OS? Is this feasible at all? Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: gdal_rasterize 1.8.0 options
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 ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: gdal_rasterize 1.8.0 options
On 05/02/2011 13:12, Even Rouault wrote: 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. Hmm. This page confirms what I mentioned in ticket #2869: Double_Stereographic is simply the ESRI alias of Oblique_Stereographic (EPSG code 9809). About r21627: I am afraid I do not understand enough about the inner workings of GDAL/OGR to make a judgement about whether this change was appropriate or not. But by the end of the day and from what I remember when looking somewhat deeper into the issue, about 2 years ago: there is one stereographic projection (known as the USGS formula or the approach given by Snyder) which translates into +proj=stere, and there is the other Double_Stereographic, aka Oblique_Stereographic which is expected to end up as +proj=sterea. Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] gdal_translate: GDAL: Should not happen. Cannot find laea_grid.asc
Hi, This is most likely a minor issue, but it looks to me that the above should not happen message always appears when gdal_translating a .vrt file that refers to an .asc file. There are 2 GDALClose() for the .asc. Hermann $ gdal_translate --debug on laea_grid.vrt laea_grid.tif \ -co compress=lzw -co tiled=yes GDAL: GDALOpen(laea_grid.vrt, this=0x236e6d0) succeeds as VRT. Input file size is 1, 1 0GDAL: GDALDefaultOverviews::OverviewScan() GDAL: GDALOpen(GTIFF_RAW:laea_grid.tif, this=0x2372180) succeeds as GTiff. GDAL: GDALDatasetCopyWholeRaster(): 9728*256 swaths, bInterleave=0 AAIGrid: Loaded SRS from laea_grid.prj GDAL: GDALOpen(laea_grid.asc, this=0x24324b0) succeeds as AAIGrid. ...10...20...30...40...50...60...70...80...90...100 - done. GDAL: GDALClose(laea_grid.tif, this=0x2372180) GDAL: GDALClose(laea_grid.vrt, this=0x236e6d0) GDAL: GDALClose(laea_grid.asc, this=0x24324b0) GDAL: GDALClose(laea_grid.asc, this=0x2371ad0) GDAL: Should not happen. Cannot find laea_grid.asc, this=0x2371ad0 in phSharedDatasetSet GDAL: GDALDeregister_GTiff() called. ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] ERROR 1: LZWDecode:Wrong length of decoded string: data probably corrupted
Hi,
Over the weekend, I tried to gdal_rasterize a 350+ MB shapefile with
149671 polygon features into a 64116 x 61850 GeoTIFF.
I first had some trouble while trying to rasterize everything at once,
because after some time, the gdal_rasterize process got killed by the
kerned (that's how I understood it).
I then tried to loop through the shapefile's CODE values, something like
the bash code [1] below. I ended up with different errors, see [2]
Any idea what is going on and how to avoid these errors? (I am using
GDAL v1.8dev from trunk)
Thanks in advance, Hermann
PS I tried without LZW compression, but this generates a 3+GB GeoTIFF,
which in return seems to lead to other (performance) problems.
[1]
for file in *.shp ; do
for code in {0..20} ; do
if [[ $code = 0 ]] ; then
gdal_rasterize --debug on --config GDAL_CACHEMAX 1000 \
-ot byte -co compress=lzw -tr 2 2 -tap -l ${file%.shp} \
-burn $code -where CODE=\'$code\' \
$file ${file%.shp}_2m.tif ${file%.shp}_2m.log 21
else
gdal_rasterize --debug on --config GDAL_CACHEMAX 1000 \
-l ${file%.shp} -burn $code -where CODE=\'$code\' \
$file ${file%.shp}_2m.tif ${file%.shp}_2m.log 21
fi
done
done
[2]
$ head -90 uk001l_london_2m.log
OGR: OGROpen(shpfiles/uk001l_london.shp/0x764610) succeeded as ESRI
Shapefile.
GDAL:
GDALDriver::Create(GTiff,uk001l_london_2m.tif,64116,61850,1,Byte,0x764ee0)
Shape: 149671 features read on layer 'uk001l_london'.
GDAL: GDALClose(uk001l_london_2m.tif, this=0x782e00)
GDAL: GDALDeregister_GTiff() called.
OGR: OGROpen(shpfiles/uk001l_london.shp/0x2080610) succeeded as ESRI
Shapefile.
GDAL: GDALOpen(uk001l_london_2m.tif, this=0x209feb0) succeeds as GTiff.
0...10...20...30...40...50...60.OGR:
OGROpen(shpfiles/uk001l_london.shp/0x10c1610) succeeded as ESRI Shapefile.
GDAL: GDALOpen(uk001l_london_2m.tif, this=0x10e0eb0) succeeds as GTiff.
.0...70..10.2080.30.90..Shape: 149671 features read on layer
'uk001l_london'.
40.GDAL: GDALClose(uk001l_london_2m.tif, this=0x209feb0)
GDAL: GDALDeregister_GTiff() called.
100 - done.
..50..OGR: OGROpen(shpfiles/uk001l_london.shp/0x225e610) succeeded as
ESRI Shapefile.
GDAL: GDALOpen(uk001l_london_2m.tif, this=0x227deb0) succeeds as GTiff.
0601020.3070.OGR:
OGROpen(shpfiles/uk001l_london.shp/0x1543610) succeeded as ESRI Shapefile.
GDAL: GDALOpen(uk001l_london_2m.tif, this=0x1562eb0) succeeds as GTiff.
0...10...2040.3080.50..4090...ERROR 1: LZWDecode:Wrong
length of decoded string: data probably corrupted at scanline 36908
ERROR 1: TIFFReadEncodedStrip() failed.
ERROR 1: IReadBlock failed at X offset 0, Y offset 36908
ERROR 1: GetBlockRef failed at X block offset 0, Y block offset 36908
Shape: 149671 features read on layer 'uk001l_london'.
..GDAL: GDALClose(uk001l_london_2m.tif, this=0x227deb0)
GDAL: GDALDeregister_GTiff() called.
50...Shape: 149671 features read on layer 'uk001l_london'.
.ERROR 1: LZWDecode:Wrong length of decoded string: data probably
corrupted at scanline 36908
ERROR 1: TIFFReadEncodedStrip() failed.
ERROR 1: IReadBlock failed at X offset 0, Y offset 36908
ERROR 1: GetBlockRef failed at X block offset 0, Y block offset 36908
Shape: 149671 features read on layer 'uk001l_london'.
GDAL: GDALClose(uk001l_london_2m.tif, this=0x10e0eb0)
GDAL: GDALDeregister_GTiff() called.
100 - done.
GDAL: GDALClose(uk001l_london_2m.tif, this=0x1562eb0)
GDAL: GDALDeregister_GTiff() called.
OGR: OGROpen(shpfiles/uk001l_london.shp/0xa21610) succeeded as ESRI
Shapefile.
GDAL: GDALOpen(uk001l_london_2m.tif, this=0xa40eb0) succeeds as GTiff.
0Warning 1: LZWDecode:LZWDecode: Strip 863 not terminated with EOI code
ERROR 1: LZWDecode:Not enough data at scanline 863 (short 17574 bytes)
ERROR 1: TIFFReadEncodedStrip() failed.
ERROR 1: IReadBlock failed at X offset 0, Y offset 863
ERROR 1: GetBlockRef failed at X block offset 0, Y block offset 863
Shape: 149671 features read on layer 'uk001l_london'.
GDAL: GDALClose(uk001l_london_2m.tif, this=0xa40eb0)
GDAL: GDALDeregister_GTiff() called.
OGR: OGROpen(shpfiles/uk001l_london.shp/0x1cd3610) succeeded as ESRI
Shapefile.
GDAL: GDALOpen(uk001l_london_2m.tif, this=0x1cf2eb0) succeeds as GTiff.
0Warning 1: LZWDecode:LZWDecode: Strip 863 not terminated with EOI code
ERROR 1: LZWDecode:Not enough data at scanline 863 (short 17574 bytes)
ERROR 1: TIFFReadEncodedStrip() failed.
ERROR 1: IReadBlock failed at X offset 0, Y offset 863
ERROR 1: GetBlockRef failed at X block offset 0, Y block offset 863
Shape: 149671 features read on layer 'uk001l_london'.
GDAL: GDALClose(uk001l_london_2m.tif, this=0x1cf2eb0)
GDAL: GDALDeregister_GTiff() called.
OGR: OGROpen(shpfiles/uk001l_london.shp/0x1cc5610) succeeded as ESRI
Shapefile.
GDAL: GDALOpen(uk001l_london_2m.tif, this=0x1ce4eb0) succeeds as GTiff.
0Warning 1: LZWDecode:LZWDecode: Strip 863 not terminated with EOI code
ERROR 1: LZWDecode:Not enough data at scanline 863 (short 17574
Re: [gdal-dev] ERROR 1: LZWDecode:Wrong length of decoded string: data probably corrupted
Even,
Thanks for the quick feedback.
Just to clarify: last week, I compiled revision 21190 from
http://svn.osgeo.org/gdal/trunk/gdal, configured
--with-libtiff=internal. I have 8G of memory, which is about 20 times
the size of the shapefile. However, I might have been processing other
shapefiles in parallel when running into the out-of-memory error. I will
double-check.
About using an uncompressed TIFF for intermediate work:
Doesn't this lead to performance issues, as the uncompressed TIFF (say:
3G) has be written to disc so many times (?)
Thanks again, Hermann
On 06/12/2010 11:19, Even Rouault wrote:
Selon Hermann Peiferpei...@gmx.eu:
Herman,
About your attempt to try to rasterize everything at once, it is likely that you
ran out of virtual memory as gdal_rasterize will load all the geometries in
memory before proceeding to the rasterization (it requires probably much more
memory than the on-disk size of the shapefile). If you can debug, you could set
up a breakpoint before the call to GDALRasterizeGeometries() at line 293 of
gdal_rasterize.cpp. If my hypothesis is right, then the process will be killed
before reaching the breakpoint. (There's a GDALRasterizeLayers() API that would
avoid loading all the OGR layer at once, but it is unused by gdal_rasterize
application, apparently because there are some options that couldn't be done
with it)
As far as the error about LZW is concerned, do you use internal libtiff or your
system libtiff ? If using external libtiff, you could try using internal libtiff
instead. And if you are already using internal libtiff, then it's likely an
evidence there's a remaining bug in how it deals with updating... But anyway I
would recommand you against using a compression scheme when doing multiple
updates of the file. Due to the way how libtiff manages strips/tiles allocation,
it will very likely to suboptimal file size because strips/tiles will be
rewritten many times. I'd rather use uncompressed TIFF for intermediate work and
as a final step do a gdal_translate -of LZW.
Best regards,
Even
Hi,
Over the weekend, I tried to gdal_rasterize a 350+ MB shapefile with
149671 polygon features into a 64116 x 61850 GeoTIFF.
I first had some trouble while trying to rasterize everything at once,
because after some time, the gdal_rasterize process got killed by the
kerned (that's how I understood it).
I then tried to loop through the shapefile's CODE values, something like
the bash code [1] below. I ended up with different errors, see [2]
Any idea what is going on and how to avoid these errors? (I am using
GDAL v1.8dev from trunk)
Thanks in advance, Hermann
PS I tried without LZW compression, but this generates a 3+GB GeoTIFF,
which in return seems to lead to other (performance) problems.
[1]
for file in *.shp ; do
for code in {0..20} ; do
if [[ $code = 0 ]] ; then
gdal_rasterize --debug on --config GDAL_CACHEMAX 1000 \
-ot byte -co compress=lzw -tr 2 2 -tap -l ${file%.shp} \
-burn $code -where CODE=\'$code\' \
$file ${file%.shp}_2m.tif ${file%.shp}_2m.log 21
else
gdal_rasterize --debug on --config GDAL_CACHEMAX 1000 \
-l ${file%.shp} -burn $code -where CODE=\'$code\' \
$file ${file%.shp}_2m.tif ${file%.shp}_2m.log 21
fi
done
done
[2]
$ head -90 uk001l_london_2m.log
OGR: OGROpen(shpfiles/uk001l_london.shp/0x764610) succeeded as ESRI
Shapefile.
GDAL:
GDALDriver::Create(GTiff,uk001l_london_2m.tif,64116,61850,1,Byte,0x764ee0)
Shape: 149671 features read on layer 'uk001l_london'.
GDAL: GDALClose(uk001l_london_2m.tif, this=0x782e00)
GDAL: GDALDeregister_GTiff() called.
OGR: OGROpen(shpfiles/uk001l_london.shp/0x2080610) succeeded as ESRI
Shapefile.
GDAL: GDALOpen(uk001l_london_2m.tif, this=0x209feb0) succeeds as GTiff.
0...10...20...30...40...50...60.OGR:
OGROpen(shpfiles/uk001l_london.shp/0x10c1610) succeeded as ESRI Shapefile.
GDAL: GDALOpen(uk001l_london_2m.tif, this=0x10e0eb0) succeeds as GTiff.
.0...70..10.2080.30.90..Shape: 149671 features read on layer
'uk001l_london'.
40.GDAL: GDALClose(uk001l_london_2m.tif, this=0x209feb0)
GDAL: GDALDeregister_GTiff() called.
100 - done.
..50..OGR: OGROpen(shpfiles/uk001l_london.shp/0x225e610) succeeded as
ESRI Shapefile.
GDAL: GDALOpen(uk001l_london_2m.tif, this=0x227deb0) succeeds as GTiff.
0601020.3070.OGR:
OGROpen(shpfiles/uk001l_london.shp/0x1543610) succeeded as ESRI Shapefile.
GDAL: GDALOpen(uk001l_london_2m.tif, this=0x1562eb0) succeeds as GTiff.
0...10...2040.3080.50..4090...ERROR 1: LZWDecode:Wrong
length of decoded string: data probably corrupted at scanline 36908
ERROR 1: TIFFReadEncodedStrip() failed.
ERROR 1: IReadBlock failed at X offset 0, Y offset 36908
ERROR 1: GetBlockRef failed at X block offset 0, Y block offset 36908
Shape: 149671 features read on layer 'uk001l_london'.
..GDAL: GDALClose(uk001l_london_2m.tif, this=0x227deb0)
GDAL: GDALDeregister_GTiff() called.
Re: [gdal-dev] ERROR 1: LZWDecode:Wrong length of decoded string: data probably corrupted
On 06/12/2010 13:15, Even Rouault wrote: Selon Hermann Peiferpei...@gmx.eu: Even, Thanks for the quick feedback. Just to clarify: last week, I compiled revision 21190 from http://svn.osgeo.org/gdal/trunk/gdal, configured --with-libtiff=internal. Hum OK so there's definitely a bug somewhere. Would be worth opening a ticket if you had a way of reproducing those errors. Ideally with a smaller shapefile. Would be interesting to see if it works better with another compression scheme, let's say DEFLATE. Thanks for the hint with DEFLATE. And yes, I am working on a 64 bit Debian/Linux system. As far as I can see: the LZW error only occurred while processing this one file, which also happens to be the biggest out of 231 shapefiles that I am handling in my current task. So it might be difficult to reproduce the error with a smaller file. 1 additional question, as I do not really understand the inner workings of gdal_rasterize: Does it help at all what I did when trying to reduce the memory demand: to loop through the geometries by some CODE value? Does gdal_rasterize always load *all* the geometries into memory, independent of the -where filter ? Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: warp, translate and world files
On 25/10/2010 11:17, Paul Meems wrote: Hi list, I've got several png files with worldfiles. I need them in another projection (28992 -- 3785). I'm using gdalwarp to reproject to a tiff file and then I use gdal_translate to create a png file again. My original png file is about 25 kB, the warped tiff is over 2MB and the translated png file is over 80 kB. I am not quite sure why you don't use -of png to save the creation of the intermediate tiff file. You might also be able to use -of byte which creates an 8 bit PNG (if this format is sufficient for your purpose). The size difference 20k vs 80k could be due to the fact that the source file is PNG8 and the generated file is PNG24. There may be other reasons, though. Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: problems to cut form large image
On 22/10/2010 18:14, Jan Tappenbeck wrote: C:\Program Files (x86)\FWTools2.4.7gdalwarp -r near -te 2608000 55900100 2613104.5424 5582270.5709 -of GTiff -o F:\wt_3a.tif D:\weissenturm\streifen_3_GK2.tif What do you expect -o to do? BTW: -r near and -of GTiff can be dropped, as both are default options. Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: ERROR 4: ... is a grib file, but no raster dataset was successfully identified.
Hello again. I am still running into this error with gdalinfo: $ gdalinfo HRES_ENS_2010101000+000.grib2 ERROR 4: HRES_ENS_2010101000+000.grib2 is a grib file, but no raster dataset was successfully identified. gdalinfo failed - unable to open 'HRES_ENS_2010101000+000.grib2'. However, other tools (wgrib2 and grib_ls from libgrib-api-tools) read the same file without any problems, see below. Any ideas? The file is here: http://transfer.eea.europa.eu/download/07ae94812cbe4823600b4d5ae64aadca Thanks in advance, Hermann $ wgrib2 HRES_ENS_2010101000+000.grib2 1:0:d=2010101000:MASSDEN:surface:anl:ENS=hi-res ctl chemical=Water Vapour 2:675793:d=2010101000:MASSDEN:500 m above ground:anl:ENS=hi-res ctl chemical=Water Vapour 3:1351586:d=2010101000:MASSDEN:1000 m above ground:anl:ENS=hi-res ctl chemical=Water Vapour 4:2027379:d=2010101000:MASSDEN:3000 m above ground:anl:ENS=hi-res ctl chemical=Water Vapour 5:2703172:d=2010101000:MASSDEN:surface:anl:ENS=hi-res ctl chemical=Nitrous Oxide 6:3378965:d=2010101000:MASSDEN:500 m above ground:anl:ENS=hi-res ctl chemical=Nitrous Oxide 7:4054758:d=2010101000:MASSDEN:1000 m above ground:anl:ENS=hi-res ctl chemical=Nitrous Oxide 8:4730551:d=2010101000:MASSDEN:3000 m above ground:anl:ENS=hi-res ctl chemical=Nitrous Oxide 9:5406344:d=2010101000:MASSDEN:surface:anl:ENS=hi-res ctl chemical=Ammonium 10:6082137:d=2010101000:MASSDEN:500 m above ground:anl:ENS=hi-res ctl chemical=Ammonium 11:6757930:d=2010101000:MASSDEN:1000 m above ground:anl:ENS=hi-res ctl chemical=Ammonium 12:7433723:d=2010101000:MASSDEN:3000 m above ground:anl:ENS=hi-res ctl chemical=Ammonium 13:8109516:d=2010101000:MASSDEN:surface:anl:ENS=hi-res ctl chemical=Nitrogen Dioxide 14:8785309:d=2010101000:MASSDEN:500 m above ground:anl:ENS=hi-res ctl chemical=Nitrogen Dioxide 15:9461102:d=2010101000:MASSDEN:1000 m above ground:anl:ENS=hi-res ctl chemical=Nitrogen Dioxide 16:10136895:d=2010101000:MASSDEN:3000 m above ground:anl:ENS=hi-res ctl chemical=Nitrogen Dioxide 17:10812688:d=2010101000:MASSDEN:surface:anl:ENS=hi-res ctl chemical=??? see code 4.230 18:11488481:d=2010101000:MASSDEN:500 m above ground:anl:ENS=hi-res ctl chemical=??? see code 4.230 19:12164274:d=2010101000:MASSDEN:1000 m above ground:anl:ENS=hi-res ctl chemical=??? see code 4.230 20:12840067:d=2010101000:MASSDEN:3000 m above ground:anl:ENS=hi-res ctl chemical=??? see code 4.230 ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: combine shapefiles?
On 13/10/2010 09:10, Jukka Rahkonen wrote: However, the .shp part of shapefile can not grow bigger than 4 gigabytes. Did you ever succeed in creating a, say: 3.9G .shp file? I do remember that some time ago, I was running into problems around 2G, which wasn't too much of a surprise, as [1] suggests that: There is a 2 GB size limit for any shapefile component file,... Hermann [1] http://webhelp.esri.com/arcgiSDEsktop/9.3/index.cfm?TopicName=Geoprocessing_considerations_for_shapefile_output ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: gdal_rasterize -tr and -te
On 05/10/2010 19:40, Even Rouault wrote: Le lundi 04 octobre 2010 17:23:47, Frank Warmerdam a écrit : gdal_translate, gdalwarp, ... Even, I would suggest -tap for target aligned pixels to go with -tr. I would think that gdalwarp, gdal_rasterize, and perhaps gdal_merge.py and gdalbuildvrt would be appropriate applications to support this switch. Foks, I've followed Frank's suggestion. Please test the patch attached to http://trac.osgeo.org/gdal/ticket/3772 (applies cleanly on latest trunk, r20770) and report if it works as you expect. Then I'll commit it. I don't see how it would be appropriate for gdal_translate. Best regards, I applied the patch and the corner coordinates seem to work fine for gdal_rasterize and gdalwarp. I haven't tried gdalbuildvrt. Perhaps someone else would like to test it? On the other hand... The new, target aligned tiffs do only have NODATA values, which seems to be related to ERROR 1: $ gdal_rasterize --debug on --config GDAL_CACHEMAX 2000 -ot byte -a_nodata 0 -co compress=lzw -tr 255 255 -tap -l lu001l_luxembourg -burn 1 -where CODE=111 lu001l_luxembourg.shp out.tif OGR: OGROpen(lu001l_luxembourg.shp/0x1978ae0) succeeded as ESRI Shapefile. GDAL: QuietDelete(out.tif) invoking Delete() GDAL: GDALOpen(out.tif, this=0x19cbe00) succeeds as GTiff. GDAL: GDALDefaultOverviews::OverviewScan() GDAL: GDALClose(out.tif, this=0x19cbe00) GDAL: GDALDriver::Create(GTiff,out.tif,224,321,1,Byte,0x1979ab0) ERROR 1: Type mismatch or improper type of arguments to = operator. GTiff: Adjusted bytes to write from 8064 to 7392. GDAL: GDALClose(out.tif, this=0x19caf40) GDAL: GDALDeregister_GTiff() called. I am not sure if this new ERROR is related to your patch or some other changes which came through checking out -r 20770. This error message is obvioulsy defined in from ogr/swq_op_general.cpp Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] ERROR 4: ... is a grib file, but no raster dataset was successfully identified.
It is here for the next 10 days: http://transfer.eea.europa.eu/download/6694e0deb608337a8d0e6c1fb0ad5068 Thanks for your time, Hermann On 04/10/2010 15:27, brian wrote: Hermann Do you have a link to the grib file you can share? Brian On Mon, 2010-10-04 at 15:16 +0200, Hermann Peifer wrote: Hi All, I am getting this error for a grib2 file: $ gdalinfo HRES_ENS_2010100300+000.grib2 ERROR 4: HRES_ENS_2010100300+000.grib2 is a grib file, but no raster dataset was successfully identified. gdalinfo failed - unable to open 'HRES_ENS_2010100300+000.grib2'. I am using gdal from trunk which has, among others, support for these formats: $ gdalinfo --formats | egrep JP|GRIB JPEG (rwv): JPEG JFIF JPEG2000 (rwv): JPEG-2000 part 1 (ISO/IEC 15444-1) GRIB (rov): GRIdded Binary (.grb) I found some older postings to the mailing list on a similar error, but whatever was mentioned there did not help in my case. Thanks for your time, Hermann ___ 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] ERROR 4: ... is a grib file, but no raster dataset was successfully identified.
Hi All, I am getting this error for a grib2 file: $ gdalinfo HRES_ENS_2010100300+000.grib2 ERROR 4: HRES_ENS_2010100300+000.grib2 is a grib file, but no raster dataset was successfully identified. gdalinfo failed - unable to open 'HRES_ENS_2010100300+000.grib2'. I am using gdal from trunk which has, among others, support for these formats: $ gdalinfo --formats | egrep JP|GRIB JPEG (rwv): JPEG JFIF JPEG2000 (rwv): JPEG-2000 part 1 (ISO/IEC 15444-1) GRIB (rov): GRIdded Binary (.grb) I found some older postings to the mailing list on a similar error, but whatever was mentioned there did not help in my case. Thanks for your time, Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: gdal_rasterize -tr and -te
On 04/10/2010 12:45, Vincent Schut wrote: On 10/01/2010 03:20 PM, Hermann Peifer wrote: On 01/10/2010 07:03, Jukka Rahkonen wrote: Hermann Peiferpeiferat gmx.eu writes: At work, we are taking this standard grid issue pretty serious, but indeed, we might be the only ones worldwide with such a business rule. You are not alone, we are reprojecting our rasters to standard grid because it helps in making seamless mosaics from the reprojected images. We are widening the target extents a little bit to all directions with a python script to suit the grid. So we are already 2 ;-) Actually, I tentatively thought there would be more. This was why I wrote to the mailing list in the first place. I am aware that this (non-)issue can be solved in a shell one-liner or perhaps X lines of some other script code. I just thought that if there is a general interest, it could perhaps be built into the library. Make that 3 :-) We also often do this, when different datasets (from several satellite sensors) need to be used together. Cause I do almost all my gdal stuff in python, I just coded this grid matching into the relevant python scripts. However, if support for it would be built into the gdal command line utilities, we would undoubtly use it often. Thanks for the time, Vincent. Even, I guess that 3 interested GDAL users (on a global scale) are not enough to make you reconsider your initial statement: This is a bit too specialized requirement... ? Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: ERROR 4: ... is a grib file, but no raster dataset was successfully identified.
Brian, Indeed, the file seems to be suspicious. See also my other mail where I used wgrib2, as suggested by Bill. I will have to contact the data provider. Thanks again for your time, Hermann On 04/10/2010 16:24, brian wrote: Hermann GDAL uses degrib to read grib files. I tryed the product in both degrib and gribdump, neither app will read it. Either the grib file is corrupt or it is a varient that neither app will read. Brian r...@octopi ~ $ degrib HRES_ENS_2010100300+000.grib2 -I ERROR: with inventory, so far: MsgNum, Byte, GRIB-Version, elem, level, reference(UTC), valid(UTC), Proj(hr) 1.0, 0, 2, (null), (null), 01/01/1970 00:00, 01/01/1970 00:00, 0.00 ERROR: In call to GRIB2Inventory. ERROR: Ran out of file in Section 1 ERROR: Problems with section 1 r...@octopi ~ $ degrib HRES_ENS_2010100300+000.grib2 -C -Csv ERROR: In call to Grib2Convert. ERROR: In call to ReadGrib2Record. Inside ReadGrib2Record.. Calling FindSectLen ERROR: Section 1 labeled as 21 r...@octopi ~ $ gribdump -v HRES_ENS_2010100300+000.grib2 Oct 04 14:21:53 gribdump[22297]: Starting Up Oct 04 14:21:53 gribdump[22297]: oversize GRIB product 2 Oct 04 14:21:53 gribdump[22297]: GRIB product: 2 Oct 04 14:21:53 gribdump[22297]: oversize GRIB product 2 Oct 04 14:21:53 gribdump[22297]: GRIB product: 4 Oct 04 14:21:53 gribdump[22297]: oversize GRIB product 2 Oct 04 14:21:53 gribdump[22297]: GRIB product: 6 Oct 04 14:21:53 gribdump[22297]: oversize GRIB product 2 Oct 04 14:21:53 gribdump[22297]: GRIB product: 8 Oct 04 14:21:53 gribdump[22297]: oversize GRIB product 2 Oct 04 14:21:53 gribdump[22297]: GRIB product: 10 Oct 04 14:21:53 gribdump[22297]: oversize GRIB product 2 Oct 04 14:21:53 gribdump[22297]: GRIB product: 12 Oct 04 14:21:53 gribdump[22297]: oversize GRIB product 2 Oct 04 14:21:53 gribdump[22297]: GRIB product: 14 Oct 04 14:21:53 gribdump[22297]: oversize GRIB product 2 Oct 04 14:21:53 gribdump[22297]: GRIB product: 16 Oct 04 14:21:53 gribdump[22297]: oversize GRIB product 2 Oct 04 14:21:53 gribdump[22297]: GRIB product: 18 Oct 04 14:21:53 gribdump[22297]: oversize GRIB product 2 Oct 04 14:21:53 gribdump[22297]: GRIB product: 20 Oct 04 14:21:53 gribdump[22297]: oversize GRIB product 2 Oct 04 14:21:53 gribdump[22297]: GRIB product: 22 Oct 04 14:21:53 gribdump[22297]: oversize GRIB product 2 Oct 04 14:21:53 gribdump[22297]: GRIB product: 24 Oct 04 14:21:53 gribdump[22297]: oversize GRIB product 2 Oct 04 14:21:53 gribdump[22297]: GRIB product: 26 Oct 04 14:21:53 gribdump[22297]: oversize GRIB product 2 Oct 04 14:21:53 gribdump[22297]: GRIB product: 28 Oct 04 14:21:53 gribdump[22297]: oversize GRIB product 2 Oct 04 14:21:53 gribdump[22297]: GRIB product: 30 Oct 04 14:21:53 gribdump[22297]: oversize GRIB product 2 Oct 04 14:21:53 gribdump[22297]: GRIB product: 32 Oct 04 14:21:53 gribdump[22297]: oversize GRIB product 2 Oct 04 14:21:53 gribdump[22297]: GRIB product: 34 Oct 04 14:21:53 gribdump[22297]: oversize GRIB product 2 Oct 04 14:21:53 gribdump[22297]: GRIB product: 36 Oct 04 14:21:53 gribdump[22297]: oversize GRIB product 2 Oct 04 14:21:53 gribdump[22297]: GRIB product: 38 Oct 04 14:21:53 gribdump[22297]: oversize GRIB product 2 Oct 04 14:21:53 gribdump[22297]: GRIB product: 40 Oct 04 14:21:53 gribdump[22297]: EOF on input Oct 04 14:21:53 gribdump[22297]: Exiting Oct 04 14:21:53 gribdump[22297]: 0 WMO msgs, 0 GRIBs decoded, 0 written On Mon, 2010-10-04 at 16:06 +0200, Hermann Peifer wrote: It is here for the next 10 days: http://transfer.eea.europa.eu/download/6694e0deb608337a8d0e6c1fb0ad5068 Thanks for your time, Hermann On 04/10/2010 15:27, brian wrote: Hermann Do you have a link to the grib file you can share? Brian On Mon, 2010-10-04 at 15:16 +0200, Hermann Peifer wrote: Hi All, I am getting this error for a grib2 file: $ gdalinfo HRES_ENS_2010100300+000.grib2 ERROR 4: HRES_ENS_2010100300+000.grib2 is a grib file, but no raster dataset was successfully identified. gdalinfo failed - unable to open 'HRES_ENS_2010100300+000.grib2'. I am using gdal from trunk which has, among others, support for these formats: $ gdalinfo --formats | egrep JP|GRIB JPEG (rwv): JPEG JFIF JPEG2000 (rwv): JPEG-2000 part 1 (ISO/IEC 15444-1) GRIB (rov): GRIdded Binary (.grb) I found some older postings to the mailing list on a similar error, but whatever was mentioned there did not help in my case. Thanks for your time, Hermann ___ 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] Re: gdal_rasterize -tr and -te
On 04/10/2010 16:31, Even Rouault wrote: This raises quite a few questions : * what is a standard grid ? Ours starts at 0,0, then goes in all directions with the given tr * what syntax would you imagine to pass to gdal utilities to define how you want the rounding ? Some creation options like standard=yes, standard_origin=0,0 * do you want to round/ceil/floor ? or so that the adjusted extent includes the exact extent ? or so that the adjusted extent is contained inside the exact extent ? The adjusted extend includes the exact extent. * which gdal utilities would be impacted ? All the ones that create new raster files: gdal_rasterize, gdal_translate, gdalwarp, ... * ... and is there a concensus among the interested users ;-) ? Hmm. Let's see. ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: gdal_rasterize -tr and -te
On 01/10/2010 07:03, Jukka Rahkonen wrote: Hermann Peiferpeiferat gmx.eu writes: At work, we are taking this standard grid issue pretty serious, but indeed, we might be the only ones worldwide with such a business rule. You are not alone, we are reprojecting our rasters to standard grid because it helps in making seamless mosaics from the reprojected images. We are widening the target extents a little bit to all directions with a python script to suit the grid. So we are already 2 ;-) Actually, I tentatively thought there would be more. This was why I wrote to the mailing list in the first place. I am aware that this (non-)issue can be solved in a shell one-liner or perhaps X lines of some other script code. I just thought that if there is a general interest, it could perhaps be built into the library. Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] gdal_rasterize -tr and -te
Hi All, GDAL = 1.8.0 comes with the helpful switches -te and -tr. If the target extent is not specified, the extent of the output file will be the extent of the vector layers. This leads usually to corner coordinates like this: Upper Left ( 4923406.374, 3088693.597) Lower Right ( 4994261.374, 3010723.597) However, what I always need are corner coordinates which follow the standard grid, which starts at 0,0 and then grows in all directions with the given pixel size. So in the above case and for e.g. a 100m raster I really want that the target extent grows to: Upper Left ( 4923500, 3088600) Lower Right ( 4994200, 3010800) There is probably no way of convincing gdal_rasterize to do the maths, based on the provided tr and te values? Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] gdal_rasterize -tr and -te
On 30/09/2010 18:50, Even Rouault wrote:
Lower Right ( 4994200, 3010800)
There is probably no way of convincing gdal_rasterize to do the maths,
based on the provided tr and te values?
No there isn't. This is a bit too specialized requirement... Your best help is
some Python GDAL scripting to create the output file before using
gdal_rasterize
OK. Thanks for replying. Currently I have a simple one-liner in my
batch-processing shell scripts, something like:
myExtent=$(ogrinfo -al -so my.shp | awk '/^Extent/{ do the maths }')
At work, we are taking this standard grid issue pretty serious, but
indeed, we might be the only ones worldwide with such a business rule.
Hermann
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: cut a .dem file
On 08/09/2010 15:26, Carmelo Terrasi wrote: Hello everybody, does anybody know how to cut a smaller dem file from another one??? I have my original file /terrain.dem /but my /image.bmp/ is smaller so I would cut a portion just to create another suitable dem file for my /image.bmp/ Is it possible to do it or am I off-board? Look at the -srcwin and -projwin options at http://gdal.org/gdal_translate.html Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: Integer error in creating bigtiff image
On 24/08/2010 21:33, Even Rouault wrote: Try reducing under 2000 and it should work better. A too big value for the -wm parameter can even sometimes affect speed negatively, and that'd probably the case for your example where you are merging a lot of individual files into a big one. Just to add that this counter-intuitive behaviour is mentioned at: http://trac.osgeo.org/gdal/wiki/UserDocs/GdalWarp#WillincreasingRAMincreasethespeedofgdalwarp You could use: --debug on, and then observe how Src and Dst windows change depending on the chosen -wm value. (There is little point in having the Src window bigger than the input raster.) Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: projection unsupported?
It looks to me that Kertau RSO Malaya Meters is also known as EPSG:3168. So for the GDAL tools, you could use the EPSG code or this list of parameters and values: http://spatialreference.org/ref/epsg/3168/proj4/ Hermann On 25/08/2010 07:35, Christoph Dohmen wrote: Dear list, I got some data to deal with which use a projection called Kertau RSO Malaya Meters. From some earlier posts to proj.4 list I got the information that this is a kind of oblique mercator projection. But I can't find anything around gdal about how to deal with this projection. Do I have to change the SRS before using gdal to deal with this data? Or will there be any change inside gdal to support this projection? Are there any hints outside I didn' get so far? Best regards, Christoph ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: Merging shapefile
On 06/08/2010 08:39, Giacomo Piva wrote: Hi all, I need to merge two shapefile, is it possible with GDAL ? Thank you See this FAQ: http://trac.osgeo.org/gdal/wiki/FAQVector#HowcanImergehundredsofShapefiles Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: Cropping raster by vector files
On 08/07/2010 23:06, michal_drozdz wrote: Hi! I have a beginers problem. How to crop collar from scanned map with GDAL? I have set of georeferenced topo maps and polygon layer with boundaries of each sheet. I used gdalwarp with -cutline and -cl parameters and it worked but as a result I got map with 0 value at the place of collars. So when I upload maps into QGIS I still have a problem with overlaping collars. How to remove it forever:) Thanks for replys! Another option would be to use the coordinates of the polygons as follows: gdal_translate -projwin ulx uly lrx lry This would remove the map collars forever. See http://www.gdal.org/gdal_translate.html Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Re: Writing non-ASCII characters to shapefile
The .cpg files we generate through ArcGIS desktop contain the string: UTF-8 Some time ago, there was a mail on this list about problems in case of conflicting information in the .cpg file compared to the Language Driver ID (LDID) in the header of a dBASE file, see: http://lists.osgeo.org/pipermail/gdal-dev/2010-May/024619.html Hermann On 06/07/2010 04:13, Francis Markham wrote: Okay, I will take that approach then. Thank you all for your help. What specific value should I write into the .cpg? The string '65001' or the string 'utf-8' or something else? -Francis On 5 July 2010 21:53, Peter Hopfgartnerpeter.hopfgart...@r3-gis.com wrote: Hi Francis, what does not portable mean? ArcMap handles UTF-8 fine, if the correct encoding is written into the .cpg file. Recent shapelib should handle this fine, too. If there is any problem with a specific GIS program, a bug report for that GIS program might be the right thing to do. Regards, Peter On Mon, 2010-07-05 at 19:18 +1000, Francis Markham wrote: At the bottom of this page, for one: http://resources.arcgis.com/content/kbase?fa=articleShowd=21106 But honestly I've found hard to find information about this. I'd be very happy to be corrected if this is not the case! Cheers, Francis On 5 July 2010 19:09, Hermann Peiferpei...@gmx.eu wrote: Francis, you wrote: I have heard that the use of UTF-8 in shapefiles is not portable. Where did you hear this? Regards, Hermann On 03/07/2010 04:40, Francis Markham wrote: Hi there, I'm trying to write data from a Microsoft Excel .xls file into a shapefile, using OGR's Python bindings in Python 2.6. This is going well, but I am having some problems when I try to write values that contain so-called smart quotes. Smart quotes are special characters, defined as characters 0x91 through 0x94 in Windows-1252 ( see http://msdn.microsoft.com/en-au/goglobal/cc305145.aspx ). What is the best way to save this data to a shapefile using OGR? I need the shapefile to be interoperable with other programs, including but not limited to ESRI products. While I assume I could simply translate these characters to standard ASCII, I would prefer not to if possible. I also haven't tested the shapefiles with data from other character encodings. I have heard that the use of UTF-8 in shapefiles is not portable. I am also aware that shapefile.cpg can store a shapefile's codepage. I don't know how to put these pieces together to create a portable solution, however. Apologies if this is a newbie question, but I can't find answers on the web. Thanks, Francis Markham ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev -- Dott. Peter Hopfgartner R3 GIS Srl - GmbH Via Johann Kravogl-Str. 2 I-39012 Meran/Merano (BZ) Email: peter.hopfgart...@r3-gis.com Tel. : +39 0473 494949 Fax : +39 0473 069902 www : http://www.r3-gis.com XING : http://www.xing.com/go/invita/8917535 ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: Writing non-ASCII characters to shapefile
In the mentioned mail from earlier this year, William Kyngesburye wrote: Maybe GDAL needs a creation option in the shapefile driver to set the LDID or to instead add a cpg file with an encoding value. As there was no reply to his mail, I assume that there is no GDAL/OGR creation option for the LDID and it is unlikely that there will be one in the near future. Hermann On 06/07/2010 13:04, Francis Markham wrote: It appears from http://lists.osgeo.org/pipermail/gdal-dev/2010-May/024619.html and http://resources.arcgis.com/content/kbase?fa=articleShowd=26015 that ESRI set the LDID to zero to indicate an unknown LDID. From my reading of the OGR sourcecode, it seems that OGR uses the shapelib default LDID of 0x57 (previously it was set to 0x3, or windows-1252). Is there any way to specify this value using the OGR API, or do I need to use shapelib directly to create my shapefiles in order to do this? Cheers, Francis On 6 July 2010 18:31, Hermann Peiferpei...@gmx.eu wrote: The .cpg files we generate through ArcGIS desktop contain the string: UTF-8 Some time ago, there was a mail on this list about problems in case of conflicting information in the .cpg file compared to the Language Driver ID (LDID) in the header of a dBASE file, see: http://lists.osgeo.org/pipermail/gdal-dev/2010-May/024619.html Hermann On 06/07/2010 04:13, Francis Markham wrote: Okay, I will take that approach then. Thank you all for your help. What specific value should I write into the .cpg? The string '65001' or the string 'utf-8' or something else? -Francis On 5 July 2010 21:53, Peter Hopfgartnerpeter.hopfgart...@r3-gis.com wrote: Hi Francis, what does not portable mean? ArcMap handles UTF-8 fine, if the correct encoding is written into the .cpg file. Recent shapelib should handle this fine, too. If there is any problem with a specific GIS program, a bug report for that GIS program might be the right thing to do. Regards, Peter On Mon, 2010-07-05 at 19:18 +1000, Francis Markham wrote: At the bottom of this page, for one: http://resources.arcgis.com/content/kbase?fa=articleShowd=21106 But honestly I've found hard to find information about this. I'd be very happy to be corrected if this is not the case! Cheers, Francis On 5 July 2010 19:09, Hermann Peiferpei...@gmx.eu wrote: Francis, you wrote: I have heard that the use of UTF-8 in shapefiles is not portable. Where did you hear this? Regards, Hermann On 03/07/2010 04:40, Francis Markham wrote: Hi there, I'm trying to write data from a Microsoft Excel .xls file into a shapefile, using OGR's Python bindings in Python 2.6. This is going well, but I am having some problems when I try to write values that contain so-called smart quotes. Smart quotes are special characters, defined as characters 0x91 through 0x94 in Windows-1252 ( see http://msdn.microsoft.com/en-au/goglobal/cc305145.aspx ). What is the best way to save this data to a shapefile using OGR? I need the shapefile to be interoperable with other programs, including but not limited to ESRI products. While I assume I could simply translate these characters to standard ASCII, I would prefer not to if possible. I also haven't tested the shapefiles with data from other character encodings. I have heard that the use of UTF-8 in shapefiles is not portable. I am also aware that shapefile.cpg can store a shapefile's codepage. I don't know how to put these pieces together to create a portable solution, however. Apologies if this is a newbie question, but I can't find answers on the web. Thanks, Francis Markham ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev -- Dott. Peter Hopfgartner R3 GIS Srl - GmbH Via Johann Kravogl-Str. 2 I-39012 Meran/Merano (BZ) Email: peter.hopfgart...@r3-gis.com Tel. : +39 0473 494949 Fax : +39 0473 069902 www : http://www.r3-gis.com XING : http://www.xing.com/go/invita/8917535 ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: Error while using gdalwarp- Values become Zero
Could it be that the top right corner of final.tif is simply outside the extent of translate.tif and filled with 0, during warping? Hermann On 30/06/2010 15:32, Monica Buescu wrote: Greetings I have TIF file that I produced using gdal_translate. I used the following expression: gdal_translate -of GTiff -a_ullr 350.25 42.75 354.75 36.75 GRIB_tests/pt.grib GRIB_tests/translate.tif All values are ok and it seems to be working. Then, I decided to gdalwarp it with this command: gdalwarp -t_srs EPSG:3763 translate.tif final.tif. This final.tif is all ok except the top right corner where all values became 0. I have tried this methodology with other final EPSG (WGS-UTM for instance) and I didn't get this problem. Thanks fr your help Best regards, Monica ___ 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] Re: gdalwarp large imagery faster?
On 24/06/2010 16:40, GeoSpatial - Kevin wrote: Hi all, I am trying to re-project a very large imagery (size 40,000 by 40,000 pixels, resolution ~120m) located at high latitudes (60N) to a standard wgs84. By searching some previous posts, I come up with the following command: gdalwarp -of EHdr -t_srs EPSG:4326 -tr 0.00 0.00 --config GDAL_CACHEMAX 200 -wm 200 -wo SKIP_NOSOURCE=YES -r bilinear -multi input_ice.tif output_ice.tif it is working now but the speed is very slow - likely this will take 12 hours... Any way to speed up or am I wrong selecting those parameters? Thanks. - multi doesn't work for me: CPLCreateThread: Fails to dummy implementation ERROR 1: CPLCreateThread() failed in ChunkAndWarpMulti() However, without the -multi switch, my input data (100m raster, 48500x41500, Type=Int16) is warped in 34 minutes or so, see here: time gdalwarp --debug on -of EHdr -t_srs EPSG:4326 -tr 0.00 0.00 --config GDAL_CACHEMAX 200 -wm 200 -wo SKIP_NOSOURCE=YES -r bilinear map2c_100.tif output.tif real34m13.134s user5m24.524s sys 0m9.973s If I raise the values for GDAL_CACHEMAX and -wm, it takes about 15 minutes. I am using GDAL 1.8dev, on a 64-bit Linux machine. Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: how to new image file or how to send data into the dataset
On 31/05/2010 10:40, mail2vajram wrote: i need to create a image file. Here i have brightness values of each pixel in a table. by using these values, how can i create an image file. You could add a few header lines to your table, so that it looks like ASCII grid format, see http://en.wikipedia.org/wiki/ESRI_grid After that, you could use gdal_translate to convert the values into an image format, see http://gdal.org/gdal_translate.html Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: OT: Need help with proj definition
On 13/05/2010 22:32, Stephen Woodbridge wrote: Hi all, I have a shapefile with the following in the prj file: GEOGCS[GCS_ITRF_2000,DATUM[D_ITRF_2000,SPHEROID[GRS_1980,6378137.0,298.257222101]],PRIMEM[Greenwich,0.0],UNIT[Degree,0.0174532925199433]] Anyone know how to convert this to a proj4 definition? Thanks, -Steve One option could be to simply convert the shapefile into GMT format. # -fid 0 keeps the GMT file small ogr2ogr -f GMT yourfile.gmt yourfile.shp -fid 0 You'll find something like this in the GMT file: # @Jp+proj=longlat +ellps=GRS80 +no_defs Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] gdal_translate -scale and negative values
Hi, I am trying to scale the values in an Int16 raster band. Basically, I want to have them divided by 100. It looks to me that negative values are not scaled at all, whatever ranges I give to -scale. What am I doing wrong? Hermann $ gdalinfo -mm infile.tif ... Band 1 Block=7200x1 Type=Int16, ColorInterp=Gray Computed Min/Max=-464.000,3122.000 NoData Value=-32768 ... $ gdal_translate infile.tif outfile.tif -scale -500 0 -5 0 Input file size is 7200, 3600 0...10...20...30...40...50...60...70...80...90...100 - done. $ gdalinfo -mm outfile.tif ... Band 1 Block=7200x1 Type=Int16, ColorInterp=Gray Computed Min/Max=-328.000,31.000 NoData Value=-32768 ... $ gdalinfo --version GDAL 1.7.1, released 2010/02/08 ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] gdal_translate -scale and negative values
didn't you notice that... Not really. I must have been blind. Thanks for your hint. This works fine. Hermann Original Message Subject: [gdal-dev] gdal_translate -scale and negative values From: Even Rouault even.roua...@mines-paris.org To: gdal-dev@lists.osgeo.org Cc: Hermann Peifer pei...@gmx.eu, Frank Warmerdam warmer...@pobox.com Date: 15/02/2010 20:53 Herman, didn't you notice that -32768 / 100 = -327.68 ~= -328 ... ;-) ? The issue is that when doing rescaling, gdal_translate doesn't currently explicetly set the source nodata value as the nodata value of the VRT source, so the pixels whose value is the nodata value are just rescaled as other pixels. However, you can get your expected result by explicetly setting the nodata value of the source : 1) gdal_translate -of VRT infile.tif outfile.vrt -scale -500 0 -5 0 2) edit outfile.vrt with your favorite text editor and add NODATA-32768/NODATA between ComplexSource and /ComplexSource 3) gdal_translate outfile.vrt outfile.tif I'm not sure why gdal_translate doesn't set the source nodata value by default. Looking at the history, it looks like it never did. Maybe Frank will have an opinion on if it would be approriate to automatically set the source nodata value as the nodata value of the ComplexSource ? (this behaviour could be controlled by an option flag if necessary) Best regards, Even Hi, I am trying to scale the values in an Int16 raster band. Basically, I want to have them divided by 100. It looks to me that negative values are not scaled at all, whatever ranges I give to -scale. What am I doing wrong? Hermann $ gdalinfo -mm infile.tif ... Band 1 Block=7200x1 Type=Int16, ColorInterp=Gray Computed Min/Max=-464.000,3122.000 NoData Value=-32768 ... $ gdal_translate infile.tif outfile.tif -scale -500 0 -5 0 Input file size is 7200, 3600 0...10...20...30...40...50...60...70...80...90...100 - done. $ gdalinfo -mm outfile.tif ... Band 1 Block=7200x1 Type=Int16, ColorInterp=Gray Computed Min/Max=-328.000,31.000 NoData Value=-32768 ... $ gdalinfo --version GDAL 1.7.1, released 2010/02/08 ___ 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] Re: no date value by gdal_merge from DEM to GTiff ?
Bruce Liang wrote: Hello, guys, I have several DEMs from USGS GTOPO30, and by merging them with gdal_merge.py, it seems the no data value can not be assigned? the command is like: gdal_merge.py -of GTiff -o Europe.tiff -n - -co PROFILE=GeoTIFF -co INTERLEAVE=PIXEL -co COMPRESS=NONE -co TILED=YES *.DEM and to check it with gdalinfo -stats Europe.tiff , i got: Band 1 Block=256x256 Type=UInt16, ColorInterp=Gray Minimum=1.000, Maximum=65535.000, Mean=30979.953, StdDev=27318.154 Metadata: STATISTICS_MINIMUM=1 STATISTICS_MAXIMUM=65535 STATISTICS_MEAN=30979.952856589 STATISTICS_STDDEV=27318.154016739 seems the no data value is now 65535 ? how to assign it correctly? No. It seems that your highest *data* value is now 65535 and there isn't any *no data* value defined. According to your GeoTIFF statistics, the average elevation in Europe is around 31000 meters above sea level. Which seems to be a bit too high, as far as I can judge. You might have also noted, that your GeoTIFF type is UInt16, which doesn't allow for negative values. Elevation data can however be negative: a good part of the Netherlands is below sea level, not to talk about the Dead Sea (-422m or so). Here is what I did the other day in order to merge GTOPO30 data: gdalwarp -srcnodata - -dstnodata - -co TILED=YES *.DEM outfile.tif In my case, the resulting GeoTIFF type was Int16 (same type as the *.DEMs), and the statistics are: Minimum=-405.000, Maximum=5417.000, Mean=359.426, StdDev=456.872 NoData Value=- Metadata: STATISTICS_MINIMUM=-405 STATISTICS_MAXIMUM=5417 STATISTICS_MEAN=359.42640750109 STATISTICS_STDDEV=456.87166274078 gdal_merge.py is an example Python script and I am not quite sure what its advantage is over gdalwarp. FrankW could perhaps shed some light on this issue. Hermann ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
