User:Hanoj/grass
Jump to navigation
Jump to search
How to use GRASS for vectorize raster cadastral boundaries.
Introduction
- http://www.czso.cz/csu/rso.nsf/i/schema_soustavy
- count of cadastral parish(+/- 13 027): http://www.czso.cz/csu/rso.nsf/i/katastralni_uzemi
- agree with OSM license
- WMS GetCapabilities, example map
- 4 layers in category "prehledky" (limited resolution):
- prehledka_kat_uz-linie = (black logo + green line)
- prehledka_kat_uz-texty = (black logo + green text + blue point)
- 2 working layer (???)
Tagging
- source=cuzk:prehledky
- boundary=administrative
- date:import_state=2009-06-21
Plan
- vectorize line of boundaries (06/2009)
- vectorize point of annotation
- OCR annotation and join to point (TODO)
- combine 1+3 and make relation? (TODO)
- hand fix bug (TODO)
1. Line vectorizing
Bug of result
- sirotci na liniich o delce cca 2m (modul r.to.vect)
- trojuhelniky o hrane cca 2m (modul r.to.vect)
- prerusene linie (po rastrovem logu)
- generalizace oproti rastrove predloze s posunem do 5 metru
Procedure
- Download WMS tiles via cuzk2tile.py (resolution 2m/px, bbox=, service WMS CUZK) (need package python 2.6, python-gdal)
old fix bug (not using) <!-# for i in *.pgw; do mv $i $i.bak; sed -e "4,4s/^2.0/-2.0/" $i.bak > $i; rm $i.bak; done-->- create or download GRASS S-JTSK region
- create GRASS new location
- run in GRASS batch script (tested in GRASS 6.4.0rc4, Ubuntu 9.04)
# clear cumulative VAR
mapvc=""
#list all PNG map tile in actual dir
for i in cuzk_km1-2-*x*.png; do
#use only non-plain tile, more than 8628 bytes
size=`stat -c "%s" $i`
if [ $size -gt 8628 ]; then
#name of maps for parts of procedure
map=`echo $i | sed -e "s/\.png//" | sed -e "s/-//g" | sed -e "s/cuzk_km12/km/"`
mapt="t"$map
mapv="v"$map
mapvc=$mapvc$mapv","
filePNG=$i
fileTIF=$map".tif"
echo "$map"
#r.in.gdal not support worldfile for PNG file; trans undirected to TIFF
gdalwarp $filePNG $fileTIF
#import to GRASS
r.in.gdal -o --overwrite --quiet input=$fileTIF output=$map
#region set to current map
g.region --quiet rast=$map
#mast set to green color with data
r.mapcalc MASK="if(r#$map == 0 && g#$map == 255 && b#$map == 0, 1, null())"
#simplify rastr
r.thin --overwrite --quiet input=$map output=$mapt
#vectorize
r.to.vect --overwrite --quiet input=$mapt output=$mapv
#remove temp raster mask, layers and files
r.mask -r
g.remove -f --quiet rast=$map,$mapt
rm $fileTIF
fi
done
#remove latest mask
r.mask -r --quiet
#region set to default, global
g.region -d -p
#merge all vector layers
v.patch --overwrite input=$mapvc output=kucr
#generalize
v.generalize --overwrite input=kucr@ku output=kucrg@ku threshold=5 \
look_ahead=7 reduction=50 slide=0.5 angle_thresh=5 degree_thresh=0 \
closeness_thresh=0.5 betweeness_thresh=0.5 alpha=2 beta=2 iterations=2
#clean topology
v.clean --overwrite input=kucrg output=kucrgc type=line \
tool=break,snap,rmdangle,chdangle,rmbridge,chbridge,rmdupl,rmdac,bpol,\
prune,rmarea,rmline,rmsa
#add table and length column
v.db.addtable map=kucrgc
v.db.addcol map=kucrgc columns="linelength integer"
v.to.db map=kucrgc type=line option=length units=me columns=linelength
#remove temp vector layers
g.remove vect=$mapvc
#export layer to ESRI Shapefile
v.out.ogr input=kucrgc dsn=kucr_jtsk.shp
4. ogr2ogr transformation S-JTSK to WGS84
ogr2ogr -s_srs "+proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813972222222 \ +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=greenwich +units=m +no_defs \ +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56" -t_srs "epsg:4326" -f "ESRI Shapefile" kucr.shp kucr_jtsk.shp
5. convert perl shp2osm.pl (require Geo::ShapeFile) to OSM
perl ./shp2osm.pl kucr > kucr.osm
6. view in josm-ng
java -Xmx512m -jar /home/enemy/.josm/josm/josm-ng.jar
2+3. Grab annotation and OCR
- získat souřadnice modrých bodů v S-JTSK jako polohu popisku každého k.ú. get_points.py => kat_body.tar.gz
- na základě bodů postahovat jednotlivé výřezy points2text.py => např. bod-2_-804738,-980170.png
- na výřezy pustit OCR
- vytvořit tabulku "x, y, číslo k.ú., název k.ú."