#AED FROM ROAD DENSITY FILE REFORMATTING FROM .TXT DATA
# Different methods for removing urban roads from AED file


cdo -f nc -input,gridfile.txt < grip4_total_dens_m_km2_oneline.asc test.nc 
cdo setmissval,9.96921e+36 test.nc test1.nc
cdo setctomiss,-9999 test1.nc test2.nc
cdo chname,var1,Rd_Density_m_km2 test2.nc test3.nc
#equalising areas grid to 
cdo remapcon,Weighted_PFT_7_8_9_ROS_discount_factor_nonNAN.nc test3.nc test4.nc
cdo mul test4.nc Areas_km2.nc test5.nc
cdo chname,Rd_Density_m_km2,Rd_Density_m_grid test5.nc grip4_total_dens_m_per_gridcell_360x720.nc

# URBAN ROAD SUBTRACTION
cdo mul ssp_2018_regrid.nc Areas_km2.nc temp1.nc
#Currently, assume xx km/km^2 of urban area
cdo mulc,5 temp1.nc temp2.nc
cdo mulc,1000 temp2.nc temp3.nc #convert to metres length 
cdo sub grip4_total_dens_m_per_gridcell_360x720.nc temp3.nc temp4.nc
ncap2 -s 'where (Rd_Density_m_grid<0) Rd_Density_m_grid=0;' temp4.nc grip4_NO_URBAN_m_per_grid_5km.nc
rm -rf temp*

# URBAN ROAD SUBTRACTION
cdo mul ssp_2018_regrid.nc Areas_km2.nc temp1.nc
#Currently, assume xx km/km^2 of urban area
cdo mulc,10 temp1.nc temp2.nc
cdo mulc,1000 temp2.nc temp3.nc #convert to metres length 
cdo sub grip4_total_dens_m_per_gridcell_360x720.nc temp3.nc temp4.nc
ncap2 -s 'where (Rd_Density_m_grid<0) Rd_Density_m_grid=0;' temp4.nc grip4_NO_URBAN_m_per_grid_10km.nc
rm -rf temp*

# URBAN ROAD SUBTRACTION
cdo mul ssp_2018_regrid.nc Areas_km2.nc temp1.nc
#Currently, assume xx km/km^2 of urban area
cdo mulc,15 temp1.nc temp2.nc
cdo mulc,1000 temp2.nc temp3.nc #convert to metres length 
cdo sub grip4_total_dens_m_per_gridcell_360x720.nc temp3.nc temp4.nc
ncap2 -s 'where (Rd_Density_m_grid<0) Rd_Density_m_grid=0;' temp4.nc grip4_NO_URBAN_m_per_grid_15km.nc
rm -rf temp*

ncap2 -s 'where(ssp1_2010<0.1) ssp1_2010 =0;' ssp_2018_regrid.nc temp1.nc
ncap2 -s 'where(ssp1_2010>0) ssp1_2010 =1;' temp1.nc gt0.1mask.nc
rm -rf temp*

ncap2 -s 'where(ssp1_2010<0.01) ssp1_2010=0;' ssp_2018_regrid.nc temp1.nc
ncap2 -s 'where(ssp1_2010>0.1) ssp1_2010=0;' temp1.nc temp2.nc
ncap2 -s 'where(ssp1_2010>0) ssp1_2010 =1;' temp2.nc gt0.01mask.nc
rm -rf temp*


cdo setrtoc,0.01,1,-1 ssp_2018_regrid.nc temp2.nc
ncap2 -s 'where(ssp1_2010=0) ssp1_2010=1;' temp1.nc temp2.nc

ncap2 -s 'where(ssp1_2010>0.01) ssp1_2010=-1;' ssp_2018_regrid.nc temp1.nc
cdo setrtoc,0,0,1 temp1.nc temp2.nc
cdo setrtoc,-1,-0.1,0 temp2.nc temp3.nc
ncap2 -s 'where(ssp1_2010>0) ssp1_2010=1;' temp3.nc lt0.01mask.nc




### BELOW INCLUDES URBAN ROADS
#ones mask
ncap2 -s 'where(Rd_Density_m_grid>1) Rd_Density_m_grid=1;' grip4_total_dens_m_per_gridcell_360x720.nc temp1.nc
cdo setrtoc,0,1,1 temp1.nc ones.nc
rm -rf temp1.nc
#AED from density
cdo div ones.nc grip4_total_dens_m_per_gridcell_360x720.nc temp1.nc
cdo mulc,500000000 temp1.nc AED_360x720.nc
cdo chname,Rd_Density_m_grid,AED AED_360x720.nc AED.nc
rm -rf temp1.nc
ncap2 -s 'where(FRI>0) FRI=0;' Weighted_PFT_7_8_9_ROS_discount_factor_nonNAN.nc temp1.nc
cdo add temp1.nc AED.nc temp2.nc
cdo chname,FRI,AED temp2.nc AED_testnew.nc
ncap2 -s where(AED>100000) AED=0 AED_testnew.nc temp1.nc
cdo setctomiss,0 temp1.nc AED_testnew_100klimit.nc


###

#BELOW EXCLUDES URBAN ROADS
ncap2 -s 'where(Rd_Density_m_grid>1) Rd_Density_m_grid=1;' grip4_NO_URBAN_m_per_grid_4km.nc temp1.nc
cdo setrtoc,0,1,1 temp1.nc ones_NU.nc
rm -rf temp1.nc
#AED from density
cdo div ones_NU.nc grip4_NO_URBAN_m_per_grid_4km.nc temp1.nc
cdo mulc,500000000 temp1.nc AED_NU_360x720.nc
cdo chname,Rd_Density_m_grid,AED AED_NU_360x720.nc AED_NU.nc
rm -rf temp1.nc
ncap2 -s 'where(FRI>0) FRI=0;' Weighted_PFT_7_8_9_ROS_discount_factor_nonNAN.nc temp1.nc
cdo add temp1.nc AED_NU.nc temp2.nc
cdo chname,FRI,AED temp2.nc AED_NU_testnew.nc
rm -rf temp1.nc
ncap2 -s 'where(AED>100000) AED=0;' AED_NU_testnew.nc temp1.nc
cdo setctomiss,0 temp1.nc AED_NO_URBAN_4km.nc
rm -rf *NU* grip4_NO_URBAN_m_per_grid.nc temp*

ncap2 -s 'where(Rd_Density_m_grid>1) Rd_Density_m_grid=1;' grip4_NO_URBAN_m_per_grid_8km.nc temp1.nc
cdo setrtoc,0,1,1 temp1.nc ones_NU.nc
rm -rf temp1.nc
#AED from density
cdo div ones_NU.nc grip4_NO_URBAN_m_per_grid_8km.nc temp1.nc
cdo mulc,500000000 temp1.nc AED_NU_360x720.nc
cdo chname,Rd_Density_m_grid,AED AED_NU_360x720.nc AED_NU.nc
rm -rf temp1.nc
ncap2 -s 'where(FRI>0) FRI=0;' Weighted_PFT_7_8_9_ROS_discount_factor_nonNAN.nc temp1.nc
cdo add temp1.nc AED_NU.nc temp2.nc
cdo chname,FRI,AED temp2.nc AED_NU_testnew.nc
rm -rf temp1.nc
ncap2 -s 'where(AED>100000) AED=0;' AED_NU_testnew.nc temp1.nc
cdo setctomiss,0 temp1.nc AED_NO_URBAN_8km.nc
rm -rf *NU* grip4_NO_URBAN_m_per_grid.nc temp*

ncap2 -s 'where(Rd_Density_m_grid>1) Rd_Density_m_grid=1;' grip4_NO_URBAN_m_per_grid_12km.nc temp1.nc
cdo setrtoc,0,1,1 temp1.nc ones_NU.nc
rm -rf temp1.nc
#AED from density
cdo div ones_NU.nc grip4_NO_URBAN_m_per_grid_12km.nc temp1.nc
cdo mulc,500000000 temp1.nc AED_NU_360x720.nc
cdo chname,Rd_Density_m_grid,AED AED_NU_360x720.nc AED_NU.nc
rm -rf temp1.nc
ncap2 -s 'where(FRI>0) FRI=0;' Weighted_PFT_7_8_9_ROS_discount_factor_nonNAN.nc temp1.nc
cdo add temp1.nc AED_NU.nc temp2.nc
cdo chname,FRI,AED temp2.nc AED_NU_testnew.nc
rm -rf temp1.nc
ncap2 -s 'where(AED>100000) AED=0;' AED_NU_testnew.nc temp1.nc
cdo setctomiss,0 temp1.nc AED_NO_URBAN_12km.nc
rm -rf *NU* grip4_NO_URBAN_m_per_grid.nc temp*

cdo mul AED_NO_URBAN_12km.nc gt0.1mask.nc temp1.nc
cdo mul AED_NO_URBAN_8km.nc gt0.01mask.nc temp2.nc
cdo mul AED_NO_URBAN_4km.nc lt0.01mask.nc temp3.nc
cdo add temp1.nc temp2.nc temp4.nc
cdo add temp4.nc temp3.nc AED_DISCRIM_SUM_5_15km.nc


rm -rf temp* AED_NO_URBAN_* AED_DISCRIM_SUM.nc *mask.nc grip4_NO_URBAN_m_per_grid_5km.nc grip4_NO_URBAN_m_per_grid_5km.nc grip4_NO_URBAN_m_per_grid_10km.nc grip4_NO_URBAN_m_per_grid_15km.nc


############# ABOVE BUT USING LINREG COEFFICIENT

cdo mul ssp_2018_regrid.nc Areas_km2.nc temp1.nc
#Currently, assume xx km/km^2 of urban area
cdo mulc,11.08 temp1.nc temp2.nc
cdo mulc,1000 temp2.nc temp3.nc #convert to metres length 
cdo sub grip4_total_dens_m_per_gridcell_360x720.nc temp3.nc temp4.nc
ncap2 -s 'where (Rd_Density_m_grid<0) Rd_Density_m_grid=0;' temp4.nc grip4_NO_URBAN_m_per_grid_SLOPE.nc
rm -rf temp*

ncap2 -s 'where(ssp1_2010<0.01) ssp1_2010 =0;' ssp_2018_regrid.nc temp1.nc
ncap2 -s 'where(ssp1_2010>0) ssp1_2010 =1;' temp1.nc gt0.01allmask.nc


ncap2 -s 'where(Rd_Density_m_grid>1) Rd_Density_m_grid=1;' grip4_NO_URBAN_m_per_grid_SLOPE.nc temp1.nc
cdo setrtoc,0,1,1 temp1.nc ones_NU.nc
rm -rf temp1.nc
#AED from density
cdo div ones_NU.nc grip4_NO_URBAN_m_per_grid_SLOPE.nc temp1.nc
cdo mulc,500000000 temp1.nc AED_NU_360x720.nc
cdo chname,Rd_Density_m_grid,AED AED_NU_360x720.nc AED_NU.nc
rm -rf temp1.nc
ncap2 -s 'where(FRI>0) FRI=0;' Weighted_PFT_7_8_9_ROS_discount_factor_nonNAN.nc temp1.nc
cdo add temp1.nc AED_NU.nc temp2.nc
cdo chname,FRI,AED temp2.nc AED_NU_testnew.nc
rm -rf temp1.nc
ncap2 -s 'where(AED>100000) AED=0;' AED_NU_testnew.nc temp1.nc
cdo setctomiss,0 temp1.nc AED_NO_URBAN_SLOPE.nc
rm -rf *NU* grip4_NO_URBAN_m_per_grid.nc temp*
cdo mul AED_NO_URBAN_SLOPE.nc gt0.01allmask.nc temp1.nc
rm -rf AED_NO_URBAN_SLOPE.nc
mv temp1.nc AED_NO_URBAN_SLOPE.nc 

