Python-NG Bindings -Complex data ReadAsArray

[tisham@…]

Reading complex datasets form Python using the Python-NG bindings produces the following error:

Wow complex
Traceback (most recent call last):
  File "testwin.py", line 17, in <module>
    databuffer = band.ReadAsArray(0,0,10,10,10,10)
  File "c:\Python25\Lib\site-packages\gdal.py", line 657, in ReadAsArray
    buf_xsize, buf_ysize, buf_obj )
  File "c:\Python25\Lib\site-packages\gdal_array.py", line 135, in BandReadAsArray
    ar = numpy.reshape(ar, [1,win_ysize,win_xsize])
  File "c:\python25\lib\site-packages\numpy-1.0.1-py2.5-win32.egg\numpy\core\fromnumeric.py", line 62, in reshape
    return reshape(newshape, order=order)
ValueError: total size of new array must be unchanged

This might be dependent on the numpy-version being used. This bug has been reproduced in Windows with the MingW and VC++ Express compilers and in Linux-Debian with stock numpy install from apt-get.

[warmerdam]

I believe there is a fundamental problem here with how complex values 
are handled by the numpy/gdal_array stuff.   The issue, I think, is that
numpy doesn't have a data type for complex values, so they essentially have
to be treated as a "doubled up" array of the base type (ie. CInt16 needs
to be treated as a doubled up array of Int16).  

This will require some careful thought. 

[warmerdam]

I have tried this script:

import gdal

ds = gdal.Open('cint_sar.tif')
x = ds.ReadAsArray()

print x

With the old bindings (fwtools) and get:

warmerda@amd64[54]% pkg/FWTools-1.1.3/bin_safe/python test_complex.py

GDAL: GDALOpen(cint_sar.tif) succeeds as GTiff.
[[  35. -32.j   -8. -12.j   25. -21.j   52.  -1.j  -47. +14.j]
 [  -1.  -3.j   24. +36.j   39.  +6.j   29. +21.j  -12. +59.j]
 [ -22. +54.j  -12. +25.j   11. -63.j  116. -16.j   63. +77.j]
 [   4. +20.j   27. +35.j  -25. +43.j   48. -35.j   83. +26.j]
 [  -4. +18.j   77. +60.j   84.+122.j   36. +15.j   33. -15.j]
 [  47. -44.j  127. -17.j  173. +29.j   48. +16.j   12. -35.j]]

So this did used to work. Still trying to find an environment with the new bindings where I can test, but I think my original guess was wrong.

[warmerdam]

The core problem was a mistaken assumption that numpy.complex64 was the same as gdalcost.GDT_CFloat64. But CFloat64 means two 64bit floats, while complex64 means a complex number totaling 64bits.

This patch fixes the problem.

Index: gdal_array.py =================================================================== --- gdal_array.py (revision 11104) +++ gdal_array.py (working copy) @@ -18,10 +18,10 @@

gdalconst.GDT_Int32 : numpy.int32, gdalconst.GDT_Float32 : numpy.float32, gdalconst.GDT_Float64 : numpy.float64,

• gdalconst.GDT_CInt16 : numpy.complex,
• gdalconst.GDT_CInt32 : numpy.complex,
• gdalconst.GDT_CFloat32 : numpy.complex,
• gdalconst.GDT_CFloat64 : numpy.complex64

+ gdalconst.GDT_CInt16 : numpy.complex64, + gdalconst.GDT_CInt32 : numpy.complex64, + gdalconst.GDT_CFloat32 : numpy.complex64, + gdalconst.GDT_CFloat64 : numpy.complex128

}

def GetArrayFilename( array ):

@@ -42,6 +42,11 @@

def flip_code(code):

if isinstance(code, type):

+ # since several things map to complex64 we must carefully select + # the opposite that is an exact match (ticket 1518) + if code == numpy.complex64: + return gdalconst.GDT_CFloat32 +

for key, value in codes.items():

if value == code:

return key

There are also some other fixes required to get gdal.Dataset.ReadAsArray() working properly. It seems it is not used in the test suite.

[warmerdam]

(Repeating last with better formatting)

The core problem was a mistaken assumption that numpy.complex64 was the same as gdalcost.GDT_CFloat64. But CFloat64 means two 64bit floats, while complex64 means a complex number totaling 64bits.

This patch fixes the problem.

Index: gdal_array.py
===================================================================
--- gdal_array.py	(revision 11104)
+++ gdal_array.py	(working copy)
@@ -18,10 +18,10 @@
             gdalconst.GDT_Int32     :   numpy.int32,
             gdalconst.GDT_Float32   :   numpy.float32,
             gdalconst.GDT_Float64   :   numpy.float64,
-            gdalconst.GDT_CInt16    :   numpy.complex,
-            gdalconst.GDT_CInt32    :   numpy.complex,
-            gdalconst.GDT_CFloat32  :   numpy.complex,
-            gdalconst.GDT_CFloat64  :   numpy.complex64
+            gdalconst.GDT_CInt16    :   numpy.complex64,
+            gdalconst.GDT_CInt32    :   numpy.complex64,
+            gdalconst.GDT_CFloat32  :   numpy.complex64,
+            gdalconst.GDT_CFloat64  :   numpy.complex128
         }
 
 def GetArrayFilename( array ):
@@ -42,6 +42,11 @@
     
 def flip_code(code):
     if isinstance(code, type):
+        # since several things map to complex64 we must carefully select
+        # the opposite that is an exact match (ticket 1518)
+        if code == numpy.complex64:
+            return gdalconst.GDT_CFloat32
+        
         for key, value in codes.items():
             if value == code:
                 return key

There are also some other fixes required to get gdal.Dataset.ReadAsArray() working properly. It seems it is not used in the test suite.

[warmerdam]

Patched in trunk (r11110).

Checked for now by autotest/gcore/numpy_rw.py.

Still needs to be migrated into 1.4 branch.

[warmerdam]

Changes, including test, now in 1.4 branch.