Skip to content

Instantly share code, notes, and snippets.

@bennyistanto
Created July 13, 2023 00:20
Show Gist options
  • Save bennyistanto/7762c768faab2571ab1c781324059765 to your computer and use it in GitHub Desktop.
Save bennyistanto/7762c768faab2571ab1c781324059765 to your computer and use it in GitHub Desktop.
Creating arrays from GeoTIFFs

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment