User:Hanoj/grass

From OpenStreetMap Wiki
Jump to: navigation, search

How to use GRASS for vectorize raster cadastral boundaries.

Introduction

Tagging

  • source=cuzk:prehledky
  • boundary=administrative
  • date:import_state=2009-06-21

Plan

  1. vectorize line of boundaries (06/2009)
  2. vectorize point of annotation
  3. OCR annotation and join to point (TODO)
  4. combine 1+3 and make relation? (TODO)
  5. 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

  1. Download WMS tiles via cuzk2tile.py (resolution 2m/px, bbox=, service WMS CUZK) (need package python 2.6, python-gdal)
  2. 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-->
  3. create or download GRASS S-JTSK region
  4. create GRASS new location
  5. 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

  1. získat souřadnice modrých bodů v S-JTSK jako polohu popisku každého k.ú. get_points.py => kat_body.tar.gz
  2. na základě bodů postahovat jednotlivé výřezy points2text.py => např. bod-2_-804738,-980170.png
  3. na výřezy pustit OCR
  4. vytvořit tabulku "x, y, číslo k.ú., název k.ú."