1

I'm using Landsat 8 imagery, and the Indicar Landsat Processing Tools to produce NDVI imagery from the Landsat 8 source. That works fine.

I then have a bounding box expressed as lat/lon coordinates, and I want to crop the NDVI image using gdalwarp:

gdalwarp -s_srs "+proj=utm +zone=32N +datum=WGS84" -t_srs WGS84 -te 8.45586 56.55605 8.62229 56.78750 LC8_197-021_20130502_LGN01_ndvi.tif cropped112.tif

The original (left) vs. result (right) looks like this: NDVI vs. warped image

I'm sure it has to do with the value of t_srs. I tried not to specify neither s_srs or t_srs which fails with: ERROR 1: Attempt to create 0x0 dataset is illegal,sizes must be larger than zero.

Probably because it expects bounding box in something different that lat/lons? Meters?

gdalinfo output for the original NDVI image:

Driver: GTiff/GeoTIFF
Files: LC8_197-021_20130502_LGN01_ndvi.tif
       LC8_197-021_20130502_LGN01_ndvi.tif.aux.xml
Size is 7921, 7711
Coordinate System is:
PROJCS["WGS 84 / UTM zone 32N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32632"]]
Origin = (364785.000000000000000,6311415.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  364785.000, 6311415.000) (  6d46'41.90"E, 56d55'35.84"N)
Lower Left  (  364785.000, 6080085.000) (  6d53'38.18"E, 54d50'58.48"N)
Upper Right (  602415.000, 6311415.000) ( 10d40'59.05"E, 56d56' 6.12"N)
Lower Right (  602415.000, 6080085.000) ( 10d35'43.57"E, 54d51'26.50"N)
Center      (  483600.000, 6195750.000) (  8d44'15.63"E, 55d54'24.87"N)
Band 1 Block=7921x1 Type=Float32, ColorInterp=Gray
  Min=-0.232 Max=0.893 
  Minimum=-0.232, Maximum=0.893, Mean=0.138, StdDev=0.210
  Metadata:
    STATISTICS_MAXIMUM=0.89341741800308
    STATISTICS_MEAN=0.1375360206514
    STATISTICS_MINIMUM=-0.23159252107143
    STATISTICS_STDDEV=0.2095535427118

gdalinfo output of the warped image:

Driver: GTiff/GeoTIFF
Files: cropped112.tif
Size is 424, 589
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (8.455859999999999,56.787500000000001)
Pixel Size = (0.000392523584906,-0.000392954159593)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (   8.4558600,  56.7875000) (  8d27'21.10"E, 56d47'15.00"N)
Lower Left  (   8.4558600,  56.5560500) (  8d27'21.10"E, 56d33'21.78"N)
Upper Right (   8.6222900,  56.7875000) (  8d37'20.24"E, 56d47'15.00"N)
Lower Right (   8.6222900,  56.5560500) (  8d37'20.24"E, 56d33'21.78"N)
Center      (   8.5390750,  56.6717750) (  8d32'20.67"E, 56d40'18.39"N)
Band 1 Block=424x4 Type=Float32, ColorInterp=Gray
3

From http://www.gdal.org/gdalwarp.html

-te xmin ymin xmax ymax: set georeferenced extents of output file to be created (in target SRS).

Without -s_srs and -t_srs the target SRS will be the same as the original SRS of the image which is EPSG:32632 as you can see from the gdalinfo report. You can also see the extent in EPSG:32632

Upper Left  (  364785.000, 6311415.000)
Lower Right (  602415.000, 6080085.000)

Your command without -s_srs and -t_srs asks for an area (8.45586 56.55605 8.62229 56.78750) in native EPSG:32632 which is

  • outside the image data
  • of size 0.17 by 0.23 meters

Small size which is less than one pixel is the reason for the "Attempt to create 0x0 dataset" error

Your latter command works because you reproject your image into EPSG:4326 and then the values of -te makes sense. It is your responsibility to study why reprojecting from 32632 into 4326 makes the image different - warped.

If your aim is to crop and keep the original SRS, use gdal_translate http://www.gdal.org/gdal_translate.html with -projwin

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

  • Works great with gdal_translate. Thanks for clearing up the concepts, makes total sense now. GIS and processing satellite imagery is still quite new to me. – MartinHN Mar 19 '15 at 19:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.