Ticket #1471 (assigned defect)

Opened 1 year ago

Last modified 10 months ago

gdalwarp should keep georeferencing info

Reported by: matt.wilkie@gov.yk.ca Assigned to: mloskot (accepted)
Priority: normal Milestone:
Component: default Version: 1.4.0
Severity: normal Keywords:
Cc: warmerdam

Description (Last modified by mloskot)

Currently gdalwarp drops the projection info if it is not specifically assigned. I think it should either keep all projection info unless directed otherwise, or, have an option to retain it.

Eg. 

gdalwarp inimg.tif outimg.tif ---> keeps prj
gdalwarp -t_srs SAME inimg.tif outimg.tif ---> keeps prj
gdalwarp -t_srs NONE inimg.tif outimg.tif ---> drops prj
gdalwarp -t_srs [a proj4 string] inimg.tif outimg.tif ---> transforms to new prj

Attachments

test.tif (1.0 MB) - added by matt.wilkie@gov.yk.ca on 02/06/07 12:40:25.
source image
gw-out.tif (1.0 MB) - added by matt.wilkie@gov.yk.ca on 02/06/07 12:41:46.
gdalwarp output with no georeferencing
fwtools.cmd (486 bytes) - added by matt.wilkie@gov.yk.ca on 02/09/07 15:00:55.
windows cmd script to run fwtools version passed on command line

Change History

02/02/07 17:26:39 changed by warmerdam 

Matt,

As far as I can see, gdalwarp does preserve the coordinate system by default.

warmerda@amd64[17]% gdalwarp utm.tif out.tif
Creating output file that is 512P x 512L.
Processing input file utm.tif.
:0...10...20...30...40...50...60...70...80...90...100 - done.
warmerda@amd64[18]% gdalinfo out.tif 
Driver: GTiff/GeoTIFF
Size is 512, 512
Coordinate System is:
PROJCS["NAD27 / UTM zone 11N",
    GEOGCS["NAD27",
        DATUM["North_American_Datum_1927",
            SPHEROID["Clarke 1866",6378206.4,294.9786982139006,
                AUTHORITY["EPSG","7008"]],
            AUTHORITY["EPSG","6267"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4267"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-117],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
...

Could you provide more detail on the situations in which the coordinate
system is being lost?

02/06/07 12:38:52 changed by matt.wilkie@gov.yk.ca 

K:\work\Env-dat.013\Albers\test>
 gdalinfo test.tif
Driver: GTiff/GeoTIFF
Size is 581, 571
Coordinate System is:
PROJCS["unnamed",
    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["Albers_Conic_Equal_Area"],
    PARAMETER["standard_parallel_1",61.666666667],
    PARAMETER["standard_parallel_2",68],
    PARAMETER["latitude_of_center",59],
    PARAMETER["longitude_of_center",-132.5],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",500000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (751020.239295186470000,630043.162140076170000)
Pixel Size = (211.949519683800250,-211.842763459462620)
Metadata:
  AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  (  751020.239,  630043.162) (127d59'33.82"W, 60d 5'22.88"N)
Lower Left  (  751020.239,  509080.944) (128d 8'30.43"W, 59d 0'14.59"N)
Upper Right (  874162.910,  630043.162) (125d47'43.31"W, 59d59'30.87"N)
Lower Right (  874162.910,  509080.944) (126d 0'58.42"W, 58d54'33.63"N)
Center      (  812591.575,  569562.053) (126d59'11.28"W, 59d30'13.33"N)
Band 1 Block=581x14 Type=Byte, ColorInterp=Red
  NoData Value=255
Band 2 Block=581x14 Type=Byte, ColorInterp=Green
  NoData Value=255
Band 3 Block=581x14 Type=Byte, ColorInterp=Blue
  NoData Value=255

K:\work\Env-dat.013\Albers\test>
 gdalwarp test.tif gw-out.tif
Creating output file that is 581P x 571L.
Processing input file test.tif.
:0...10...20...30...40...50...60...70...80...90...100 - done.

K:\work\Env-dat.013\Albers\test>
 gdalinfo gw-out.tif
Driver: GTiff/GeoTIFF
Size is 581, 571
Coordinate System is `'
Origin = (751020.239295186470000,630043.162140076170000)
Pixel Size = (211.897074926247290,-211.897074926247290)
Corner Coordinates:
Upper Left  (  751020.239,  630043.162)
Lower Left  (  751020.239,  509049.932)
Upper Right (  874132.440,  630043.162)
Lower Right (  874132.440,  509049.932)
Center      (  812576.340,  569546.547)
Band 1 Block=581x14 Type=Byte, ColorInterp=Red
Band 2 Block=581x14 Type=Byte, ColorInterp=Green
Band 3 Block=581x14 Type=Byte, ColorInterp=Blue

02/06/07 12:40:25 changed by matt.wilkie@gov.yk.ca

  • attachment test.tif added.

source image

02/06/07 12:41:46 changed by matt.wilkie@gov.yk.ca

  • attachment gw-out.tif added.

gdalwarp output with no georeferencing

02/06/07 15:09:09 changed by warmerdam 

Matt,

I made a quick try with this file in my local dev setup, and it carries
over the coordinate system ok.  So now I am starting to suspect something
to do with software version or the environment (ie. finding support .csv files
and so forth).  Could you be specific about what version of GDAL you are 
using, and how it was built?  Perhaps you could try the same thing with 
FWTools?

02/06/07 17:17:41 changed by matt.wilkie@gov.yk.ca 

When the bug report was initially filed I used "GDAL 1.4.0.0, FWTools 1.1.4, released 2007/01/16" (windows xp). I just tried 1.20 with the same results. I also tried opening fwtools shell and deleting all environment variables except:

ALLUSERSPROFILE=C:\Documents and Settings\All Users
APPDATA=C:\Documents and Settings\mattw\Application Data
FWTOOLS_DIR=c:\local\fwtools
GDAL_DATA=c:\local\fwtools\data
GDAL_DRIVER_PATH=c:\local\fwtools\gdal_plugins
GEOTIFF_CSV=c:\local\fwtools\data
Path=c:\local\fwtools\bin;c:\local\fwtools\python

...and then deleted all of these except for PATH, to no avail.

contents of c:\local\fwtools\data:

2005/07/12  05:50 PM           131,088 alias.csv
2005/07/12  05:50 PM           285,199 area.csv
2005/07/12  05:50 PM           150,002 change.csv
2005/07/12  05:50 PM             8,511 codes.csv
2005/07/12  05:50 PM             2,134 coordinate_axis.csv
2005/07/12  05:50 PM             6,389 coordinate_axis_name.csv
2005/07/12  05:50 PM           534,293 coordinate_operation.csv
2005/07/12  05:50 PM            90,617 coordinate_operation_method.csv
2005/07/12  05:50 PM            19,499 coordinate_operation_parameter.csv
2005/07/12  05:50 PM           245,852 coordinate_operation_parameter_value.csv
2005/07/12  05:50 PM             4,419 coordinate_operation_path.csv
2005/07/12  05:50 PM           595,129 coordinate_reference_system.csv
2005/07/12  05:50 PM             7,863 coordinate_system.csv
2007/01/15  09:41 PM            11,328 cubewerx_extra.wkt
2005/07/12  05:50 PM            74,102 datum.csv
2005/07/12  05:50 PM            33,612 deprecation.csv
2007/01/15  09:41 PM           345,488 ecw_cs.dat
2007/01/15  09:41 PM            12,335 ellipsoid.csv
2007/01/15  09:41 PM                50 epsg.wkt
2007/01/15  09:41 PM           206,014 esri_extra.wkt
2007/01/15  09:41 PM            36,097 gcs.csv
2007/01/15  09:41 PM             2,037 gdalicon.png
2007/01/15  09:41 PM            13,084 GDALLogoBW.svg
2007/01/15  09:41 PM            12,367 GDALLogoColor.svg
2007/01/15  09:41 PM            12,367 GDALLogoGS.svg
2007/01/15  09:41 PM           150,148 gdal_datum.csv
2005/07/12  05:50 PM             1,514 naming_system.csv
2007/01/15  09:41 PM           400,251 pcs.csv
2007/01/15  09:41 PM             1,599 prime_meridian.csv
2007/01/15  09:41 PM           164,847 projop_wparm.csv
2007/01/15  09:41 PM            13,304 s57agencies.csv
2007/01/15  09:41 PM             7,254 s57attributes.csv
2007/01/15  09:41 PM            18,115 s57attributes_aml.csv
2007/01/15  09:41 PM             8,745 s57attributes_iw.csv
2007/01/15  09:41 PM            20,885 s57expectedinput.csv
2007/01/15  09:41 PM            31,212 s57objectclasses.csv
2007/01/15  09:41 PM            63,970 s57objectclasses_aml.csv
2007/01/15  09:41 PM            37,016 s57objectclasses_iw.csv
2007/01/15  09:41 PM             9,216 seed_2d.dgn
2007/01/15  09:41 PM             2,048 seed_3d.dgn
2007/01/15  09:41 PM            10,315 stateplane.csv
2007/01/15  09:41 PM            18,366 unit_of_measure.csv
2005/07/12  05:50 PM             1,524 version_history.csv
              43 File(s)      3,800,205 bytes

02/09/07 14:59:12 changed by matt.wilkie@gov.yk.ca 

Something changed between fwtools 1.0.6 and 1.0.7. I installed (simultaneously) win32 fwtools 1.0.6 through 1.1.0 and they all fail to carry over the georeferencing except 1.0.6.

------
 fwtools.cmd 1.0.6
Attempting to start FWTools 1.0.6 shell
C:\WINDOWS\SYSTEM32\cmd.exe /D /T:80 /K ".\FWTools1.0.6\setfw.bat"
E:\scratch>
gdalinfo --version
GDAL 1.3.2.0, FWTools 1.0.6, released 2006/08/24

E:\scratch>
cd \temp\gwarp-no-georef\

E:\temp\gwarp-no-georef>
gdalwarp utm.tif v1.0.6.tif
Creating output file that is 1101P x 1084L.
Processing input file utm.tif.
:0...10...20...30...40...50...60...70...80...90...100 - done.
Open GDAL Datasets:
  1 N DriverIsNULL 512x512x0

E:\temp\gwarp-no-georef>
gdalinfo v1.0.6.tif
Driver: GTiff/GeoTIFF
Size is 1101, 1084
Coordinate System is:
PROJCS["NAD83 / UTM zone 9N",
    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["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-129],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","26909"]]
Origin = (555661.751911,6655216.415796)
Pixel Size = (105.86363614,-105.86363614)
Metadata:
  AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  (  555661.752, 6655216.416) (128d 0'4.09"W, 60d 1'49.42"N)
Lower Left  (  555661.752, 6540460.234) (128d 1'52.12"W, 59d 0'0.16"N)
Upper Right (  672217.615, 6655216.416) (125d54'42.70"W, 59d59'53.03"N)
Lower Right (  672217.615, 6540460.234) (126d 0'16.11"W, 58d58'8.44"N)
Center      (  613939.684, 6597838.325) (126d59'13.69"W, 59d30'12.42"N)
Band 1 Block=1101x7 Type=Byte, ColorInterp=Red
Band 2 Block=1101x7 Type=Byte, ColorInterp=Green
Band 3 Block=1101x7 Type=Byte, ColorInterp=Blue

E:\temp\gwarp-no-georef>
exit

E:\scratch>fwtools.cmd 1.0.7
Attempting to start FWTools 1.0.7 shell
C:\WINDOWS\SYSTEM32\cmd.exe /D /T:80 /K ".\FWTools1.0.7\setfw.bat"
E:\scratch>
 cd \temp\gwarp-no-georef\

E:\temp\gwarp-no-georef>
 gdalwarp utm.tif v1.0.7.tif
Creating output file that is 1101P x 1084L.
Processing input file utm.tif.
:0...10...20...30...40...50...60...70...80...90...100 - done.
Open GDAL Datasets:
  1 N DriverIsNULL 512x512x0

E:\temp\gwarp-no-georef>
 gdalinfo v1.0.7.tif
Driver: GTiff/GeoTIFF
Size is 1101, 1084
Coordinate System is `'
Origin = (555661.751911,6655216.415796)
Pixel Size = (105.86363614,-105.86363614)
Corner Coordinates:
Upper Left  (  555661.752, 6655216.416)
Lower Left  (  555661.752, 6540460.234)
Upper Right (  672217.615, 6655216.416)
Lower Right (  672217.615, 6540460.234)
Center      (  613939.684, 6597838.325)
Band 1 Block=1101x7 Type=Byte, ColorInterp=Red
Band 2 Block=1101x7 Type=Byte, ColorInterp=Green
Band 3 Block=1101x7 Type=Byte, ColorInterp=Blue

E:\temp\gwarp-no-georef>
 gdalinfo --version
GDAL 1.3.2.0, FWTools 1.0.7, released 2006/10/20

02/09/07 15:00:55 changed by matt.wilkie@gov.yk.ca

  • attachment fwtools.cmd added.

windows cmd script to run fwtools version passed on command line

03/26/07 21:04:06 changed by warmerdam

  • status changed from assigned to new.
  • description changed.
  • cc set to warmerdam.
  • priority changed from highest to normal.
  • owner changed from warmerdam to mloskot.
  • severity changed from blocker to normal.

Mateusz,

Can you try and reproduce this?

09/13/07 12:28:58 changed by mloskot

  • status changed from new to assigned.
  • description changed.

09/13/07 14:17:06 changed by mloskot

I can not reproduce this problem with GDAL from FWTools 1.3.4 (under Windows) nor with GDAL from SVN (under Linux). Here are steps I've taken to test it:

  • warping test.tif 
    D:\dev\gdal\bugs\1471>gdalwarp test.tif out.tif
Creating output file that is 581P x 571L.
Processing input file test.tif.
:0...10...20...30...40...50...60...70...80...90...100 - done.
  • info reported for the output file 
    D:\dev\gdal\bugs\1471>gdalinfo out.tif
Driver: GTiff/GeoTIFF
Files: out.tif
Size is 581, 571
Coordinate System is:
PROJCS["unnamed",
    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["Albers_Conic_Equal_Area"],
    PARAMETER["standard_parallel_1",61.666666667],
    PARAMETER["standard_parallel_2",68],
    PARAMETER["latitude_of_center",59],
    PARAMETER["longitude_of_center",-132.5],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",500000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (751020.239295186470000,630043.162140076170000)
Pixel Size = (211.897074926247290,-211.897074926247290)
Metadata:
  AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  (  751020.239,  630043.162) (127d59'33.82"W, 60d 5'22.88"N)
Lower Left  (  751020.239,  509049.932) (128d 8'30.56"W, 59d 0'13.58"N)
Upper Right (  874132.440,  630043.162) (125d47'45.26"W, 59d59'30.98"N)
Lower Right (  874132.440,  509049.932) (126d 1'0.50"W, 58d54'32.73"N)
Center      (  812576.340,  569546.547) (126d59'12.33"W, 59d30'12.87"N)
Band 1 Block=581x14 Type=Byte, ColorInterp=Red
Band 2 Block=581x14 Type=Byte, ColorInterp=Green
Band 3 Block=581x14 Type=Byte, ColorInterp=Blue
  • original file 
    D:\dev\gdal\bugs\1471>gdalinfo test.tif
Driver: GTiff/GeoTIFF
Files: test.tif
Size is 581, 571
Coordinate System is:
PROJCS["unnamed",
    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["Albers_Conic_Equal_Area"],
    PARAMETER["standard_parallel_1",61.666666667],
    PARAMETER["standard_parallel_2",68],
    PARAMETER["latitude_of_center",59],
    PARAMETER["longitude_of_center",-132.5],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",500000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (751020.239295186470000,630043.162140076170000)
Pixel Size = (211.949519683800250,-211.842763459462620)
Metadata:
  AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  (  751020.239,  630043.162) (127d59'33.82"W, 60d 5'22.88"N)
Lower Left  (  751020.239,  509080.944) (128d 8'30.43"W, 59d 0'14.59"N)
Upper Right (  874162.910,  630043.162) (125d47'43.31"W, 59d59'30.87"N)
Lower Right (  874162.910,  509080.944) (126d 0'58.42"W, 58d54'33.63"N)
Center      (  812591.575,  569562.053) (126d59'11.28"W, 59d30'13.33"N)
Band 1 Block=581x14 Type=Byte, ColorInterp=Red
  NoData Value=255
Band 2 Block=581x14 Type=Byte, ColorInterp=Green
  NoData Value=255
Band 3 Block=581x14 Type=Byte, ColorInterp=Blue
  NoData Value=255

The only difference is that pixel size in the original file (test.tif) is a square size: 

Pixel Size = (211.949519683800250,-211.842763459462620)

So, gdalwarp corrects the size to square in the out.tif file, according to the semantic of GDALSuggestedWarpOutput function: 

Pixel Size = (211.897074926247290,-211.897074926247290)

I'm not sure if it's reasonable to test old GDAL versions. Is it?

09/14/07 16:53:55 changed by maphew

Today I can't reproduce the problem with either fwtools 1.3.4 or 1.3.6 on windows. I don't know what I'm doing differently. I guess I'll just have to wait for it show it's face again and then study it right away!