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 yOrigin = transform pixelWidth = transform pixelHeight = -transform
Just need some concrete steps (with gdal Python functions).