2

I've a Sentinel 2 jp2 in UTM32 with this metadata:

Files: T32UQD_20200421T102021_TCI_60m.jp2
Size is 1830, 1830
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,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AXIS["Latitude",NORTH],
        AXIS["Longitude",EAST],
        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"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32632"]]
Origin = (699960.000000000000000,5900040.000000000000000)
Pixel Size = (60.000000000000000,-60.000000000000000)
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  699960.000, 5900040.000) ( 11d59'40.71"E, 53d12'43.17"N)
Lower Left  (  699960.000, 5790240.000) ( 11d55'40.15"E, 52d13'34.41"N)
Upper Right (  809760.000, 5900040.000) ( 13d38' 3.06"E, 53d 9'33.92"N)
Lower Right (  809760.000, 5790240.000) ( 13d31'51.72"E, 52d10'31.74"N)
Center      (  754860.000, 5845140.000) ( 12d46'18.91"E, 52d41'45.89"N)

I tried to warp this full jp2 by gdalwarp from UTM32 to EPSG:4326

gdalwarp T32UQD_20200421T102021_TCI_60m.jp2 -t_srs EPSG:4326 full_warp.tif

I also make a small subset from the original jp2 file with gdal_translate:

gdal_translate -projwin 798326.68347 5795956.80269 799609.577326 5794397.95642 T32UQD_20200421T102021_TCI_60m.jp2 crop.tif

And I project this small subset from UTM32 to EPSG:4326

 gdalwarp crop.tif -t_srs EPSG:4326 small_warp.tif -dstnodata 0

Then I overlay the full_warp.tif and small_warp.tif by QGIS 2.14.9. However, I realized that the small_warp.tif is not aligned on the full_warp.tif.

On the left hand side is the full_warp.tif and on the right hand side is small_warp.tif overlays on full_warp.tif.

enter image description here

I cannot find out the reason here? Can someone give an idea?

Update 1:

The gdalinfo of the crop.tif by gdal_translate show it is aligned (I don't see the geo bounds and pixel resolutions with odd numbers)

gdalinfo crop.tif 
Driver: GTiff/GeoTIFF
Files: crop.tif
Size is 21, 26
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 = (798300.000000000000000,5796000.000000000000000)
Pixel Size = (60.000000000000000,-60.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  798300.000, 5796000.000) ( 13d22' 8.43"E, 52d14' 0.35"N)
Lower Left  (  798300.000, 5794440.000) ( 13d22' 3.48"E, 52d13'10.01"N)
Upper Right (  799560.000, 5796000.000) ( 13d23'14.66"E, 52d13'57.89"N)
Lower Right (  799560.000, 5794440.000) ( 13d23' 9.69"E, 52d13' 7.55"N)
Center      (  798930.000, 5795220.000) ( 13d22'39.06"E, 52d13'33.95"N)

Update 2: I tried to use gdalwarp with -tr (values come from full_warp.tif) by Python Gdal: dataset.GetGeoTransform() and -tap options

gdalwarp crop.tif -t_srs EPSG:4326 small_warp.tif -dstnodata 0 -tap -tr 0.0007161230037719291 -0.0007161230037719291

The result is still misaligned by row and column as below when overlaying small_warp.tif on full_warp.tif

enter image description here

10
  • Did you make sure that the edges of your crop aligned exactly on the pixel boundary? – Ian Turton May 12 '20 at 14:20
  • @IanTurton I think this is done by gdal_translate under neath? It will align my given bounding box in UTM32 to what it wants to calculate. I never thought gdal user needs to align the bbox. – Bằng Rikimaru May 12 '20 at 14:25
  • I'm pretty sure the pixels start at the edge of the requested bbox and continue across from there - that is how rasters work – Ian Turton May 12 '20 at 14:31
  • @IanTurton: I checked the result of gdal_translate for the given bbox (crop.tif) by gdalinfo. It is aligned with the original JP2 file. Coordinates in UTM 32. Clipped file: Upper Left ( 798300.000, 5796000.000). Original file: Upper Left ( 699960.000, 5900040.000) – Bằng Rikimaru May 12 '20 at 14:40
  • 2
    I recommend also to set the output resolution with the -tr switch. Otherwise the outputs from gdalwarp may lead to slightly different pixel sizes. – user30184 May 12 '20 at 15:38
1

I wanted to post the suggestions from Even Rouault from gdal-dev mailing list:

When you don't specify any target extent or resolution, gdalwarp guesses them from the input raster georeferencing and target SRS. It doesn't try to align on any particular resolution

You may want to use the -tr and -tap switches to specify a target resolution, and constraint the extents to be aligned on it.

Another reason for the misalignment is that you use the default nearest resampling. Using -r cubic could likely improve things

To get perfectly overlapping images, you may also need to specify -et 0

From his idea (same as https://gis.stackexchange.com/users/30184/user30184 with -tap option), I tried these commands and results are aligned on Qgis

gdalwarp -tr 0.0007161230037719291 -0.0007161230037719291 T32UQD_20200421T102021_TCI_60m.jp2 full_warp.tif -t_srs EPSG:4326 -dstnodata 0 -tap

gdalwarp -tr 0.0007161230037719291 -0.0007161230037719291 crop.tif small_warp.tif -t_srs EPSG:4326 -dstnodata 0 -tap

Illustration: https://i.imgur.com/1TwmVMN.png

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.