2

I've made a python script which tries to reproduce what gdal_translate does (given a bounding box BBOX in EPSG:4326 and an UTM32 .jp2 file. Then, translates the BBOX from EPSG:4326 to EPSG:32632 (UTM32 North) and subsets the result from the geotiff file.

Gdalinfo of the Sentinel 2 UTM32 file:

gdalinfo 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)

Below is the python script:

from osgeo import gdal
import sys
import ogr, osr
import math
import os

utm32_file = "T32UQD_20200421T102021_TCI_60m.jp2"
gdalwarp_crop_file = "gdal_cropped.tiff"

# GIVEN BBOX in EPSG:4326
X_MIN, Y_MIN, X_MAX, Y_MAX = 12.46826, 52.58852, 12.48516, 52.60034

def project(pointX, pointY):
    # Spatial Reference System
    inputEPSG = 4326
    outputEPSG = 32632

    # create a geometry from coordinates
    point = ogr.Geometry(ogr.wkbPoint)
    point.AddPoint(pointX, pointY)

    # create coordinate transformation
    inSpatialRef = osr.SpatialReference()
    inSpatialRef.ImportFromEPSG(inputEPSG)

    outSpatialRef = osr.SpatialReference()
    outSpatialRef.ImportFromEPSG(outputEPSG)

    coordTransform = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)

    # transform point
    point.Transform(coordTransform)

    # print point in EPSG 4326
    return point.GetX(), point.GetY()

xmin, ymin = project(X_MIN, Y_MIN)
xmax, ymax = project(X_MAX, Y_MAX)

# coordinates to get pixel values for (as tuples of points)
points = [(xmin, ymin), (xmax, ymax)]

# open the raster file
ds = gdal.Open(utm32_file)

if ds is None:
    print('Could not open the raster file')
    sys.exit(1)
else:
    print('The raster file was opened satisfactorily')

# get georeference info
transform = ds.GetGeoTransform() # (-2493045.0, 30.0, 0.0, 3310005.0, 0.0, -30.0)
xOrigin = transform[0] # -2493045.0
yOrigin = transform[3] # 3310005.0
pixelWidth = transform[1] # 30.0
pixelHeight = transform[5] # -30.0

band = ds.GetRasterBand(1) # 1-based index

data = band.ReadAsArray()

grid_x = 0
grid_y = 0

# loop through the coordinates
for point in points:
    x = point[0]
    y = point[1]

    xOffset = int((x - xOrigin) / pixelWidth)
    yOffset = int((y - yOrigin) / pixelHeight)

    if grid_x == 0:
        grid_x = xOffset
        grid_y = yOffset
    else:
        grid_x = abs(xOffset - grid_x) + 1
        grid_y = abs(yOffset - grid_y) + 1

gdal_subset = str(points[0][0]) + " " + str(points[1][1]) + " " + str(points[1][0]) + " " + str(points[0][1])
command = "gdal_translate -projwin " + gdal_subset + " " + utm32_file + " " + gdalwarp_crop_file + " >> nohup.out"
os.system(command)

# PRINT GRID SIZES (width x height)
print "python:    ", grid_x, grid_y

ds = gdal.Open(gdalwarp_crop_file)
print "gdal:      ", ds.RasterXSize, ds.RasterYSize

The problem is the affine transformation in the script sometimes returns correct grid sizes (width x height) as gdal_translate, but in many cases it returns more few pixels than gdal_translate. For example output of the python script is:

python:     19 x 24
gdal:       18 x 23

What could be the problem here to return the same result (grid sizes) as gdal_translate?

1

Your original question was how to calculate the same number of pixels as gdal_translate. I honestly do not know the answer to that question other than by looking at the gdal_translate code, where there are several nuances in cropping the source raster. You could compare the gdal_translate code with your code and possibly step through both codes side by side in a debugger.

I guess my approach to solve this problem would not be to attempt to rewrite gdal_translate, but rather to use the gdal.Translate python method as described in the answer to this question.

1
  • Thanks for correcting my question. Indeed, I will need to copy the way which gdal_translate in C++ does for -projwin. Because, Python code is just a temporary step for me to check the results before applying it to a complex java application which does the same thing. In this application which does not support GDAL version 2.x. May 10 '20 at 7:22

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.