process.py

import rasterio as rio
from rasterio import features
from affine import Affine
import pandas as pd
import numpy as np
import geopandas as gpd
import numexpr as ne

# Input Albers tile
x = 20.0
y = -30.0
pix_res = 25.0
epsg = 3577


# 1- LC proc
ds = rio.open(f"biodiv_workflow/ga_ls_landcover_class_cyear_2_1-0-0_au_x{x}y{y}_2015-01-01_level4.tif")
lc = ds.read(1)


# 2- VE proc
gdf = gpd.GeoDataFrame.from_file("biodiv_workflow/QSC_Extracted_Data_20220307_151321841000-69916/data.gdb/")
gdf = gdf.to_crs(epsg)

trans = Affine(pix_res, 0.0, x*100000, 0.0, -1*pix_res, (y+1)*100000)
plant = features.rasterize(gdf[gdf['COVER']=='plantation'].geometry.values, out_shape=(4000,4000), fill=0, transform=trans, default_value=1, dtype='uint8')
nremn = features.rasterize(gdf[gdf['COVER']=='non-remnant'].geometry.values, out_shape=(4000,4000), fill=0, transform=trans, default_value=1, dtype='uint8')
ocean = features.rasterize(gdf[gdf['COVER']=='ocean'].geometry.values, out_shape=(4000,4000), fill=0, transform=trans, default_value=1, dtype='uint8')

ve = np.ones((4000,4000), dtype='uint8')
ve[plant==1] = 25
ve[nremn==1] = 25
ve[ocean==1] = 50


# 3- LU proc
df = pd.read_csv("biodiv_workflow/FBA_LU_reclass.csv", index_col=0)

gdf = gpd.GeoDataFrame.from_file("biodiv_workflow/QLD_LANDUSE_June_2019/QLD_LANDUSE_June_2019.gdb/")
gdf = gdf.to_crs(epsg)

lut_values = df.index.unique().values

vals = []
for idx in gdf.index:
    val = gdf.iloc[idx]['QLUMP_Code']
    if val in lut_values:
        val = df.loc[[gdf.iloc[idx]['QLUMP_Code']]]['out_code'].values[0]
    else:
        val = -1
    vals.append(val)

gdf['LUT'] = vals

vals = gdf['LUT'].unique()
vals = vals[vals>0]

lu = np.zeros((4000,4000), dtype='uint16')

for val in vals:
    lu += features.rasterize(gdf[gdf['LUT']==val].geometry.values, out_shape=(4000,4000), fill=0, transform=trans, default_value=val, dtype='uint16')


# 4- Derive final product

#rem veg on grazing and other landuses except nature conservation, or native forestry
s1 = ne.evaluate('1*(((lu!=1100)&(lu!=1200)&(lu!=1200)&(lu!=2200)&(lu!=3600))&((ve!=25)&((ve!=99)&(ve!=0))))')
#grazing native pastures
s2 = ne.evaluate('2*((lu!=2100)&(((lc>3)&(lc<8))|((lc>16)&(lc<21))|(lc==84)|(lc==41)|((lc>44)&(lc<55)))&((ve==25)|(ve==99)|(ve==0)))')
s3 = ne.evaluate('3*((lu==2100)&(((lc>7)&(lc<12))|((lc>20)&(lc<25))|(lc==42)|((lc>56)&(lc<70)))&((ve==25)|(ve==99)|(ve==0)))')
s4 = ne.evaluate('4*((lu==2100)&(((lc>11)&(lc<14))|((lc>24)&(lc<27))|((lc>71)&(lc<78)))&((ve==25)|(ve==99)|(ve==0)))')
#native forestry
s22 = ne.evaluate('22*(((lu==1200)|(lu==2200)|(lu==3600))&((ve!=25)&(ve!=99)&(ve!=0)))')
...

plt.imshow(s22)