rem ILWIS Script
//The parameter %1 refers to the name of the factor map. It should be less than 7 characters long.
// Make sure that each map has a domain with the same name

//FIRST WE WILL DELETE EXISTING RESULT FILES
// the crosstable s%1.tbt
//The attribute table %1.tbt
// and we make a new attribute table

del s%1.*
del w%1.*
del %1.tbt
crtbl %1 %1

//NOW WE CROSS THE FACTOR MAP WITH THE ACTIVITY MAP
// The landslide map should be called ACTIVE and should have either 0 or 1 values. 1 values mean landslides. 
// The cross table is called s%1

s%1=TableCross(%1.mpr,active.mpr,IgnoreUndefs)
calc s%1.tbt

//Now we calculate one column in the cross table to indicate only the pixels with landslides.

Tabcalc s%1 npixact=iff(active=1,NPix,0)

//NOW WE USE AGGREGATION FUNCTION, WITH OR WITHOUT A KEY TO CALCULATE:
//NCLASS = number of pixels in the class. We sum the values from columns Npix and group them by %1
//nslclass = number of pixels with landslides in the class.We sum the values from columns Npixact and group them by %1
//nmap = number of pixels with landslides in the map. We sum the values from columns Npix and don't group them 
//nslide = number of pixels with landslide in the map. We sum the values from columns Npixact and don't group them 
//THE RESULTS ARE NOT STORED IN THE CROSS TABLE S%1 BUT IN THE ATTRIBUTE TABLE %1

Tabcalc s%1 %1.nclass = ColumnJoinSum(s%1.tbt,Npix,%1,1)
Tabcalc s%1 %1.nslclass = ColumnJoinSum(s%1.tbt,Npixact,%1,1)
Tabcalc s%1 %1.nmap = ColumnJoinSum(s%1.tbt,Npix,,1)
Tabcalc s%1 %1.nslide = ColumnJoinSum(s%1.tbt,Npixact,,1)

//NOW WE CALCULATE THE FOUR VALUES NPIX1 - NPIX4 AS INDICATED IN THE EXERCISE BOOK. THIS IS DONE IN THE ATTRIBUTE TABLE
// We correct for the situation when Npix1 - Npix3 might be 0 pixels, and change it into 1 pixel 

Tabcalc %1 npix1 =IFF((nslclass>0),nslclass,1)
Tabcalc %1 npix2 = IFF((nslide-nslclass)=0,1,nslide-nslclass)
Tabcalc %1 npix3 = IFF((nclass-nslclass)=0,1,nclass-nslclass)
Tabcalc %1 npix4 = nmap-nslide-nclass+nslclass

//NOW WE CALCULATE THE WEIGHTS IN THE ATTRIBUTE TABLE
Tabcalc %1 wplus {dom=value.dom; vr=-10:10:0.00001} = LN((npix1/(npix1+npix2))/(npix3/(npix3+npix4)))
Tabcalc %1 wminus {dom=value.dom; vr=-10:10:0.000001}  = LN((npix2/(npix1+npix2))/(npix4/(npix3+npix4)))

//NOW WE CALCULATE THE CONTRAST FACTOR
Tabcalc %1 Cw = wplus-wminus

//NOW WE CALCULATE THE FINAL WEIGHT
//The final weight is the sum of the positive weight and the negative weights of the other classes
Tabcalc %1 WminSum=aggsum(wminus)
Tabcalc %1 Wmap=wplus+Wminsum-Wminus

//NOW WE MAKE AN ATTRIBUTE MAP OF THE FINAL WEIGHTS
w%1.mpr = MapAttribute(%1,%1.Wmap)
calc w%1.mpr 