Recently, I have been working on algorithms which need 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 resolutions (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 raster 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 converts them into a 2d index.

Hope this helps :)