Re: [gdal-dev] gdalwarp resampling_method option

2017-10-05 Thread Mike Toews 
Hi,

I requested this feature 2 years ago:
https://trac.osgeo.org/gdal/ticket/6015

I probably had some workaround in Python, but I can't recall any details.

Cheers,
Mike

On 6 October 2017 at 00:59, sujankoirala2  wrote:
> Hi there,
> I was wondering if it would be possible to get the sum instead of the mean
> or (any other regridding method). In our data, in many cases we need to
> aggregate mass to coarser resolution. In theory, we can just multiply the
> output of average method with number of pixels. But, due to differing size
> of valid pixels (non nodata values) in each window, it won't produce an
> accurate result.
> I would be grateful if anyone can point me to 1) an alternate way to get the
> number of valid pixels in each regridding window) or 2) get the sum within
> the regridding window or 3) the location in the source where average method
> of resampling is implemented so that I can edit the source and build from
> scratch.
> Thank you.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Improve clipping of shapefiles

2017-09-11 Thread Mike Toews 
On 11 September 2017 at 21:09, Paul Meems  wrote:
> I have a large shapefile with over 2.8 million shapes (fishnet) and I have a
> border file with only 1 multipolygon.
>
> I'm trying to clip the fishnet with the border.
> Using code is takes about 5 min. using command line it takes even longer.
>
> This is my command:
> ogr2ogr fishnetClipped.shp fishnet.shp -clipsrc border.shp
> ...
> Can I somehow improve the speed?

If the fishnet shapefile does not have a spatial index, try creating
one. There are a few ways to do this:
http://www.gdal.org/drv_shapefile.html

I'm not sure if -clipsrc datasource will use the spatial index, but
it's worth trying.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] gdal on freebsd

2016-11-07 Thread Mike Toews 
On 7 November 2016 at 21:33, Jyrki Von Karkki
 wrote:
> Hi all
>
> I'm having a bit of an issue with the postgresql/postgis driver on freebsd.
> It simply does not appear installed.

How did you install GDAL? The PostgreSQL driver is not always
available since it not compiled by default.

If you are compiling GDAL from source, see --with-pg=ARG from the
configure script. Also be sure to have the development headers for
libpq.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Simple circle in OGR

2016-11-02 Thread Mike Toews 
On 3 November 2016 at 09:46, Renison  wrote:
> Hello all, I am trying to create a simple circle geometry in OGR. My approach
> so far has been something like this:
>
> SpatialReference defaultGeoSys = new SpatialReference("");
> defaultGeoSys.SetWellKnownGeogCS("EPSG:4326");
>
> //define a circle using point and buffer
> string pointlineString = "POINT (" + longitude + " " + latitude
> +")";
> var circle = OGR.Ogr.CreateGeometryFromWkt(ref pointlineString,
> defaultGeoSys);
> circle.Buffer(MilesToMeters(radiusInMiles), 2000);
>
> But the circle.Buffer(params,params) does not seem to create a circle. All I
> have is a point and a radius. Is there a way to achieve what I want?

Buffer operates in Cartesian space, and uses the same units as the
SRS, therefore these are degrees.

One strategy for a simple Point is to create a Cartesian buffer from
(0, 0) and project it to EPSG:2193 from a custom azimuthal equidistant
projection centred at the point location, e.g. a PROJ.4 string
"+proj=aeqd +lat_0=" + latitude + " +lon_0=" + longitude + " +x_0=0
+y_0=0"

See http://gis.stackexchange.com/a/121539/1872 for this recipe in R.
You can make something similar with OSR.

Cheers,
-Mike
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] VRT raster for reading formatted files?

2016-06-25 Thread Mike Toews 
Hi all,

I've been using the VRTRawRasterBand feature of the VRT format to read
unformatted (so-called "BINARY") 2D grid data, which has been very
successful to adapt to obscure formats.

However, I would like to know if there is a similar feature of the VRT
format to read formatted values (so-called "ASCII") from a
SourceFilename.

There are already a few GDAL drivers that do this already,
specifically AAIGrid / GRASSASCIIGrid and GSAG. Possibly others. But
it would be great to see if there is a generic approach with VRT.

An example of a SourceFilename is at:
https://github.com/mwtoews/MODFLOW-2005/blob/release/test-run/etsdrt.bas

With this example there are two formatted 2D arrays (integer and
float), which of course start on different line numbers (4 and 17).
Other parts of this file have a specific format, and are out of scope.
Similar to VRTRawRasterBand, all of other metadata, including grid
size, are provided elsewhere.

Cheers,
-Mike
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Hard coded Pi constants

2016-03-03 Thread Mike Toews 
This is not just a MSVC issue. M_PI and other constants are not part
of the ISO standard, so if you try:
$ gcc --std=c99 small_program_with_M_PI.c
you will get "error: ‘M_PI’ undeclared"

I'm not sure if you can compile GDAL with `--std=c99`, but I don't see
any compelling reason to remove #ifndef blocks for these.

On 4 March 2016 at 06:55, Tanuj Kumar  wrote:
> This is with regard to the bug listed on
> https://trac.osgeo.org/gdal/ticket/6388
>
> The fix is to remove the hard coded pi constants, but apparently, using M_PI
> from math.h as it is breaks compilation on Visual Studio
> http://stackoverflow.com/questions/26065359/m-pi-flagged-as-undeclared-identifier
> (I haven't checked it yet myself; that would need me to switch to windows,
> and I'm on linux right now)
>
> However, we could also declare #define _USE_MATH_CONSTANTS, and then use
> M_PI without declaring it(And not just M_PI, all the other math constants
> listed on https://msdn.microsoft.com/en-us/library/4hwaceh6.aspx).
>
> Could someone please check whether math constants work in later versions of
> visual studio(Both with and without #define _USE_MATH_CONSTANTS)?
>
> Thank you
> Tanuj
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] RFC 61: Call for vote on adoption

2016-02-10 Thread Mike Toews 
On 11 February 2016 at 01:31, Even Rouault  wrote:
> Le mercredi 10 février 2016 13:05:20, Peter Halls a écrit :
>>   ESRI handle this in a non-intuitive way:  XYM is supported, but Z
>> always has a Measure, so is XYZM!   The formal definition is here:
>> https://www.*esri*.com/library/whitepapers/pdfs/*shapefile*.pdf   (1998)
>>
>>   Shape Types having Z are defined on pp19ff where it states:
>>
>> "Shape Types inX,Y,Z Space
>> Shapes of this type have an optional coordinate
>
> Woo, really  That's a genuine scoop.
>
> Unless I'm seriously mistaken, shapelib / OGR has produced XYZ shapefiles
> without M values for more than 20 years. I'm surprised we wouldn't have heard
> about that if such shapefiles couldn't be read. Or perhaps ESRI software is
> robust to missing M too

I've looked into this a few years ago, because I found differences in
XYZ files written by GDAL and Esri software. There are two styles of
(e.g.) PointZ files, where the coordinates are written with either
three or four doubles (the last with "no data" if M is not used). The
style of PointZ is determined by the content length provided in the
record header as 16-bit words (either 14 or 18). Therefore, the specs
are correct when they describe M as an "optional coordinate" (p.15).

For XYZ shapefiles, GDAL always write three doubles, whereas Esri
write four doubles with "no data" as the last coordinate. As a result,
Esri's XYZ shp files are slightly larger than GDAL's.

Occasionally I use software that requires PointZ shapefiles to have
four doubles, so my existing workaround is to use Esri's
CopyFeatures_management tool to rewrite the shapefile.

Could "ogr2ogr -dim 4" be added to this RFC to write XYZ[M] Shapefiles
that always use four floats for coordinates, as Esri software does?


On 11 February 2016 at 02:05, Even Rouault  wrote:
> Would someone looking at this discussion and having access to ESRI software be
> willing to generate tiny (meaning just one single shape, with the smallest
> number of vertices) shapefiles of type PointZ, ArcZ, PolygonZ and MultiPointZ,
> and *without* M values ? So as to enhance the test suite and see in particular
> which value is set exactly as the nodata value in the M component

I can generate these today, and will attach a zip file to
https://trac.osgeo.org/gdal/wiki/MeasuredGeometriesInDrivers
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] RFC 61: Call for vote on adoption

2016-02-10 Thread Mike Toews 
On 11 February 2016 at 02:05, Even Rouault  wrote:
> Would someone looking at this discussion and having access to ESRI software be
> willing to generate tiny (meaning just one single shape, with the smallest
> number of vertices) shapefiles of type PointZ, ArcZ, PolygonZ and MultiPointZ,
> and *without* M values ? So as to enhance the test suite and see in particular
> which value is set exactly as the nodata value in the M component

Done, see 
https://trac.osgeo.org/gdal/attachment/wiki/MeasuredGeometriesInDrivers/

I created two geometry examples per shapefile, and added an "Expected"
attribute to describe the WKT of the expected geometry. While the M
values were editable (despite not requesting this creation option in
ArcCatalog), I left these as "NaN".
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] EPSG:6933

2015-11-10 Thread Mike Toews 
On 11 November 2015 at 02:33, Even Rouault  wrote:
> Long answer:
> This CRS probably didn't exist in the EPSG database v8.5 that was used the
> last time the proj/GDAL derived EPSG database was generated. And I see that
> the ellipsoidal Lamber Cylindrical Equal Area projection method 9835 was also
> unknown of GDAL. I've fixed GDAL to recognize it so next time the proj/GDAL
> EPSG database is regenerated EPSG:6933 will be there.

Are EPSG:6931 and EPSG:6932 supported as well?

The EASE-Grid Data website [1] indicates they are supported by GDAL
with correct version of PROJ.4 (4.8.0), but I can only get the
original EASE-Grid codes working with gdalsrsinfo / ImportFromEPSG.

[1] https://nsidc.org/data/ease/versions.html
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Gdal_calc.py with more than 26 input files

2015-07-16 Thread Mike Toews 
On 17 July 2015 at 11:13, Even Rouault even.roua...@spatialys.com wrote:
 But when you have more than 26 input files, it is dubious that you really want
 to individually identify them with a particular letter/name. As in one of the
 examples, you likely just want to apply a global operation on them, like sum,
 mean,  So they would be better passed as a lis
 t so as to be able to build a 3 dimension numpy array, so you can do things
 like:

 --input foo.tif bar.tif ... --calc rasters.mean(axis=0)

 And you could also still reference a particular raster with rasters[i] syntax

 (To be clear: the above is not something currently available)

It's actually not too far off. For instance, stack several single-band
var*.tif rasters in a multi-band VRT, then calculate the mean along
the 0th axis:

gdalbuildvrt -separate var_stacked.vrt var*.tif
gdal_calc.py -A var_stacked.vrt --outfile=avg_var.tif --calc=A.mean(axis=0)

However, the result avg_var.tif is not as expected due to the use of
blocks. Is this ticket worthy?
-Mike
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] gdal_rasterize with text attributes.

2015-06-27 Thread Mike Toews 
Rasters can only store numeric types, not strings. But you could
either rasterize the primary key, or reclassify the strings into a
so-called lookup table of ID and string (e.g. 1=green, 2=blue,
etc.). And then you could look-up the text from a raster using the
common integer ID.

-Mike

On 27 June 2015 at 11:06, Ricardo Oliveira oliveira.rica...@hotmail.com wrote:
 Hello everyone,

 I'm trying to do a simple vector to raster conversion using an attribute
 field stored as string. Am I wrong or does gdal_rasterize does not accept
 text fields as attributes?

 Thanks for any insight.

 Ricardo O.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Design for sub-second accuracy in OGR ?

2015-04-06 Thread Mike Toews 
On 6 April 2015 at 09:39, Even Rouault even.roua...@spatialys.com wrote:
 I should have mentionned what currently exists indeed :

 struct {
 GInt16  Year;
 GByte   Month;
 GByte   Day;
 GByte   Hour;
 GByte   Minute;
 GByte   Second;
 GByte   TZFlag; /* 0=unknown, 1=localtime(ambiguous),
100=GMT, 104=GMT+1, 80=GMT-5, etc */
 } Date;

In regards to TZFlag, it should also be used to optionally convey
daylight savings information, where known. E.g. see the 'is_dst' in
struct tm [1]. I'm not sure if a GMT value can adequately convey this,
but correct me if it can.

[1] http://www.cplusplus.com/reference/ctime/tm/
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] OGR: How to access Coded Value Domains in File Geodatabases?

2015-01-06 Thread Mike Toews 
Hi Stefan,

I found this issue too and wrote an enhancement ticket:
http://trac.osgeo.org/gdal/ticket/5741

The only workaround appears to be lots of work. The GDB_Items table
(a0004.gdbtable) can be read with OGR, where the code/value pairs
are in attributes of XML code, which can be used to either hash the
values or do an SQL JOIN.

-Mike

On 7 January 2015 at 14:20, Stefan Keller sfkel...@gmail.com wrote:
 Hi,

 ArcGIS offers Coded Value Domains (short: coded domain) to specify
 a valid set of values (code+description) for an attribute [1].

 Any clues on how to access this information using ogr2ogr, i.e. how to
 export this to a separate table when converting e.g. from a ESRI File
 Geodatabase (OpenFileGDB, .gdb directories) to Spatialite or
 GeoPackage?

 Yours, S.

 [1] 
 http://resources.arcgis.com/en/help/main/10.1/index.html#//001s000100
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] gdal_calc.py: produce median raster

2014-11-26 Thread Mike Toews 
Hi Simen,

Try this:

gdal_calc.py -A rgb.tif --outfile ouput.tif --calc median(A, axis=0)

The axis=0 parameter is to perform a median calculation along the
first dimension, which is the band dimension. Lean more here:
http://docs.scipy.org/doc/numpy/reference/generated/numpy.median.html

-Mike

On 24 November 2014 at 23:14, Simen Langseth simlan...@gmail.com wrote:
 I have one GeoTiff file with 3 bands. I want to produce a raster file
 computing pixel wise median value for the three raster bands. How can I do
 that?

 I tried as follows:

 gdal_calc.py -A infile --allBands --outfile outfile --calc median(A)

 But could not get any out file.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Call for discussion on RFC 51: RasterIO() improvements : resampling and progress callback

2014-11-25 Thread Mike Toews 
On 25 November 2014 at 05:34, Even Rouault even.roua...@spatialys.com wrote:
 Otherwise I have not much to say about this RFC. Blurring has two r's,

 Fixed

Gaussian blur also has another complexity, since it's output
dimensions depend on the fftconvolve mode (using scipy.signal docs
http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.fftconvolve.html).

For example, to blur a DTM and maintain the same output dimensions
and geotransform, the mode would need to use a 'valid' convolution
with a padded input array. This padded array would also have a mode
(similarly from numpy docs
http://docs.scipy.org/doc/numpy/reference/generated/numpy.pad.html),
such as 'symmetric'.
http://gis.stackexchange.com/a/10467/1872

Someone else may want to blur a DTM without using a padded array, and
just return a smaller array with modified geotransform.

And a different blur would expand the output dimension and modify
the geotransform using a 'full' fftconvolve mode, using a 'constant'=0
padded array. An example of this is to blur a polygon that barely
extends to the shape's extent, but maintain the blur beyond the
initial extent. The constant value would also need to be supplied, as
a value of 0 for numpy is default.
http://gis.stackexchange.com/a/104338/1872

So in summary, I'm not sure how GDAL could best cover all combinations
of padding and convolution modes, except to offer a few prepared
use-cases, or somehow provide extra parameters. Or just provide a
single use case, similar to Gaussian blur routines in various raster
editing software.

-Mike
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Questions about polygonization

2014-08-19 Thread Mike Toews 
Hi Jack,

I'm not sure if I understand the description of the shapes you see.
Can you provide the WKT for them? I seem to get two polygons that look
like they should:

import numpy as np
import rasterio.features
from shapely.geometry import shape
ar = np.ones((6, 5), 'B') * 11
ar[2, 2] = 12
gt = [0, 1, 0, 6, 0, -1]
for (sh, v) in rasterio.features.shapes(ar, transform=gt):
print('%s: %s' % (v, shape(sh)))

11: POLYGON ((0 6, 0 0, 5 0, 5 6, 0 6), (2 4, 3 4, 3 3, 2 3, 2 4))
12: POLYGON ((2 4, 2 3, 3 3, 3 4, 2 4))

-Mike
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] [BULK] Re: Dissolve shapefile using GDAL/OGR

2014-07-14 Thread Mike Toews 
The path strings are not escaped, so dissolve[13] is '\a' or '\x07'
and not 'a', which makes the path invalid.

For Windows paths, the best practice is to either use forward slashes
(like unix), or use raw escape with r as a string prefix, e.g.

dissolve = r'C:\OSGeo4W64\apps...'

On 12 July 2014 02:21, kalu671 kalu...@gmail.com wrote:
 Hello Dan,
 I tried to use the syntax and I am getting the error:
 The filename, directory name, or volume label syntax is incorrect.

 my syntax is as below:
 inFile = D:\Dir1\Polygons.shp
 outFile = D:\Dir2\Dissolved_polygons.shp

 dissolve = 'C:\OSGeo4W64\apps\saga\saga_cmd shapes_polygons Polygon
 Dissolve -POLYGONS = ' + inFile + ' -DISSOLVED = ' + outFile + ' -FIELD_1=
 Value -DISSOLVE= 0'

 subprocess.call(dissolve, shell=True)

 Would you please suggest what did I do wrong?

 Thanks
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Add to PythonGotchas: saving and closing raster datasets

2014-05-27 Thread Mike Toews 
Thanks Even and Sean for the input. I'll proceed to add this to the
wiki when I have an extended moment.

I think I'll rephrase this gotcha to also apply to OGR datasources,
where the API reads you should always close any opened datasource
with OGR_DS_Destroy() that will ensure all data is correctly
flushed[1]

I think it is worth to mention that raster datasets can use
FlushCache() and vector datasources can use SyncToDisk() to save the
data intermittently without closing, but is driver-specific and not
required before closing.

As for save and close, I think I'll use `del ds` as the recommended
method, with mention that any copies of the references to the
dataset/source also need to be dereferenced.

-Mike

[1] http://www.gdal.org/ogr__api_8h.html#a9d845a6cf6652756925530418905471a
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] compiling without writing permission - problem gdal-1.10.1

2014-05-06 Thread Mike Toews 
On 6 May 2014 21:01, Margherita Di Leo direg...@gmail.com wrote:
 I have a problem compiling gdal-1.10.1 on a Linux system for which I'm not
 administrator.
 In the configure I use the --prefix and --bindir options. Everything goes
 smoothly and it confirms that Libraries have been installed in: the folder
 that I selected (for which i have writing permission). After that, though,
 something goes wrong:
...
 How comes? What am I missing?

Unfortunately this is a bug: http://trac.osgeo.org/gdal/ticket/4563

There are a few workarounds, such as installing the Python bindings
manually with the appropriate --prefix:
$ cd swig/python
$ python setup.py install --prefix=/to/where/you/want/it

or edit swig/python/GNUmakefile (line 71):

$(PYTHON) setup.py install --prefix=$(DESTDIR)$(prefix)

-Mike
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Add to PythonGotchas: saving and closing raster datasets

2014-05-03 Thread Mike Toews 
Hi,

This seems to be a common gotcha: http://gis.stackexchange.com/q/93212/1872

I've been caught by it before, and I'm certain many others have too.
This is a gotcha since WriteArray doesn't write the array to disk,
but rather it is written when the dataset object is dereferenced. I've
drafted up an entry for the wiki, which could be placed before
Certain objects contain a Destroy() method, but you should never use
it.

Furthermore, I'm unsure of the best way to demonstrate how to
dereference Python objects. The wiki has two different styles:
obj = None
or
del obj

I'm unsure of which is regarded best practice, but encourage
consistency. Please edit as necessary or reply with suggestions.

-Mike


=== Saving and closing datasets ===

To save and close raster datasets, the object needs to be
dereferenced, such as setting it to `None`. It is not written with
{{{WriteArray}}}. For example:
{{{
from osgeo import gdal
driver = gdal.GetDriverByName('GTiff')
dst_ds = driver.Create('new.tif', 10, 15)
band = dst_ds.GetRasterBand(1)
arr = band.ReadAsArray()  # raster values are all zero
arr[2,4:] = 50  # modify some data
band.WriteArray(arr)  # raster file still unmodified
band = None  # dereference band to avoid gotcha described previously
dst_ds = None  # save, close
}}}
The last dereference to the raster dataset writes the data
modifications and closes the raster file. The objects may also be
dereferenced using: {{{del band, dst_ds}}}
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] ogr.GeometryTypeToName oddness

2013-04-16 Thread Mike Toews 
With the Python module for GDAL 1.9.2 for Python 2.7, Windows 64-bit,
downloaded from http://www.lfd.uci.edu/~gohlke/pythonlibs/#gdal I see
this oddness:

from osgeo import ogr
for name in dir(ogr):
if name.startswith('wkb'):
i = getattr(ogr, name)
print('%s : %d : %r' % (name, i, ogr.GeometryTypeToName(i)))

wkb25Bit : -2147483648 : 'Unknown (any)'
wkb25DBit : -2147483648 : 'Unknown (any)'
wkbGeometryCollection : 7 : '3D Geometry Collection'
wkbGeometryCollection25D : -2147483641 : 'Geometry Collection'
wkbLineString : 2 : '3D Line String'
wkbLineString25D : -2147483646 : 'Line String'
wkbLinearRing : 101 : 'Unrecognised: 101'
wkbMultiLineString : 5 : '3D Multi Line String'
wkbMultiLineString25D : -2147483643 : 'Multi Line String'
wkbMultiPoint : 4 : '3D Multi Point'
wkbMultiPoint25D : -2147483644 : 'Multi Point'
wkbMultiPolygon : 6 : '3D Multi Polygon'
wkbMultiPolygon25D : -2147483642 : 'Multi Polygon'
wkbNDR : 1 : '3D Point'
wkbNone : 100 : 'None'
wkbPoint : 1 : '3D Point'
wkbPoint25D : -2147483647 : 'Point'
wkbPolygon : 3 : '3D Polygon'
wkbPolygon25D : -2147483645 : 'Polygon'
wkbUnknown : 0 : '3D Unknown (any)'
wkbXDR : 0 : '3D Unknown (any)'

A few observations:

 - The 2D geometries (e.g. wkbLineString) are reported as 3D
geometries by ogr.GeometryTypeToName(), but this is false.
Furthermore, other OGR documents refer to these geometries as so
called 2.5D geometries, and not true 3D geometries.

- The 2.5D geometries (e.g. wkbLineString25D) are simply reported as
their geometry type, and no mention of dimensions.

It seems the 3D  name prepending logic is reversed.

-Mike
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Capture stderr in interactive Python shell

2012-08-05 Thread Mike Toews 
I would like to capture or show GDAL error messages in an interactive
Python shell. For instance, if trying to open an non-existing raster:

from osgeo import gdal
ds = gdal.Open('noexist.tif')

shows nothing in an interactive shell (PythonWin, IDLE), but shows
messages when used in a shell (cmd.exe or Bash):

ERROR 4: `noexist.tif' does not exist in the file system, and is not
recognised as a supported dataset name.

It appears this output is sent to stderr. How can this output also
appear in an interactive Python shell, such as PythonWin or IDLE?

Thanks,
-Mike
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Lengthy GDALRasterAttributeTable output

2011-06-01 Thread Mike Toews 
I'm struggling to understand why both the aux.xml file and output from
gdalinfo are really bulky in file size for my GeoTIFFs. First, output
from gdalinfo on a 375 KB GeoTIFF file has 196664 lines:

Driver: GTiff/GeoTIFF
Files: HorizB.tif
   HorizB.aux
Size is 317, 301
Coordinate System is:
...skip lines, not relevant...
Band 1 Block=317x6 Type=Float32, ColorInterp=Gray
  Min=195.724 Max=842.538
  Minimum=195.724, Maximum=842.538, Mean=440.591, StdDev=127.420
  Metadata:
STATISTICS_MINIMUM=195.72373962402
STATISTICS_MAXIMUM=842.53814697266
STATISTICS_MEAN=440.59141722129
STATISTICS_MEDIAN=434.39051816566
STATISTICS_MODE=559.99110636022
STATISTICS_STDDEV=127.419575971
STATISTICS_HISTONUMBINS=65536
STATISTICS_HISTOMIN=195.72373962402
STATISTICS_HISTOMAX=842.53814697266
LAYER_TYPE=athematic
STATISTICS_HISTOBINVALUES=1|0|0|0|... there are 136875 characters on
line 48 ...0|0|1|
GDALRasterAttributeTable Row0Min=195.7237396240234
BinSize=0.009869755204831507
  FieldDefn index=0
NameHistogram/Name
Type1/Type
Usage0/Usage
  /FieldDefn
  Row index=0
F1/F
  /Row
  Row index=1
F0/F
  /Row
... skip about 196600 lines ...
 Row index=65535
F1/F
  /Row
/GDALRasterAttributeTable

Note: this file was made in RGDAL, and it is just a stand-alone TIFF
file, no aux.xml (I might have deleted it a while ago).

For the second part, I'm using CreateCopy in Python to produce a
similar dataset (same dimensions, also 1 band). The new GeoTIFF file
is 381 KB, but the aux.xml file is nearly 10x larger at 3704 KB
(196632 lines). The driver also produces the output in my console
(yes, all three lines):
Warning 1: Lost metadata writing to GeoTIFF ... too large to fit in tag.
Warning 1: Lost metadata writing to GeoTIFF ... too large to fit in tag.
Warning 1: Lost metadata writing to GeoTIFF ... too large to fit in tag.

And excerpts of the file look like this:
PAMDataset
  PAMRasterBand band=1
GDALRasterAttributeTable Row0Min=195.7237396240234
BinSize=0.009869755204831507
  FieldDefn index=0
NameHistogram/Name
Type1/Type
Usage0/Usage
  /FieldDefn
  Row index=0
F1/F
  /Row
... skip about 196600 lines ...
  Row index=65535
F1/F
  /Row
/GDALRasterAttributeTable
Metadata
  MDI key=STATISTICS_MINIMUM195.72373962402/MDI
  MDI key=STATISTICS_MAXIMUM842.53814697266/MDI
  MDI key=STATISTICS_MEAN440.59141722129/MDI
  MDI key=STATISTICS_MEDIAN434.39051816566/MDI
  MDI key=STATISTICS_MODE559.99110636022/MDI
  MDI key=STATISTICS_STDDEV127.419575971/MDI
  MDI key=STATISTICS_HISTONUMBINS65536/MDI
  MDI key=STATISTICS_HISTOMIN195.72373962402/MDI
  MDI key=STATISTICS_HISTOMAX842.53814697266/MDI
  MDI key=LAYER_TYPEathematic/MDI
  MDI key=STATI|0|1|... there are 136890 characters on line
196629...0|0|1|/MDI
/Metadata
  /PAMRasterBand
/PAMDataset

I can recreate the gdalinfo output with GDAL 1.6.0 and 1.8.0, and I'm
using via Python GDAL version 1.6.0, both on MS Windows. Is there a
way to suppress the creation of the aux.xml file? Is there any way to
fix the GeoTIFF file to stop reporting meaningless histograms?

-Mike
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Adding field to Shapefile via Python

2011-03-20 Thread Mike Toews 
Ok, I found the issue. Rather than adding a new field def'n:
layer_defn.AddFieldDefn(new_field)

it needs to instead be created on the layer:
layer.CreateField(new_field)

The documentation is misleading, as I thought I couldn't add a field
to an existing shapefile if there were features in it. I would
recommend adding a See also link from AddFieldDefn to:
http://www.gdal.org/ogr/ogr__api_8h.html#ab585ef1166c61c4819f7fd46ee4a275

-Mike

On 18 March 2011 18:45, Chaitanya kumar CH chaitanya...@gmail.com wrote:

 Mike,

 A new field cannot be added to a feature definition while there are any 
 features based on that feature definition.
 See the docs for OGRFeatureDefn::AddFieldDefn()
 http://www.gdal.org/ogr/classOGRFeatureDefn.html#40e681d8464b42f1a1fac655f16ac3dd

 On Fri, Mar 18, 2011 at 11:01 AM, Mike Toews mwto...@gmail.com wrote:

 In a Python script, I am updating an exiting shapefile with data by
 adding a column. However, my changes don't seem to be saved. Here is
 what I have:

 from osgeo import ogr
 obs_file = 'myfile.shp'

 source = ogr.Open(obs_file, 1)
 layer = source.GetLayer()
 layer_defn = layer.GetLayerDefn()

 new_field = ogr.FieldDefn('MYFLD', ogr.OFTInteger)
 layer_defn.AddFieldDefn(new_field)

 source = None

 However, when I check the resulting file, I do not see the new field.

 I have also tried inserting the following (before source=None)

 feature = layer.GetNextFeature()
 feature.GetField('EXIST') # an existing float field with data
 feature.SetField('EXIST', 111.1) # works
 feature.SetField('MYFLD', 123) # does nothing
 layer.SetFeature(feature)

 Which successfully updates the data in the existing field, but not the
 added field.

 Where did I go wrong? Thanks in advance.

 -Mike
 ___
 gdal-dev mailing list
 gdal-dev@lists.osgeo.org
 http://lists.osgeo.org/mailman/listinfo/gdal-dev



 --
 Best regards,
 Chaitanya kumar CH.
 /tʃaɪθənjə/ /kʊmɑr/
 +91-9494447584
 17.2416N 80.1426E
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Adding field to Shapefile via Python

2011-03-17 Thread Mike Toews 
In a Python script, I am updating an exiting shapefile with data by
adding a column. However, my changes don't seem to be saved. Here is
what I have:

from osgeo import ogr
obs_file = 'myfile.shp'

source = ogr.Open(obs_file, 1)
layer = source.GetLayer()
layer_defn = layer.GetLayerDefn()

new_field = ogr.FieldDefn('MYFLD', ogr.OFTInteger)
layer_defn.AddFieldDefn(new_field)

source = None

However, when I check the resulting file, I do not see the new field.

I have also tried inserting the following (before source=None)

feature = layer.GetNextFeature()
feature.GetField('EXIST') # an existing float field with data
feature.SetField('EXIST', 111.1) # works
feature.SetField('MYFLD', 123) # does nothing
layer.SetFeature(feature)

Which successfully updates the data in the existing field, but not the
added field.

Where did I go wrong? Thanks in advance.

-Mike
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Polygon topology

2011-03-08 Thread Mike Toews 
You could use the `interiors` length in Shapely:

from osgeo import ogr
from shapely.wkb import loads
source = ogr.Open('my_polygons.shp')
layer = source.GetLayer()
feature = layer.GetNextFeature()
num = 0
while feature:
g = loads(feature.GetGeometryRef().ExportToWkb())
if g.geom_type == 'Polygon':
g.num_holes = len(g.interiors)
elif g.geom_type == 'MultiPolygon':
g.num_holes = sum([len(mp.interiors) for mp in g])
else:
raise Exception('Unexpected geom_type: '+g.geom_type)
print g.geom_type, num, 'with', g.num_holes, 'hole'+('s' if
g.num_holes != 1 else '')
feature = layer.GetNextFeature()
num += 1


-Mike

On 8 March 2011 21:57, Simon Lyngby Kokkendorff sil...@gmail.com wrote:
 Hi List,

   I am using ogr via the python bindings to construct various polygons.
 Here's just a simple question, to which someone might have some input. Is
 there anyway to determine the topology, i.e. the number of holes, of a
 polygon without e.g. having to export to WKT and examining the output
 string. There doesn't seem to be any methods exposed in the bindings to do
 this directly?

   Cheers,
   Simon Kokkendorff,
   National Survey and Cadastre of Denmark


 ___
 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] Polygon topology

2011-03-08 Thread Mike Toews 
On 8 March 2011 23:22, Even Rouault even.roua...@mines-paris.org wrote:
 Not exactly. In fact you have to use the Geometry.GetGeometryCount() that
 returns 1 (the exterior ring) + the number of interior rings. So
 polygon.GetGeometryCount() - 1 should return the number of interior rings

I initially thought so too, except that assumption doesn't work for
MultiPolygon geometry types.

-Mike
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Cannot configure SpatiaLite

2010-07-24 Thread Mike Toews 
Hi all,

I cannot seem to configure GDAL 1.7.2 with SpatiaLite support. My
system is Ubuntu 8.04 Hardy server LTS. I've installed SpatiaLite
version 2.3.1 from source (both lib and tools), and it appears to work
normally. I'm not sure if it matters, but I installed both sqlite3 and
the development library from normal apt packages. I've even run
ldconfig a few times.

Here is my configure command, and some excerpts of the output:

$ ./configure --with-spatialite=/usr/local/lib/libspatialite.so
...
checking for SpatiaLite... checking for spatialite_init in -lspatialite... yes
disabled
checking for SQLite3 library = 3.0.0... yes
...
GDAL is now configured for i686-pc-linux-gnu
...
  SQLite support:yes
  SpatiaLite support:no


I don't understand the output -lspatialite... yes disabled when I
specified the path. Thanks for any help.

-Mike
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Re: Cannot configure SpatiaLite

2010-07-24 Thread Mike Toews 
Ok I figured it out after inspecting the code[1]. This works:

$ ./configure --with-spatialite=/usr/local

or (if SpatiaLite installed by apt in Ubuntu 10.04)

$ ./configure --with-spatialite=/usr

Why isn't the presence of -lspatialite tested and added automatically,
like with all other formats?

-Mike

[1] http://trac.osgeo.org/gdal/browser/trunk/gdal/configure.in#L2088

On 24 July 2010 09:26, Mike Toews mwto...@gmail.com wrote:
 Hi all,

 I cannot seem to configure GDAL 1.7.2 with SpatiaLite support. My
 system is Ubuntu 8.04 Hardy server LTS. I've installed SpatiaLite
 version 2.3.1 from source (both lib and tools), and it appears to work
 normally. I'm not sure if it matters, but I installed both sqlite3 and
 the development library from normal apt packages. I've even run
 ldconfig a few times.

 Here is my configure command, and some excerpts of the output:

 $ ./configure --with-spatialite=/usr/local/lib/libspatialite.so
 ...
 checking for SpatiaLite... checking for spatialite_init in -lspatialite... yes
 disabled
 checking for SQLite3 library = 3.0.0... yes
 ...
 GDAL is now configured for i686-pc-linux-gnu
 ...
  SQLite support:            yes
  SpatiaLite support:        no


 I don't understand the output -lspatialite... yes disabled when I
 specified the path. Thanks for any help.

 -Mike

___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev