Opened 16 years ago

Closed 16 years ago

#2463 closed defect (invalid)

projwin fails when asking reversed rasters

Reported by: regodon Owned by: warmerdam
Priority: normal Milestone:
Component: Utilities Version: unspecified
Severity: normal Keywords:
Cc:

Description

gdal_translate documentation says projwin parameters are:

-projwin ulx uly lrx lry:
    Selects a subwindow from the source image for copying (like -srcwin) but with the corners given in georeferenced coordinates.

We are having problems with "reverse" rasters which start on lower left corner and end at upper right corner (with a positive Y steps instead of negative ones)

# /usr/local/FWTools-2.0.6/bin_safe/gdal_translate -of AAIgrid -projwin 108882 4842103 110382 4843703   meteologica.tiff /tmp/meteologica.txt
Input file size is 2390, 1750
Computed -srcwin 346 37 3 -2 from projected window.
-srcwin 346 37 3 -2 falls outside raster size of 2390x1750
or is otherwise illegal.

It seems -projwin is trying to convert to cells and getting negative values and failing.

This is our raster:

[root@gregal FWTools-2.0.6]# /usr/local/FWTools-2.0.6/bin_safe/gdalinfo meteologica.tiff
Driver: GTiff/GeoTIFF
Files: meteologica.tiff
Size is 2390, 1750
Coordinate System is:
PROJCS["UTM Zone 30, Northern Hemisphere",
    GEOGCS["International 1909 (Hayford)",
        DATUM["unknown",
            SPHEROID["unnamed",6378388,297.000000000005]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-3],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-64500.000000000000000,4861000.000000000000000)
Pixel Size = (500.000000000000000,-500.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  -64500.000, 4861000.000) ( 10d 0'10.81"W, 43d41'10.38"N)
Lower Left  (  -64500.000, 3986000.000) (  9d14'51.84"W, 35d51'18.04"N)
Upper Right ( 1130500.000, 4861000.000) (  4d48'52.79"E, 43d37'59.93"N)
Lower Right ( 1130500.000, 3986000.000) (  3d58'24.96"E, 35d48'53.65"N)
Center      (  533000.000, 4423500.000) (  2d36'49.11"W, 39d57'37.14"N)
Band 1 Block=2390x1 Type=Int16, ColorInterp=Gray
  NoData Value=-9999

If we change ymin and ymax, gdal_traslate works fine with this raster, but fails with other rasters. Since we have lots of rasters and we don't know if they are from north to south or viceversa, we can't guess how to call projwin.

Change History (4)

comment:1 by warmerdam, 16 years ago

Resolution: invalid
Status: newclosed

Francisco,

The raster you show is a normal "north up" raster but the specified -projwin argument has a upper-left-y value that is south of the lower-right-y. This is selecting a flipped window which is not supported by gdal_translate.

If some of your rasters are north-up and some are south-up, you are just going to need to detect this in your processing and specify projected windows accordingly.

sorry

comment:2 by warmerdam, 16 years ago

Component: defaultUtilities

comment:3 by regodon, 16 years ago

Resolution: invalid
Status: closedreopened

I'm sorry, you are right. I included the wrong raster. This is the right one:

# /usr/local/FWTools-2.0.6/bin_safe/gdal_translate -of AAIgrid -projwin -8.254 43.083 -8.204 43.033 predicciones5000-hirlam-precipitacion-20080710-h18.tiff /tmp/meteologica.txt
Input file size is 281, 171
Computed -srcwin 25 152 1 0 from projected window.
-srcwin 25 152 1 0 falls outside raster size of 281x171
or is otherwise illegal.

We generated this GeoTiff from a raw+vrt. This is the VRT file:

<VRTDataset rasterXSize="281" rasterYSize="171">
  <SRS>+proj=latlong +ellps=intl</SRS>
  <GeoTransform>-9.525, 0.05, 0.0, 35.475, 0.0, 0.05</GeoTransform>
  <VRTRasterBand dataType="Int16" band="1" subClass="VRTRawRasterBand">
          <SourceFilename relativeToVRT=\"1\">file.raw</SourceFilename>
  </VRTRasterBand>
</VRTDataset>

Gdalinfo of generated tiff:

# /usr/local/FWTools-2.0.6/bin_safe/gdalinfo predicciones5000-hirlam-precipitacion-20080710-h18.tiff
Driver: GTiff/GeoTIFF
Files: predicciones5000-hirlam-precipitacion-20080710-h18.tiff
Size is 281, 171
Coordinate System is:
GEOGCS["International 1909 (Hayford)",
    DATUM["unknown",
        SPHEROID["unnamed",6378388,297.000000000005]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]
Origin = (-9.525000000000000,35.475000000000001)
Pixel Size = (0.050000000000000,0.050000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=DEFLATE
Corner Coordinates:
Upper Left  (  -9.5250000,  35.4750000) (  9d31'30.00"W, 35d28'30.00"N)
Lower Left  (  -9.5250000,  44.0250000) (  9d31'30.00"W, 44d 1'30.00"N)
Upper Right (   4.5250000,  35.4750000) (  4d31'30.00"E, 35d28'30.00"N)
Lower Right (   4.5250000,  44.0250000) (  4d31'30.00"E, 44d 1'30.00"N)
Center      (  -2.5000000,  39.7500000) (  2d30'0.00"W, 39d45'0.00"N)
Band 1 Block=281x14 Type=Int16, ColorInterp=Gray

I tryed with different SRS fields, forcing +north, several epsg codes without success.

comment:4 by warmerdam, 16 years ago

Resolution: invalid
Status: reopenedclosed

Francisco,

Using "-projwin -8.254 43.033 -8.204 43.083" with this VRT file seems to work fine for me. It is important that the corner coordinates given to -projwin be the upper-left and lower-right corners, not the lower-left and upper-right as you have given in this case.

Note that upper/lower are determined by image orientation, not north/southness.

I still see no need for changes in the code here, and if you want a particular window out of rasters with gdal_translate then you will need to know their orientation. Certainly gdal_translate will not flip a raster for you.

I would note that gdalwarp may be more suitable for your purpose.

Note: See TracTickets for help on using tickets.