User Tools

Site Tools


wiki:python:segmentwatershed

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

wiki:python:segmentwatershed [2018/06/02 18:05] (current)
Line 1: Line 1:
 +====== Segmenting the Global Human Settlement Layer ======
 +
 +{{ :​wiki:​python:​screenshot_2018-06-02_07-03-31.png?​400 | }}
 +
 +== Download the file and setting up working directory ​ ==
 +
 +<code bash >
 +mkdir $HOME/​GHSLseg
 +cd $HOME/​GHSLseg
 +wget http://​www.spatial-ecology.net/​ost4sem/​exercise/​GHSLseg/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0.tif
 +
 +</​code>​
 +
 +== Thresholding the GHSL and clumping each bin-unit ==
 +
 +<code bash| sc01_bin_clump.sh>​
 +
 +
 +export DIR=$HOME/​GHSLseg
 +
 +oft-calc -ot Byte $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0.tif $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin.tif <<EOF
 +1
 +#1 0.05 + 10 * 1 -
 +EOF
 +rm  $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin.tif.eq
 +
 +# colored ​ the bin file 
 +pkcreatect -min 0 -max 9 > /​tmp/​color.txt ​
 +pkcreatect ​  -co COMPRESS=DEFLATE -co ZLEVEL=9 ​ -ct  /​tmp/​color.txt ​ -i   ​$DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin.tif -o   ​${DIR}/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin_ct.tif
 +gdal_edit.py -a_nodata -1   ​${DIR}/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin_ct.tif
 +mv ${DIR}/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin_ct.tif ${DIR}/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin.tif
 +
 +echo start the clump operation ​
 +
 +rm -fr ${DIR}/​grassdb
 +mkdir ${DIR}/​grassdb
 +
 +echo  1 2 3 4 5 6 7 8 9   | xargs -n 1 -P 2 bash -c $' ​
 +
 +MIN=$( echo $1 - 0.5 | bc )
 +BIN=$1
 +
 +echo masking the bin 
 +
 +pkgetmask ​ -co COMPRESS=DEFLATE -co ZLEVEL=9 ​ -min $MIN  -max 9.5  -data 1 -nodata 0 -ct  /​tmp/​color.txt ​ -i ${DIR}/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin.tif -o ${DIR}/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin$BIN.tif ​
 +gdal_edit.py -a_nodata 0   ​${DIR}/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin$BIN.tif  ​
 +
 +# create the new location and exit
 +
 +rm -fr  ${DIR}/​grassdb/​loc_clump$BIN
 +
 +grass72 -e -text  -c  ${DIR}/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin$BIN.tif ${DIR}/​grassdb/​loc_clump$BIN ​
 +
 +# set up grass variables
 + 
 +echo "​GISDBASE:​ ${DIR}/​grassdb" ​                 >  $HOME/​.grass7/​rc$$
 +echo "​LOCATION_NAME:​ loc_clump$BIN" ​             >> $HOME/​.grass7/​rc$$
 +echo "​MAPSET:​ PERMANENT" ​                        >>​ $HOME/​.grass7/​rc$$
 +echo "GUI: text" ​                                >>​ $HOME/​.grass7/​rc$$
 +echo "​GRASS_GUI:​ wxpython" ​                      >>​ $HOME/​.grass7/​rc$$
 +
 +export GISBASE=/​usr/​lib/​grass72
 +export PATH=$PATH:​$GISBASE/​bin:​$GISBASE/​scripts
 +export LD_LIBRARY_PATH="​$GISBASE/​lib"​
 +export GISRC=$HOME/​.grass7/​rc$$
 +export GRASS_LD_LIBRARY_PATH="​$LD_LIBRARY_PATH"​
 +export PYTHONPATH="​$GISBASE/​etc/​python:​$PYTHONPATH"​
 +export MANPATH=$MANPATH:​$GISBASE/​man
 +
 +export GIS_LOCK=$$
 +
 +g.gisenv ​
 +r.in.gdal input=${DIR}/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin$BIN.tif output=GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin$BIN ​ --o
 +
 +r.clump -d  --overwrite ​   input=GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin$BIN ​    ​output=GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin${BIN}_clump
 +r.colors -r map=GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin$BIN
 +
 +r.out.gdal -m  nodata=0 --overwrite -f -c createopt="​COMPRESS=DEFLATE,​ZLEVEL=9"​ format=GTiff type=UInt32 ​ input=GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin${BIN}_clump ​ output=${DIR}/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin${BIN}_clump.tif ​
 +rm -rf ${DIR}/​grassdb/​loc_clump$BIN ​ ${DIR}/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin${BIN}_clump.tif.aux.xml ​
 +
 +' _  ​
 +
 +rm -fr ${DIR}/​grassdb
 +
 +</​code>​
 +
 +{{ :​wiki:​python:​bin.png?​400 |}}
 +
 +== Re-class the bin-unit to have consecutive ID values among the all levels. == 
 +
 +<code bash| sc02_bin_clump_reclass.sh>​
 +
 +# re-class the bin clump and then overlay ​
 +export DIR=$HOME/​GHSLseg
 +
 +# will add to 1 ... so the overall will start from 1 
 +
 +lastmaxb=0 ​  
 +
 +for BIN in 1 2 3 4 5 6 7 8 9 ; do 
 +
 +echo hist for BIN $BIN
 +pkstat -hist -i ${DIR}/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin${BIN}_clump.tif | grep -v " 0" | awk -v lastmaxb=$lastmaxb '{ if ($1==0) { print $1,0} else {lastmaxb=1+lastmaxb;​ print $1,​lastmaxb}}'​ > $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin${BIN}_clump.txt
 +
 +lastmaxb=$(tail -1 $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin${BIN}_clump.txt | awk '{ print $2 }')
 +
 +echo reclass for BIN $BIN
 +pkreclass -ot UInt32 -code $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin${BIN}_clump.txt ​ -co COMPRESS=DEFLATE -co ZLEVEL=9 -i ${DIR}/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin${BIN}_clump.tif -o ${DIR}/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin${BIN}_clump.tif
 +
 +gdal_edit.py -a_nodata 0 ${DIR}/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin${BIN}_clump.tif
 +rm $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin${BIN}_clump.txt
 +done 
 +
 +# combine all the bin clump and get the highst value 
 +
 +gdalbuildvrt -separate -overwrite ${DIR}/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin_clump.vrt ${DIR}/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin?​_clump.tif
 +
 +oft-calc -of GTiff -ot UInt32 ${DIR}/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin_clump.vrt ${DIR}/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin_clump_tmp.tif <<EOF
 +1
 +#1 #2 #3 #4 #5 #6 #7 #8 #9 M M M M M M M M
 +EOF
 +gdal_translate -a_nodata 0 -ot UInt32 -co COMPRESS=DEFLATE -co ZLEVEL=9 ${DIR}/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin_clump_tmp.tif ${DIR}/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin1-9_clump.tif
 +
 +rm ${DIR}/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin_clump_tmp.tif ${DIR}/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin_clump.vrt $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin_clump_tmp.tif.eq
 +
 +</​code>​
 +
 +== Identify the peaks ==
 +
 +<code bash| sc03_bin_peak_clump.sh>​
 +
 +export DIR=$HOME/​GHSLseg
 +
 +pkcreatect -min 0 -max 1 > /​tmp/​color.txt ​
 +pkcreatect -co COMPRESS=DEFLATE -co ZLEVEL=9 -ot Byte -nodata -1 -min 0 -max 1 -i $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin9.tif -o  $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin9_clean.tif
 +gdal_edit.py -a_nodata -1  $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin9_clean.tif
 +
 +# change the nodata to -1 to allow the gdacalc to consider ​
 +echo  0 1 2 3 4 5 6 7 8 9 | xargs -n 1 -P 8 bash -c $'
 +gdal_edit.py -a_nodata -1  $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin$1.tif
 +' _
 +
 +echo 2 3 4 5 6 7 8 9 | xargs -n 1 -P 2 bash -c $'
 +
 +    bin1=$1
 +    bin2=$(bc <<<​ "​$bin1-1"​)
 +    echo $bin1
 +        rm -f $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin${bin1}+${bin2}.tif ​
 +        gdal_calc.py --calc="​A+B"​ --co="​COMPRESS=LZW"​ --overwrite --NoDataValue=-1 \
 +            -A $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin${bin1}.tif \
 +            -B $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin${bin2}.tif \
 +            --outfile=$DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin${bin1}+${bin2}.tif
 +
 +        oft-stat -mm --noavg --nostd \
 +            -i  $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin${bin1}+${bin2}.tif \
 +            -um $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin${bin2}_clump.tif \
 +            -o  $DIR/​zonal_stats_${bin2}.txt
 +
 +        awk \'​{if($4==1) { print $1,1 }  else { print $1,0 }}\' $DIR/​zonal_stats_${bin2}.txt > $DIR/​code_${bin2}.txt
 +
 +        pkreclass -co COMPRESS=DEFLATE -co ZLEVEL=9 ​ -ot Byte -nodata -1 -ct /​tmp/​color.txt \
 +            -i $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin${bin2}_clump.tif \
 +            -o $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin${bin2}_clean.tif \
 +            --code $DIR/​code_${bin2}.txt
 +
 +        rm -f $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin${bin1}+${bin2}.tif
 +        rm -f $DIR/​zonal_stats_${bin2}.txt
 +        rm -f $DIR/​code_${bin2}.txt
 +
 +' _
 +
 +#  only search for peaks with bin values larger or equal to 2 
 +  ​
 +gdalbuildvrt -overwrite ​ -separate $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_peak.vrt $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin{2,​3,​4,​5,​6,​7,​8,​9}_clean.tif
 +
 +oft-calc -ot Byte $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_peak.vrt $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_peak_tmp.tif << EOF
 +1
 +#1 #2 #3 #4 #5 #6 #7 #8 + + + + + + +
 +EOF
 +
 +gdal_translate ​  -co COMPRESS=DEFLATE -co ZLEVEL=9 ​ -ot Byte   ​$DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_peak_tmp.tif $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_peak.tif
 +# rm  $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_peak_tmp.tif $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_peak.vrt $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_peak_tmp.tif.eq ​
 +
 +pksetmask -co COMPRESS=DEFLATE -co ZLEVEL=9 -m $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_peak.tif ​ -msknodata 0 -nodata 0 -i $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin.tif ​ -o $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin_peak.tif
 +
 +pksetmask -co COMPRESS=DEFLATE -co ZLEVEL=9 -m $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_peak.tif ​ -msknodata 0 -nodata 0 -i $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_bin1-9_clump.tif ​ -o $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_peak_clump.tif ​
 +
 +</​code> ​
 +
 +== Built up cost function using the GHSL and peaks == 
 +
 +<code bash| sc04_built-up_cost.sh>​
 +
 +
 +DIR=/​home/​selv/​spatial-ecology-codes/​wiki/​watershedsegment
 +
 +rm -f $DIR/​grassdb ; mkdir $DIR/​grassdb
 +
 +gdal_edit.py -a_nodata 0  $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_peak.tif
 +grass72 -e -text  -c  $DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_peak.tif ​  ​$DIR/​grassdb/​loc_cost
 +
 +# set up grass variables
 + 
 +echo "​GISDBASE:​ ${DIR}/​grassdb" ​        > ​ $HOME/​.grass7/​rc$$
 +echo "​LOCATION_NAME:​ loc_cost" ​         >> $HOME/​.grass7/​rc$$
 +echo "​MAPSET:​ PERMANENT" ​               >> $HOME/​.grass7/​rc$$
 +echo "GUI: text" ​                       >> $HOME/​.grass7/​rc$$
 +echo "​GRASS_GUI:​ wxpython" ​             >> $HOME/​.grass7/​rc$$
 +
 +export GISBASE=/​usr/​lib/​grass72
 +export PATH=$PATH:​$GISBASE/​bin:​$GISBASE/​scripts
 +export LD_LIBRARY_PATH="​$GISBASE/​lib"​
 +export GISRC=$HOME/​.grass7/​rc$$
 +export GRASS_LD_LIBRARY_PATH="​$LD_LIBRARY_PATH"​
 +export PYTHONPATH="​$GISBASE/​etc/​python:​$PYTHONPATH"​
 +export MANPATH=$MANPATH:​$GISBASE/​man
 +
 +export GIS_LOCK=$$
 +
 +g.gisenv ​
 +
 +r.in.gdal input=$DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_peak.tif output=peak ​ --o
 +r.in.gdal input=$DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0.tif ​     output=builtup ​   --overwrite ​  ​memory=2047 ​
 +
 +# complementary builtup
 +r.mapcalc " builtup_comp =  ( 1 - builtup )   " ​ --overwrite
 +
 +echo start to calculate the cost
 +
 +r.cost ​ -k input=builtup_comp output=builtup_cost start_raster=peak null_cost=-1 --overwrite ​
 +r.colors ​ -r map=builtup_comp
 +
 +r.out.gdal --overwrite nodata=-1 -c -f createopt="​COMPRESS=DEFLATE,​ZLEVEL=9"​ type=Float32 format=GTiff ​ input=builtup_cost ​ output=$DIR/​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_cost.tif
 +
 +rm -rf $DIR/​grassdb
 +
 +</​code>​
 +
 +== Segmentation analyse ==  ​
 +
 +
 +<code bash| sc05_watershed.sh>​
 +
 +# python -m pip install SimpleITK
 +# itk installation ​ https://​itk.org/​Wiki/​SimpleITK/​GettingStarted/​Visual_guide_to_building_on_Linux#​Step_4:​_Build_SimpleITK ​
 +
 +export DIR=$HOME/​GHSLseg
 +rm -f  $DIR/​watershed_line_nogeo.tif ​ $DIR/​watershed_poly_nogeo.tif  ​
 +
 +# export PATH="/​home/​selv/​anaconda3/​bin:​$PATH"​
 +
 +cd $DIR
 +python <<EOF
 +import os
 +
 +import SimpleITK as sitk
 +print("​importing image"​)
 +img  = sitk.ReadImage("​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_cost.tif"​ , sitk.sitkFloat32 ​ )  ​
 +
 +# to check img.GetPixelIDTypeAsString() 32-bit unsigned integer  ​
 +peak = sitk.ReadImage("​GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_peaka.tif"​)
 +
 +marker_img ​ = sitk.ConnectedComponent(peak)
 +
 +print("​start watershed"​)
 +ws_line ​ = sitk.MorphologicalWatershedFromMarkers( img, marker_img , markWatershedLine=True, ​ fullyConnected=True)
 +sitk.WriteImage( sitk.Cast( ws_line ​ ,  sitk.sitkFloat32 ​ ),        "​watershed_line_nogeo.tif"​ )
 +del(ws_line)
 +
 +ws_poly ​ = sitk.MorphologicalWatershedFromMarkers( img, marker_img , markWatershedLine=False,​ fullyConnected=True)
 +sitk.WriteImage( sitk.Cast( ws_poly ​ ,  sitk.sitkFloat32 ​ ),        "​watershed_poly_nogeo.tif"​ )
 +del(ws_poly)
 +EOF
 +
 +gdal_edit.py ​ -a_srs EPSG:​4326 ​ -a_ullr $(getCorners4Gtranslate GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_cost.tif ​ ) watershed_line_nogeo.tif ​
 +gdal_edit.py ​ -a_srs EPSG:​4326 ​ -a_ullr $(getCorners4Gtranslate GHS_BUILT_LDS2014_GLOBE_R2016A_54009_1k_v1_0_cost.tif ​ ) watershed_poly_nogeo.tif  ​
 +
 +</​code>​
 +
  
wiki/python/segmentwatershed.txt ยท Last modified: 2018/06/02 18:05 (external edit)