[Gfoss] v.to.rast
stefano romanelli
romanelli a lamma.rete.toscana.it
Mer 7 Gen 2009 11:44:17 CET
Riesumo questo vecchio post (28 novembre 2008)
> sto utilizzando il comando v.to.rast
> sto convertendo un vettoriale in un raster di 1 km di pixel...
> il raster ha valori 0 e 1 (presenza e assenza bosco).
> come posso fare ad ottenere il pixel soltanto quando il vettoriale
> copre una determinata percentuale?
per:
A) fornire una possibile soluzione al problema: si tratta di pensare al
raster come un grigliato di poligoni vettoriali di lato 1 Km (v.mkgrid)
e di fare una overlay (operator = or). In questo modo dopo opportuni
passaggi è possibile stabilire per ogni quadrato la quantità di area
bosco/non bosco e rasterizzare di conseguenza. Sotto il dettaglio della
metodologia
B) essendo un neofita di grass ho alcune difficoltà concettuali per
quanto riguarda sopratutto il possibile legame del vettoriale a più
tabelle (layer da 1 a x), pertanto vorrei sapere se le soluzioni
adottate nella procedura di sotto sono corrette o se ci sono passaggi
più rapidi
C) ritengo di avere un problema di visualizzazione usando d.mon
sovrapponendo il raster col vettoriale: in pratica se rasterizzo il
grigliato fatto con v.mkgrid, il raster ottenuto e il vettoriale di
origine visualizzato come boundary non si sovrappongono perfettamente;
eppure per come è fatto il v.mkgrid e per quanto riportato nel manuale
sul raster "These values describe the lines which bound the map at its
edges. These lines do NOT pass through the center of the grid cells at
the edge of the map, but along the edge of the map itself." i due
dovrebbero sovrapporsi perfettamente. A riprova se rivettorializzo il
raster i due vettoriali si sovrappongono perfettamente (vedi #NOTA BENE#
sotto). Dove sbaglio?
################################## PROIEZIONE #########################
projection: 1 (UTM)
zone: 32
datum: eur50
ellipsoid: international
north: 4925313
south: 4678313
west: 554799
east: 771799
nsres: 1000
ewres: 1000
rows: 247
cols: 217
################################### CORINE #############################
#scarico il corine della toscana da:
#http://www.centrointerregionale-gis.it/script/corinedownload.asp
#la proiezione è (vedi sopra)
#importo il corine solo parte boschi (terzo livello tutti i 3xx)
v.in.ogr -o dsn=programs/shape/corine/ output=cor_tos layer=toscana
where='code > 299 and code < 400'
#dissolvo...
v.extract -d --o list=1-4756 input=cor_tos output=cor_ext type=area
new=3
#elimino la categoria
v.category input=cor_ext output=cor_del option=del
#riassegno la categoria per la tabella
v.category input=cor_del output=cor_cat type=centroid option=add step=1
#creo la tabella
echo "CREATE TABLE cor_cat (cat integer, tipol integer)" | db.execute
#connetto la tabella alla mappa, ma siamo sicuri che tutti questi
#passaggi siano necessari o posso aggirare in altro modo?
v.db.connect map=cor_cat table=cor_cat
#metto la cat
v.to.db map=cor_cat type=centroid option=cat
#aggiorno il campo tipol col valore del bosco
echo "UPDATE cor_cat SET tipol=3" | db.execute
################################ GRIGLIA ###############################
#creo il grigliato
v.mkgrid map=msk_tos grid=247,217 position=region
#aggiungo la colonna di univoci
echo "ALTER TABLE msk_tos ADD COLUMN id integer" | db.execute
#aggiorno il campo id col valore di cat
echo "UPDATE msk_tos SET id=cat" | db.execute
################################### OVERLAY ############################
#unisco i due layer
v.overlay ainput=cor_cat binput=msk_tos output=cor_msk operator=or
#necessità di fare un ulteriore dissolve sennò sballo il calcolo delle
#aree infatti a causa dei bordi frastagliati e dell'area +tosto grossa
#del pixel il bosco entra ed esce più volte nello stesso pixel come da
#esempio
db.select table=cor_msk sql='SELECT * FROM cor_msk WHERE a_tipol=3
ORDER BY b_cat'
.................................
.................................
1175|386|3|31134|144|103|31134
37318|386|3|31134|144|103|31134
82179|386|3|31134|144|103|31134
37061|386|3|31135|144|104|31135
37320|386|3|31135|144|104|31135
81994|386|3|31135|144|104|31135
82181|386|3|31135|144|104|31135
37062|386|3|31136|144|105|31136
81997|386|3|31136|144|105|31136
.................................
.................................
#fine esempio
#aggiungo una colonna dove sommo 10.000.000 a b_id dove a_tipol = 3
# e 20.000.000 a b_id dove id_tipol IS NULL
#in questo modo tengo traccia del bosco/non bosco (a_id_tipol=3/IS NULL)
#e mantengo anche l'informazione sull'id univoco di ciacuna cella di
#msk_tos (b_id)
echo "ALTER TABLE cor_msk ADD COLUMN mill integer" | db.execute
echo "UPDATE cor_msk SET mill=b_id+10000000 WHERE a_tipol=3" |
db.execute
echo "UPDATE cor_msk SET mill=b_id+20000000 WHERE a_tipol IS NULL" |
db.execute
#ora posso dissolvere
v.dissolve input=cor_msk output=cmk_dis column=mill
#ricreo la tabella di cmk_dis
echo "CREATE TABLE cmk_dis (cat integer, id integer, area double)" |
db.execute
#connetto la tabella alla mappa
v.db.connect map=cmk_dis table=cmk_dis
#rimetto la cat
v.to.db map=cmk_dis type=centroid option=cat
#riottengo l'id univoco levando o 10 o 20 milioni
echo "UPDATE cmk_dis SET id=cat-10000000 WHERE cat<20000000" |
db.execute
echo "UPDATE cmk_dis SET id=cat-20000000 WHERE cat>=20000000" |
db.execute
#controllo (vedi esempio precedente)
db.select table=cmk_dis sql="select * from cmk_dis where cat=10031135"
--cat |id |area
--10031135|31135|
#ok ora c'è solo una riga con id 31135 e tipologia 3=bosco (che
#rimetto in seguito)
db.select table=cmk_dis sql="select * from cmk_dis where id=31135"
cat |id |area
10031135|31135|
20031135|31135|
#in questo caso invece vedo che nel grigliato 31135 c'è sia
#bosco che non bosco
#fine controllo
#ora posso calcolare le aree in Km così è ho già il risultato tra 0 e 1
v.to.db map=cmk_dis option=area units=k column=area
#rimetto anche la tipologia bosco non bosco anche se potevo tenere la
#cat > o < 20 milioni
echo "ALTER TABLE cmk_dis ADD COLUMN tip integer" | db.execute
echo "UPDATE cmk_dis SET tip=3 WHERE cat<20000000" | db.execute
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
#il VALORE AGGIUNTO di tale metodologia è che posso riclassficare i
#bordi delle aree a piacimento e non solo con una sola soglia ad esempio
#> 30% a seconda dell'uso e delle finalità che abbiamo. Cioè tenderei a
#dare un valore diverso al pixel in funzione della percentuale di
#superficie coperta da bosco e non solo 0-1
#Una delle cose principali da tenere in considerazione è la SUPERFICIE
#SPECIFICA: un'area di 10 x 10 pixel consta di 100 pixel di cui 36 di
#bordo che è cosa però ben diversa da 4 aree di 5 x 5 pixel che constano
#in totale di 100 pixel ma 16 x 4 = 64 pixel di bordo
#Se ad esempio devo incrociare la maschera bosco non bosco con una time
#series di immagini NDVI posso avere delle statistiche completamente
#sballate se ho aree boscate fortemente frastagliate e limitate di
#superficie (cosa abbastanza veritiera per i boschi in Italia)
#Un punto di partenza per valutare le diverse soglie può essere quella
#di aggiungere un campo alla tabella e aggiornarlo con v.to.db
#option=compact che valuta la compattezza dell'area, ma questo
#attualmente esula da questa prova
#Nel caso della Toscana i pixel di bordo non sono affatto trascurabili.
#Nella ripartizione in 3 categorie come nel proseguio sottostante il
#numero di pixel con copertura del bosco da 30 all'80% sono circa pari
#al numero di pixel con copertura bosco da 80 a 100% (vedi RASTER MAP
#CATEGORY REPORT più sotto)
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
#aggiungo il campo recl per mettere i valori che assegnerò ai pixel
#bosco della maschera
echo "ALTER TABLE cmk_dis ADD COLUMN recl integer" | db.execute
echo "UPDATE cmk_dis SET recl=1 WHERE (tip=3 AND area >= 0.8)" |
db.execute
echo "UPDATE cmk_dis SET recl=2 WHERE (tip=3 AND area >= 0.5 AND area <
0.8)" | db.execute
echo "UPDATE cmk_dis SET recl=3 WHERE (tip=3 AND area >= 0.3 AND area <
0.5)" | db.execute
#################### PREPARAZIONE PER IL RASTER ########################
#creo una tabella con db.copy mettendo solo i valori dove è presente il
#bosco in questo modo ho solo valori univoci di cat
db.copy from_table=cmk_dis to_table=msk_rec where='recl >= 1'
#connetto la tabella esportata alla msk_tos
v.db.connect map=msk_tos table=msk_rec key=id layer=2
#elimino le categorie per riassegnarla al layer 2
#ma è la maniera giusta di fare?
v.category input=msk_tos output=msk_del option=del
#creo la MAPPA VETTORIALE FINALE per trasformarla poi in raster
v.category input=msk_del output=msk_fin option=add type=centroid layer=2
step=1
#aggiungo la colonna per i colori
echo "ALTER TABLE msk_fin_2 ADD COLUMN grassrgb varchar(11)" |
db.execute
#Aggiungo i colori
echo "UPDATE msk_fin_2 SET grassrgb='000:255:000' WHERE recl = 1" |
db.execute
echo "UPDATE msk_fin_2 SET grassrgb='255:255:000' WHERE recl = 2" |
db.execute
echo "UPDATE msk_fin_2 SET grassrgb='255:000:000' WHERE recl = 3" |
db.execute
################################ RASTER ################################
#raster finale
v.to.rast input=msk_fin output=msk_fin use=attr type=area layer=2
column=recl value=0 rgbcolumn=grassrgb
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
r.report map=msk_fin units=k
+-----------------------------------------------------------------------------+
| RASTER MAP CATEGORY REPORT
|
|LOCATION: tosc_utm_ed50 Thu Jan 1 22:18:19 2009|
|-----------------------------------------------------------------------------|
| north: 4925313
| east:771799
|
|REGION south: 4678313
west:554799 |
| res: 1000 res: 1000
|
|-----------------------------------------------------------------------------|
|MASK:none
|
|-----------------------------------------------------------------------------|
|MAP: Labels (msk_fin in tos_bos)
|
|-----------------------------------------------------------------------------|
| Category Information |
square |
|#|description |
kilometers|
|-----------------------------------------------------------------------------|
|1| . . . . . . . . . . . . . . . . . . . . . . . . . . . .| 7399.000|
|2| . . . . . . . . . . . . . . . . . . . . . . . . . . . .| 4158.000|
|3| . . . . . . . . . . . . . . . . . . . . . . . . . . . .| 2822.000|
|*|no data. . . . . . . . . . . . . . . . . . . . . . . . .|39,220.000|
|-----------------------------------------------------------------------------|
|TOTAL |53,599.000|
+-----------------------------------------------------------------------------+
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
###############################################################################
###############################################################################
################################## NOTA BENE ###########################
###############################################################################
###############################################################################
#se in un monitor visualizzo il raster msk_fin e il vettoriale msk_tos
#come boundary i due file non si sovrappongono correttamente.
#Si dovrebbe trattare di un errore di visualizzazione (uno
#shift) di uno dei due - raster o vettoriale - rispetto all'altra
#tipologia, infatti:
d.mon start=x1
d.rast map=msk_fin
d.vect map=msk_tos type=boundary
d.zoom (vari zoom)
#mostra una non perfetta sovrapposizione tra vettoriale e raster
#(generato da una elaborazione di tale vettoriale) ma se:
d.erase
g.region -d
#riporto il raster nuovamente a vettoriale
r.to.vect input=msk_fin output=controllo feature=area
d.vect -c map=controllo width=2 color=red
d.vect map=msk_tos type=boundary
d.zoom (vari zoom)
#in questo caso i due vettoriali si sovrappongono perfettamente!
###############################################################################
###############################################################################
################################ FINE NOTA #############################
###############################################################################
###############################################################################
--ho lavorato con
GRASS 6.3.1svn
Debian Lenny aggiornata al 2009-01-02 montato su mvware che gira su
Debian Stable aggiornata al 2009-01-02
Maggiori informazioni sulla lista
Gfoss