[gdal-dev] python gdal gdal.ReprojectImage produces no data

[Jachym Cepicky]

Hi

I would like to write "custom" version of the gdalwarp tool in Python
using gdal.ReprojectImage() function, but the only result I get, is just
a empty "black" image.

The code would look like this:

# -------     python code start  -------
from osgeo import gdal

# let's have georeferenced geotiff
inds = gdal.Open("filename.tif")

# it *is* georeferenced
inds.GetProjection()
# 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.2572235629972,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]]'

# create the driver
driver = gdal.GetDriverByName("GTiff")

# create new output file
outds = driver.Create("outfile.tif",inds.RasterXSize, inds.RasterYSize)

# do the transformation - actually no transformation at all
gdal.ReprojectImage(inds, outds, inds.GetProjection(), inds.GetProjection())

# close
outds = None
inds = None
# -------     python code end  -------

As result, the new created file should be the same, like to original
file was.

When I run gdalinfo on the original file, I get 

Driver: GTiff/GeoTIFF
Files: /tmp/elevationwgs.tif
Size is 709, 386
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.2572235629972,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-103.871109718285695,44.501735270606410)
Pixel Size = (0.000340671212621,-0.000340671212621)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-103.8711097,  44.5017353) (103d52'15.99"W, 44d30'6.25"N)
Lower Left  (-103.8711097,  44.3702362) (103d52'15.99"W, 44d22'12.85"N)
Upper Right (-103.6295738,  44.5017353) (103d37'46.47"W, 44d30'6.25"N)
Lower Right (-103.6295738,  44.3702362) (103d37'46.47"W, 44d22'12.85"N)
Center      (-103.7503418,  44.4359857) (103d45'1.23"W, 44d26'9.55"N)
Band 1 Block=709x5 Type=UInt16, ColorInterp=Palette
  Color Table (RGB with 65536 entries)
    0: 255,255,255,255
    1: 255,255,255,255
    [...]
    65535: 0,0,0,255

When I run gdalinfo on the new created file file, I get 

Driver: GTiff/GeoTIFF
Files: /tmp/pokus12
Size is 709, 386
Coordinate System is `'
Image Structure Metadata:
  INTERLEAVE=BAND
  Corner Coordinates:
  Upper Left  (    0.0,    0.0)
  Lower Left  (    0.0,  386.0)
  Upper Right (  709.0,    0.0)
  Lower Right (  709.0,  386.0)
  Center      (  354.5,  193.0)
  Band 1 Block=709x11 Type=Byte, ColorInterp=Gray


That was all :-( Projection and coordinates can be set using
dsout.SetProjection and dsout.SetGeoTransform - that would not be a
problem. The problem are the data, which were not copied.

I'm I missing something? Could this be a problem of my GDAL
installation?

gdalinfo --version
GDAL 1.6.1, released 2009/05/11

(gdal.__version__  is not existing in python)

Thanks

Jachym
-- 
Jachym Cepicky
e-mail: jachym.cepicky gmail com
URL: http://les-ejk.cz
GPG: http://www.les-ejk.cz/pgp/JachymCepicky.pgp
Key fingerprint: 0C6D 0EAE 76BD 506C F299  ED8A C8AB 74B8 08D4 E08F
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 197 bytes
Desc: Digital signature
Url : http://lists.osgeo.org/pipermail/gdal-dev/attachments/20090729/961437bc/attachment.bin