Code for converting LULC image from Bhuvan WMS service into an xarray dataset and saving it as a netcdf file

lulc_conversion.py

# the lulc class codes will be in the 'class_new' band. What those codes represent can be seen in the url below
	# https://bhuvan-vec2.nrsc.gov.in/bhuvan/wms?REQUEST=GetLegendGraphic&VERSION=1.0.0&FORMAT=application/json&WIDTH=20&HEIGHT=20&LAYER=lulc:KL_LULC50K_1112
	# use the 'red_new', 'green_new' and 'blue_new' bands to visualise the data
	# Be aware there will be a bhuvan watermark, usually bottom-right, tried to encode it as white (255,255,255) 
		# and lulc class code 99

from rasterio.io import MemoryFile
import requests 
import json 
import rioxarray
from colormath.color_objects import sRGBColor, LabColor
from colormath.color_conversions import convert_color
from colormath.color_diff import delta_e_cie2000
import numpy as np
from geopy.distance import distance as geopy_dist

def download_lulc_data(state, bbox, period, save_path):

	'''
	This function gets 1:50,000 or 25m resolution land use/land cover imagery \n
	from bhuvan's wms service, estimates the underlying data, and \n
	saves that data as a netcdf file.

	Parameters:
	state -- location of area of interest. Use 2-letter codes like KL for Kerala
	bbox -- bounding box for area of interest: [min_lon,min_lat,max_lon,max_lat] 
	period -- 3 options are '2005-06', '2011-12' and '2015-16'
	save_path -- where netcdf will be saved 
	
	Returns the data as an xarray dataset
	'''

	wms_url = 'https://bhuvan-vec2.nrsc.gov.in/bhuvan/wms'

	# estimating how many 25m pixels high and wide the area is 
	south_west = (bbox[1],bbox[0])
	south_east = (bbox[1],bbox[2])
	north_east = (bbox[3],bbox[2])
	height = round(geopy_dist(south_east, north_east).meters/25)
	width = round(geopy_dist(south_west, south_east).meters/25)
	
	bbox_string = ','.join([str(x) for x in bbox])

	# querying bhuvan wms with specific parameters
	if period == '2005-06':
		layer_name = 'lulc:' + state + '_LULC50K_0506'
		style_name = layer_name
	elif period == '2011-12':
		layer_name = 'lulc:' + state + '_LULC50K_1112'
		style_name = 'lulc:LULC50K_1112_NEW'
	elif period == '2015-16':
		layer_name = 'lulc:' + state + '_LULC50K_1516'
		style_name = 'lulc:LULC50K_1516_NEW'
	
	params_main = {
		'REQUEST': 'GetMap',
		'LAYERS': layer_name,
		'STYLES': style_name,
		'CRS': 'EPSG:4326',
		'BBOX': bbox_string,
		'WIDTH': width, 
		'HEIGHT': height,
		'FORMAT': 'image/geotiff', 
		}
	response_main = requests.get(wms_url, params=params_main)

	# converting response to xarray dataset that's easier to work with 
	img_darray = rioxarray.open_rasterio(MemoryFile(response_main.content)\
				.open(driver='GTiff'))
	img_dset = img_darray.to_dataset(dim='band')
	img_dset = img_dset.rename_vars({1:'red', 2:'green', 3:'blue'})


	# getting legend json that allows us to estimate underlying data
	params_legend = {
	'REQUEST': 'GetLegendGraphic',
	'LAYER': layer_name,
	'STYLE': style_name,
	'FORMAT': 'application/json',
				}
	response_legend = requests.get(wms_url, params=params_legend)
	legend_dict = json.loads(response_legend.content)

	# converting legend json into a simpler dict
	legend_dict_simple = {}
	for rule in legend_dict['Legend'][0]['rules']:
		bhuvan_webcode = int(rule['name'])
		bhuvan_class = rule['title']
		bhuvan_color_hex = rule['symbolizers'][0]['Polygon']['fill']
		bhuvan_color_rgb =\
			tuple(int(bhuvan_color_hex.lstrip('#')[i:i+2], 16) for i in (0, 2, 4))	
		legend_dict_simple[bhuvan_color_rgb] = {} # rgb tuples are dict keys here 
		legend_dict_simple[bhuvan_color_rgb]['bhuvan_webcode'] = bhuvan_webcode
		legend_dict_simple[bhuvan_color_rgb]['bhuvan_class'] = bhuvan_class
		
	# adding explicit 'no data' class to dict, not there in bhuvan classes otherwise
	legend_dict_simple[(255,255,255)] =\
		{'bhuvan_webcode':99, 'bhuvan_class':'No Data'}
		
	# converting colors from rgb colorspace to cie colorspace 
		# conversion necessary to determine closest color in legend json 
	palette_colors_rgb = tuple(legend_dict_simple.keys())
	palette_colors_srgb = [sRGBColor(x[0]/255, x[1]/255, x[2]/255) \
				for x in palette_colors_rgb]
	palette_colors_lab = [convert_color(x, LabColor) for x in palette_colors_srgb]

	# defining function to find closest color in legend json 
	def get_closest_palette_entry(red_val, green_val, blue_val):

		val_tuple = (red_val, green_val, blue_val)
		if val_tuple in palette_colors_rgb:
			closest_palette_rgb_tuple = val_tuple
		# encoding watermark colors as white
		elif ((red_val == green_val == blue_val) and (red_val != 225)):
			closest_palette_rgb_tuple = (255,255,255)
		else:
			rgb_tuple_srgb = sRGBColor(val_tuple[0]/255, val_tuple[1]/255, 
									val_tuple[2]/255)
			rgb_tuple_lab = convert_color(rgb_tuple_srgb, LabColor)
			dist_list =\
				[delta_e_cie2000(rgb_tuple_lab, x) for x in palette_colors_lab]
			min_pos = dist_list.index(min(dist_list))
			closest_palette_rgb_tuple = palette_colors_rgb[min_pos]

		closest_bhuvan_webcode =\
			legend_dict_simple[closest_palette_rgb_tuple]['bhuvan_webcode']

		return closest_palette_rgb_tuple[0], closest_palette_rgb_tuple[1],\
				closest_palette_rgb_tuple[2], closest_bhuvan_webcode

	# applying function to every pixel
	new_ufunc = np.frompyfunc(get_closest_palette_entry,3,4)
	closest_palette_entries = new_ufunc(img_dset['red'], img_dset['green'], 
						img_dset['blue'])

	# putting r,g,b values and estimated landuse/landcover class in dataset 
		# as new variables 
	img_dset['red_new'] = closest_palette_entries[0]
	img_dset['green_new'] = closest_palette_entries[1]
	img_dset['blue_new'] = closest_palette_entries[2]
	img_dset['class_new'] = closest_palette_entries[3]

	img_dset  = img_dset.astype('uint8')
	img_dset.to_netcdf(save_path)
	
	return img_dset

# example of usage, put in your own requirements
state = 'KL'
bbox = (76.5368557, 9.5704093, 76.5883541, 9.628581)
period = '2011-12'
save_path = 'kl_2011_12_lulc.nc'

xarray_dset = download_lulc_data(state, bbox, period, save_path)

# open saved netcdf with xarray_dset = xr.open_dataset('kl_2011_12_lulc.nc', decode_coords='all')