[gdal-dev] Question about gdal2wktraster

Tue Jul 14 03:37:50 EDT 2009

Jorge Arévalo wrote:
> Hi Vincent,
>
> 2009/7/13 Vincent Schut <schut at sarvision.nl>:
>   
>> Hi Jorge,
>>
>> I'm not sure what projection your cea.tif was in (it seems you've not pasted
>> the full wkt output from gdalinfo on cea.tif?), but you say in your message
>> that it is supposed to be UTM. For UTM, it is very easy to get the right
>> EPSG code:
>>     
>
> Right. The whole output is:
>
> Driver: GTiff/GeoTIFF
> Files: gsoc09/test_data/tiff_files/cea.tif
> Size is 514, 515
> Coordinate System is:
> PROJCS["unnamed",
>     GEOGCS["NAD27",
>         DATUM["North_American_Datum_1927",
>             SPHEROID["Clarke 1866",6378206.4,294.9786982138982,
>                 AUTHORITY["EPSG","7008"]],
>             AUTHORITY["EPSG","6267"]],
>         PRIMEM["Greenwich",0],
>         UNIT["degree",0.0174532925199433],
>         AUTHORITY["EPSG","4267"]],
>     PROJECTION["Cylindrical_Equal_Area"],
>     PARAMETER["standard_parallel_1",33.75],
>     PARAMETER["central_meridian",-117.333333333333],
>     PARAMETER["false_easting",0],
>     PARAMETER["false_northing",0],
>     UNIT["metre",1,
>         AUTHORITY["EPSG","9001"]]]
> Origin = (-28493.166784412522247,4255884.543802191503346)
> Pixel Size = (60.022136983193739,-60.022136983193739)
> Metadata:
>   AREA_OR_POINT=Area
> Image Structure Metadata:
>   INTERLEAVE=BAND
> Corner Coordinates:
> Upper Left  (  -28493.167, 4255884.544) (117d38'28.21"W, 33d54'13.08"N)
> Lower Left  (  -28493.167, 4224973.143) (117d38'28.21"W, 33d37'30.66"N)
> Upper Right (    2358.212, 4255884.544) (117d18'28.28"W, 33d54'13.08"N)
> Lower Right (    2358.212, 4224973.143) (117d18'28.28"W, 33d37'30.66"N)
> Center      (  -13067.478, 4240428.844) (117d28'28.24"W, 33d45'51.47"N)
> Band 1 Block=514x15 Type=Byte, ColorInterp=Gray
>
>   
Right, I'm really no projections expert, and this could very well be 
another way to describe an UTM projection, but is does not seem to be 
the standard WKT for an UTM zone projection. At least the PROJCS is 
'unnamed', and should otherwise be something like 'UTM 31N', I think... 
Your filename is 'cea' and your PROJECTION string is 'Cylinderical Equal 
Area', while for an UTM projection it should read 'Transverse Mercator'. 
Are you sure the image is in UTM, and not in CEA projection? (For a fast 
overview of projection types, see e.g. 
http://www.colorado.edu/geography/gcraft/notes/mapproj/mapproj.html)
If you're sure it's UTM, you can just use an UTM EPSG code. However, 
have you considered that it might not be you neither your projection 
info that is faulty, but the gdal wktraster driver? As gdalinfo does 
report the correct coordinates on the input file, I think something goes 
wrong either importing the projection stuff into the pg tables, or at 
reading/translating it back into gdal when running gdalinfo on it. I'd 
say, rephrase your question and ask it again at the gdal list. Hopefully 
someone more knowledgeable on wktraster (which I know next to nothing 
about) will have something to say.

> I'm not sure of the UTM zone of my data. How could I can get it? From
> gdalinfo output?
>   
Nope. By googling :-)
Search for 'utm zones', click on the images tab, and choose one you like 
best :-)

Btw, it's considered good practise to not only reply to the sender of an 
email (me), but CC also to the gdal list, so others might chime in 
and/or learn useful stuff from our conversation. As a bonus it will be 
indexed by search engines and you'll be able to find your answer later 
again by searching the web :-)

cheers,
Vincent.
> Best regards
> Jorge
>
>   
>> Vincent.
>>
>> Jorge Arévalo wrote:
>>     
>>> Hello,
>>>
>>> Context: GDAL WKT Raster driver
>>>
>>> I've loaded the tif image
>>> ftp://ftp.remotesensing.org/geotiff/samples/gdal_eg/cea.tif on PostGIS
>>> using gdal2wktraster script using this line: gdal2wktraster.py -r
>>> cea.tif -t usa_mountain_one_band -s 4267 -b 1 -k 514x15 -I -M -o
>>> usa_raster_one_band.sql
>>>
>>> I used the srid 4267 because the information given for "gdalinfo
>>> cea.tif" (Maybe it is my error).
>>>
>>> The problem is: If I execute gdalinfo over the table created in
>>> PostGIS using my driver, I get this:
>>>
>>> (...)
>>> GEOGCS["NAD27",
>>>    DATUM["North_American_Datum_1927",
>>>        SPHEROID["Clarke 1866",6378206.4,294.9786982138982,
>>>            AUTHORITY["EPSG","7008"]],
>>>        AUTHORITY["EPSG","6267"]],
>>>    PRIMEM["Greenwich",0,
>>>        AUTHORITY["EPSG","8901"]],
>>>    UNIT["degree",0.01745329251994328,
>>>        AUTHORITY["EPSG","9122"]],
>>>    AUTHORITY["EPSG","4267"]]
>>> (...)
>>> Corner Coordinates:
>>> Upper Left  (  -28493.168, 4224973.000) (Invalid angle,Invalid angle)
>>> Lower Left  (  -28493.168, 4194061.600) (Invalid angle,Invalid angle)
>>> Upper Right (    2358.210, 4224973.000) (Invalid angle,Invalid angle)
>>> Lower Right (    2358.210, 4194061.600) (Invalid angle,Invalid angle)
>>> Center      (  -13067.479, 4209517.300) (Invalid angle,Invalid angle)
>>> (...)
>>>
>>> Look at "Invalid angle error". If I execute gdalinfo over the tif
>>> image, I get this:
>>> (...)
>>> PROJCS["unnamed",
>>>    GEOGCS["NAD27",
>>>        DATUM["North_American_Datum_1927",
>>>            SPHEROID["Clarke 1866",6378206.4,294.9786982138982,
>>>                AUTHORITY["EPSG","7008"]],
>>>            AUTHORITY["EPSG","6267"]],
>>>        PRIMEM["Greenwich",0],
>>>        UNIT["degree",0.0174532925199433],
>>>        AUTHORITY["EPSG","4267"]],
>>>    PROJECTION["Cylindrical_Equal_Area"],
>>>    PARAMETER["standard_parallel_1",33.75],
>>>    PARAMETER["central_meridian",-117.333333333333],
>>>    PARAMETER["false_easting",0],
>>>    PARAMETER["false_northing",0],
>>>    UNIT["metre",1,
>>>        AUTHORITY["EPSG","9001"]]]
>>> (...)
>>> Corner Coordinates:
>>> Upper Left  (  -28493.167, 4255884.544) (117d38'28.21"W, 33d54'13.08"N)
>>> Lower Left  (  -28493.167, 4224973.143) (117d38'28.21"W, 33d37'30.66"N)
>>> Upper Right (    2358.212, 4255884.544) (117d18'28.28"W, 33d54'13.08"N)
>>> Lower Right (    2358.212, 4224973.143) (117d18'28.28"W, 33d37'30.66"N)
>>> Center      (  -13067.478, 4240428.844) (117d28'28.24"W, 33d45'51.47"N)
>>> (...)
>>>
>>> The important thing is the "UNIT["metre",1,AUTHORITY["EPSG","9001"]]]
>>> part of the string that doesn't appear in my driver. For this reason,
>>> the transformation between utm (the original coords of the tif image)
>>> and lat,long doesn't take place. Then, the UTM coordinates are
>>> considered as invalid angles, of course.
>>>
>>> I think that my great mistake, as I said, was to choose 4267 as srid
>>> (I wasn't sure). Am I right? If yes, What srid should I use? If not,
>>> maybe the error is in my driver's code...
>>>
>>> Many thanks in advance
>>> Best regards,
>>>
>>> Jorge
>>> _______________________________________________
>>> gdal-dev mailing list
>>> gdal-dev at lists.osgeo.org
>>> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>>>
>>>       
>>     