From OpenStreetMap Wiki
Jump to: navigation, search

How to use GRASS for vectorize raster cadastral boundaries.



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


  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
  • trojuhelniky o hrane cca 2m (modul
  • prerusene linie (po rastrovem logu)
  • generalizace oproti rastrove predloze s posunem do 5 metru


  1. Download WMS tiles via (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

#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/"`
               echo "$map" 

      not support worldfile for PNG file; trans undirected to TIFF
               gdalwarp $filePNG $fileTIF

               #import to GRASS
      -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

      --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


#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

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 \

#add table and length column
v.db.addtable map=kucrgc 
v.db.addcol map=kucrgc columns="linelength integer" 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 (require Geo::ShapeFile) to OSM

perl ./ 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.ú. => kat_body.tar.gz
  2. na základě bodů postahovat jednotlivé výřezy => např. bod-2_-804738,-980170.png
  3. na výřezy pustit OCR
  4. vytvořit tabulku "x, y, číslo k.ú., název k.ú."