[R-sig-Geo] Probability map of presence(1)-absence(0) data including trend and covariate

[TM]

Hi,
I want to create a biodiversity species (probability) map of Spain from 
presence(1)-absence(0) data (approx. 100 sites, and 3 species) and also I 
have data from 400 metereological stations in Spain (Coordinates, 
temperature, rain, hight, etc). I have been looking at using indicator 
kriging on each species layer.

Can I transform this data somehow to still be able to krige? (i.e can I use 
Universal Kriging on binary data?)

If some form of kriging is appropriate, I presume I should model the 
variogram to the same (overall) distance for each species. Is it ok to have 
a different lag distance/number of lags for each species or all species? Are 
there general rules for specifying lag distances?

I also have some metereological data which I could use as a covariate. Is 
Cokriging appropriate to use when one dataset is binary and the other is 
continuous? Does the data need to be co-located?

Are there better statistical methods (in R or other software) to use when 
interpolating presence-absence data, including with a trend and covariate? 
Any references or code in R etc would be much appreciated.

At last, How can I add time dimension to probability map, because I have 
also temperature and rain time series?

Apologies for quite so many questions. Thanks very much

Toni
Departament of Statistics
University of Barcelona