Programme

Bonjour,

Voila, je bataille sur un problème depuis deux semaines, j’ai donc décidé de venir demander votre aide :

J’ai créer ce programme afin de calculer des statistiques zonales (moyennes, etc…) sur des bassins versants (polygones).
Ce programme (ci-dessous) marche très bien lorsque je le croise avec un raster au km, mais des lors que je veux le faire fonctionner avec un raster au 100m (malgré que j’apporte les modifications adéquates), les deux dernières lignes du programme ne s’effectue pas.

Je vais vous mettre ce programme, et les différents programmes dont il se sert pour certaines opérations. J’espère avoir été a peu prés clair dans mon explication

Merci à tous

1 : Programme de Stat Zonales :

#source("C:/Documents and Settings/nom/Bureau/mnt50m_ingnpaca/TRAITEMENTS SIG.r")

# ZONE D'ETUDE : créaction d'un raster "unitaire" équivalent à la grille de travail (domaine d'étude : km ROND)
nbX=1856 ; nbY=1805 ; xll=845.9 ; yll=1841.9 ; dX=0.1 ; dY=0.1
domaine=list(nbX=nbX,nbY=nbY,xll=xll,yll=yll,dX=dX,dY=dY,val=array(1,dim=c(nbX,nbY)))
memory.limit(4000) ; memory.size(max=T)

#chemin entier fichiershape contenant bv
f.shp="C:/Documents and Settings/nom/Bureau/mnt50m_ingnpaca/BV040506/bvtest2.shp"

#chemin sortie bv, creation grille sur les BV
x=fc.XYbassin(f.shp,rep.sortie="C:/Documents and Settings/nom/Bureau/mnt50m_ingnpaca/bv456/",graph=F,ID=1,IDsurf=2,pas=100)

#grille raster(ascii)
grille=fc.read.raster("C:/Documents and Settings/nom/Bureau/mnt50m_ingnpaca/drain456bis.asc")
grille=fc.raster.rond(grille)

#codage des BV
CODBV2=fc.CODBV.raster(domaine,x)

#Stat zonales
y=fc.extract.raster(domaine,grille) #; fc.plot.raster(x,"PJ10 SHYREG")

y=fc.statzonale(y$val,CODBV2,f.sortie="C:/Documents and Settings/nom/Bureau/mnt50m_ingnpaca/PATATE.txt",lacune=T)

2. Le texte souligné en rouge, sont les modifications apportés pour un raster au 100 m , voir ci-dessous sur le programme qui est utilisé et que je charge a la première ligne :

fc.read.raster = function (f.raster,graph=F) {
info=scan(f.raster,list("",""),nmax=6)
a=as.numeric(info2[1:5])
nbX=a[1] ; nbY=a[2] ; xll=a[3]/1000 ; yll=a[4]/1000 ; dX=a[5]/1000 ; dY=a[5]/1000
z=array(scan(f.raster,skip=6,na.strings=info2[6]),dim=c(nbX,nbY)) ; z=z[,nbY:1]
z[which(z==as.numeric(info2[6]),arr.ind=T)]=NA # par exemple sur NODATA=-99.00 et -99 dans les données
raster=list(nbX=nbX,nbY=nbY,xll=xll,yll=yll,dX=dX,dY=dY,val=z)
if (graph) {
x=seq(xll+dX/2,xll+dX/2+(nbX-1)*dX,dX)
y=seq(yll+dY/2,yll+dX/2+(nbY-1)*dY,dY)
windows() ; filled.contour(x,y,z,color = terrain.colors,asp=1,sub=f.raster)
}

Et voici le programme de Stat zonales qui échoue à la fin :

EXTRACTION des valeurs de raster2 sur le domaine de raster1
fc.extract.raster=function(raster1,raster2) {
# Extraction sur la grille de travail
if (raster1$dX ==1 & raster2$dX !=1) raster2=fc.resolution.1x1(raster2)
if (raster2$dX !=1 & raster1$dX!=raster2$dX) break
xa=seq(raster1$xll+raster1$dX/2,by=raster1$dX,length=raster1$nbX)
ya=seq(raster1$yll+raster1$dY/2,by=raster1$dY,length=raster1$nbY) # raster d'extraction (km ronds)
xb=seq(raster2$xll+raster2$dX/2,by=raster2$dX,length=raster2$nbX)
yb=seq(raster2$yll+raster2$dY/2,by=raster2$dY,length=raster2$nbY) # raster à extraire (km ronds)
xbis=seq(raster2$xll+raster2$dX/2,by=raster2$dX,length=raster2$nbX)
ybis=seq(raster2$yll+raster2$dY/2,by=raster2$dY,length=raster2$nbY) # raster à extraire (km ronds)
xc=intersect(trunc(xa),trunc(xb)) ; yc=intersect(trunc(ya),trunc(yb)) # zone commune
decal.xa=xa[1]-trunc(xa[1]) ; decal.xb=xb[1]-trunc(xb[1])
decal.ya=ya[1]-trunc(ya[1]) ; decal.yb=yb[1]-trunc(yb[1])
raster1$val=raster1$val*NA
raster1$val[(xc-raster1$xll-raster1$dX/2)/raster1$dX+1+decal.xa,(yc-raster1$yll-raster1$dY/2)/raster1$dY+1+decal.ya]=raster2$val[(xc-raster2$xll-raster2$dX/2)/raster2$dX+1+decal.xb,(yc-raster2$yll-raster2$dY/2)/raster2$dY+1+decal.yb]
raster1
}