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