From below script:
ds = gdal.Open('../SPI/mar/orig_geotiff/mar_cli_chirps_precip_month1_198101.tif') # Data location
a = ds.ReadAsArray()
nlat,nlon = np.shape(a)
b = ds.GetGeoTransform() #bbox, interval
lon = np.arange(nlon)*b[1]+b[0]
lat = np.arange(nlat)*b[5]+b[3]
basedate = dt.datetime(1980,1,1,0,0,0)
The purpose of these lines is to create the arrays of longitude (lon
) and latitude (lat
) coordinates for the dataset.
In a GeoTIFF file, the geographic information (like the extent, the coordinate reference system, and the resolution) is stored in the metadata of the file. This is accessed using the GetGeoTransform()
method, which returns a tuple of six coefficients that define the affine transformation from the pixel space (row, column) to the geospatial coordinate space (x, y).
The returned tuple is structured like so:
(b[0], b[1], b[2],
b[3], b[4], b[5])
Where:
b[0]
= top left x coordinate (longitude of the top-left corner of the image)b[1]
= w-e pixel resolution (width of one pixel in geographic units - often longitude)b[2]
= rotation, 0 if image is "north up"b[3]
= top left y coordinate (latitude of the top-left corner of the image)b[4]
= rotation, 0 if image is "north up"b[5]
= n-s pixel resolution (height of one pixel in geographic units - often latitude; this is usually negative)
In the script you provided, b = ds.GetGeoTransform()
is retrieving this transformation data, then the following lines:
lon = np.arange(nlon)*b[1]+b[0]
lat = np.arange(nlat)*b[5]+b[3]
are creating arrays of the longitude (lon
) and latitude (lat
) coordinates for each pixel in the dataset. Here nlon
and nlat
are the dimensions of the dataset in pixels.
The np.arange(nlon)
creates an array [0, 1, 2, ..., nlon-1] and np.arange(nlat)
creates an array [0, 1, 2, ..., nlat-1].
Then for each element i
in the np.arange(nlon)
array, the longitude is calculated as b[0] + i*b[1]
and for each element j
in the np.arange(nlat)
array, the latitude is calculated as b[3] + j*b[5]
. This gives us the coordinates of the center of each pixel.