1

I have a GeoTIFF file which is in UTM 32651 and I need to develop a program, using GDAL library version 1.x (e.g., not using gdalwarp/gdal_translate binary tools) which takes a subset bounding box in EPSG:4326 and applies this subset on the GeoTIFF file (i.e: crop the GeoTIFF file), after that it projects the cropped file from UTM 32651 to EPSG:4326.

For example, I have this input GeoTIFF file in UTM 32651

Size is 9371, 12381
Coordinate System is:
PROJCS["WGS 84 / UTM zone 51N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32651"]]
Origin = (200085.000000000000000,2823015.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Upper Left  (  200085.000, 2823015.000) (120d 1' 0.12"E, 25d29'38.49"N)
Lower Left  (  200085.000, 2451585.000) (120d 5'32.63"E, 22d 8'36.74"N)
Upper Right (  481215.000, 2823015.000) (122d48'46.95"E, 25d31'27.38"N)
Lower Right (  481215.000, 2451585.000) (122d49' 4.06"E, 22d10' 9.72"N)
Center      (  340650.000, 2637300.000) (121d26' 7.36"E, 23d50'21.01"N)

And I have a requesting bounding box in EPSG:4326 like this:

Lat_Min=23.5546875   Long_Min=120.849609375
Lat_Max=23.642578125 Long_Max=120.9375

I knew that I could translate a pair of Long, Lat (EPSG:4326) to X Y (EPSG:32651) e.g:

gdaltransform -s_srs EPSG:4326 -t_srs EPSG:32651

However, I feel like I'm missing something about affine transformation here as the bounding box is a square and I cannot just project each pair of Long, Lat in EPSG:4326 to EPSG:32651 and crop the input GeoTIFF image by these transformed coordinates and after that warp the cropped GeoTIFF from UTM EPSG:32651 to EPSG:4326.

So my question is, how can I use GDAL (e.g: gdal python library) to work on this case properly? I knew so far about opening GeoTIFF file with GDAL Python, e.g:

from osgeo import gdal

driver = gdal.GetDriverByName('GTiff')
filename = "rescale.tiff" #path to raster in UTM 32651
dataset = gdal.Open(filename)
band = dataset.GetRasterBand(1)

cols = dataset.RasterXSize
rows = dataset.RasterYSize

transform = dataset.GetGeoTransform()

xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = -transform[5]

Just need some concrete steps (with gdal Python functions).

Vince
  • 20,017
  • 15
  • 45
  • 64
Bằng Rikimaru
  • 483
  • 1
  • 3
  • 13
  • @user30184 thanks for quick reply, the program I need to develop can use GDAL version 1.11.x only (not 2.1 to call gdal_translate as in the answer you showed). – Bằng Rikimaru Feb 16 '19 at 18:08
  • 1
    If you could use a current GDAL version then this answer would suit you https://gis.stackexchange.com/questions/42584/how-to-call-gdal-translate-from-python-code/237411#237411. You should really consider to upgrade because GDAL 1.11 has been end-of-life for quite a long time already but I understand the problems with developing on top of old existing systems. – user30184 Feb 16 '19 at 18:10
  • 1
    Hi, I am not familiar with what functions are available in gdal 1.11.x, is the gdal.ReprojectImage available? If so I can potentially help. – Jascha Muller Mar 12 '19 at 10:30

0 Answers0