#Now we look at spatial trends
 parana.trendest<-lm(parana$data~parana$coords) #linear trend
 summary(parana.trendest)

#Take out the linear trend and look at the residuals
 parana.detrend <- as.geodata(cbind(parana$coords,parana.trendest$residuals))
 plot.geodata(parana.detrend)

 parana2 <- parana$coords^2
 parana.trendest2<-lm(parana$data~parana$coords*parana2) #quadratic trend
 parana.detrend2 <- as.geodata(cbind(parana$coords,parana.trendest2$residuals))
 plot.geodata(parana.detrend2)

#look at the variogram of the residuals from the quadratic fit
 parana.cloud2 <- variog(parana.detrend2,option="cloud")
 plot(parana.cloud2)
 parana.bin2 <- variog(parana.detrend2,uvec = max(parana.cloud2$u)*seq(0,1,l=11))
 plot(parana.bin2)
 env.mc2 <- variog.mc.env(parana.detrend2,obj.var=parana.bin2)
 plot(parana.bin2,envelope = env.mc2)