Using Geotiff as datasource via gdal

Recently, I have been working on algorithms which needs elevation data as well as Land Cover data, with world coverage. Google has an excellent elevation API however free usage comes with a limit.

While searching, I came across a dataset in geotiff format for landcover as well as a processed version of world elevation.  Elevation data comes in various resolution (250m,500m,1km), landcover is 500m .

So how do we read it ? In python its quite easy to use osgo/gdal library and read all bands. Geotiff is a rster file and values can be packed in every band (which is basically a 2d array).

landUsageLookup = [
    "Water",  # 1
    "Evergreen Needle leaf Forest",  # 2
    "Evergreen Broadleaf Forest",  # 3
    "Deciduous Needle leaf Forest",  # 4
    "Deciduous Broadleaf Forest",  # 5
    "Mixed Forests",  # 6
    "Closed Shrublands",  # 7
    "Open Shrublands",  # 8
    "Woody Savannas",  # 9
    "Savannas",  # 10
    "Grasslands",  # 11
    "Permanent Wetland",  # 12
    "Croplands",  # 13
    "Urban and Built-Up",  # 14
    "Cropland/Natural Vegetation Mosaic",  # 15
    "Snow and Ice",  # 16
    "Barren or Sparsely Vegetated"  # 17
]

def getLandUsage(latitude, longitude):
    ds = gdal.Open("data/LCType.tif") 
    gt = ds.GetGeoTransform()
     
    band = ds.GetRasterBand(1)
    px = int((longitude - gt[0]) / gt[1])  # x pixel
    py = int((latitude - gt[3]) / gt[5])  # y pixel
    structval = band.ReadRaster(px, py, 1, 1, buf_type=gdal.GDT_UInt16)
    lc = struct.unpack('h', structval)[0]
    return lc, landUsageLookup[lc]

The code sample above takes in Lat/Lng in wgs84  and convert them into 2d index.

Hope this helps :)

Varun Pant

Using Geotiff as datasource via gdal   by  Varun Pant