Raster map generator

script/classify-mbes-data_in_separate_files.bash

@@ -5,7 +5,7 @@ if [[ $# -ne 1 ]]; then
 	exit 1
 fi
 
-if [[ "$(basename $(pwd))" -ne "BenthicClassifier" ]]; then
+if [[ "$(basename $(pwd))" != "BenthicClassifier" ]]; then
 	echo "please execute script from project root"
 	exit 1
 fi
@@ -22,6 +22,7 @@ do
 	fbname=${fname%.*}
 	#echo $dir/classifications/$fbname.xyzc
 	echo x, y, z, boosting class, gmm class > $dir/classifications/$fbname.csv
-	python3 src/apply_both_models.py trained.model gmm_trained.model $FILE >> $dir/classifications/$fbname.csv
+	#python3 src/apply_both_models.py trained.model gmm_trained.model $FILE >> $dir/classifications/$fbname.csv
+	python3 src/apply-model.py trained.model $FILE >> $dir/classifications/$fbname.csv
 done
 

script/classify-mbes-data_in_single_file.bash

@@ -5,7 +5,7 @@ if [[ $# -ne 1 ]]; then
 	exit 1
 fi
 
-if [[ "$(basename $(pwd))" -ne "BenthicClassifier" ]]; then
+if [[ "$(basename $(pwd))" != "BenthicClassifier" ]]; then
 	echo "please execute script from project root"
 	exit 1
 fi
@@ -15,7 +15,7 @@ dir=$1
 echo $dir
 mkdir -p $dir/classification
 
-echo x, y, z, boosting class, gmm class > $dir/classification/classified_both_model.xyzc
+echo x, y, z, boosting class, gmm class > $dir/classification/classified_both_model.csv
 for FILE in $dir*.hackel;
 do
 	echo $FILE

script/gen_hackel.bash

@@ -5,7 +5,7 @@ if [[ $# -ne 3 ]]; then
 	exit 1
 fi
 
-if [[ "$(basename $(pwd))" -ne "BenthicClassifier" ]]; then
+if [[ "$(basename "$PWD")" != "BenthicClassifier" ]]; then
 	echo "please execute script from project root"
 	exit 1
 fi

script/gen_raster.bash

@@ -1,94 +1,25 @@
 #!/bin/bash
 
-gen_vrt_gbx(){
-	classificationsDir=$1
-	xyzFile=$2
-	vrtDir=$3
-	vrtFilename=$4
-#	echo "1 " $xyzFile
-#	echo "2 " ${vrtFilename%.*}
-	output=$vrtDir/$vrtFilename
-	echo "create virtual file format = " $output
-
-
-	echo "<OGRVRTDataSource>
-		<OGRVRTLayer name=\"${vrtFilename%.*}\">
-		    <SrcDataSource>$classificationsDir/$xyzFile</SrcDataSource>
-		    <GeometryType>wkbPoint</GeometryType>
-		    <GeometryField encoding=\"PointFromColumns\" x=\"x\" y=\"y\" z=\"boosting class\"/>
-		</OGRVRTLayer>
-	</OGRVRTDataSource>" > $output
-
-}
-
-gen_vrt_gmm(){
-	classificationsDir=$1
-	xyzFile=$2
-	vrtDir=$3
-	vrtFilename=$4
-	echo "1 " $xyzFile
-	echo "2 " ${vrtFilename%.*}
-	output=$vrtDir/$vrtFilename
-	echo "create virtual file format = " $output
-
-
-	echo "<OGRVRTDataSource>
-		<OGRVRTLayer name=\"${vrtFilename%.*}\">
-		    <SrcDataSource>$classificationsDir/$xyzFile</SrcDataSource>
-		    <GeometryType>wkbPoint</GeometryType>
-		    <GeometryField encoding=\"PointFromColumns\" x=\"x\" y=\"y\" z=\"gmm class\"/>
-		</OGRVRTLayer>
-	</OGRVRTDataSource>" > $output
-
-}
-
-if [[ $# -ne 4 ]]; then
-	echo "usage : bash gen_raster.bash classificationsDirectory resolutionX resolutionY epsgCode"
+if [[ "$(basename "$PWD")" != "BenthicClassifier" ]]; then
+	echo "please execute script from project root"
 	exit 1
 fi
 
-if [[ "$(basename $(pwd))" -ne "BenthicClassifier" ]]; then
-	echo "please execute script from project root"
+if [[ $# -ne 2 ]]; then
+	echo "usage : bash hackel_batch.bash directory outputDirectory"
 	exit 1
 fi
 
-
-
 dir=$1
-vrtGbxDir=$dir/vrts_gbx
-vrtGmmDir=$dir/vrts_gmm
-rasterGbxDir=$dir/rasters_gbx
-rasterGmmDir=$dir/rasters_gmm
-mkdir -p $rasterGbxDir
-mkdir -p $rasterGmmDir
-mkdir -p $vrtGbxDir
-mkdir -p $vrtGmmDir
-
-resolutionX=$2
-resolutionY=$3
-epsg=$4
-
-# gbx
-for FILE in $(ls $dir/*.csv);
-do
-	FILE=$(basename "$FILE")
-	vrtFilename=${FILE%.*}.vrt
-	echo $vrtFilename
-	gen_vrt_gbx $dir $FILE $vrtGbxDir $vrtFilename
-	gdal_grid -outsize $resolutionX $resolutionY -a_srs EPSG:$epsg -a invdist -of GTiff -ot Float64 -l ${vrtFilename%.*} $vrtGbxDir/$vrtFilename $rasterGbxDir/$FILE.tiff --config GDAL_NUM_THREADS ALL_CPUS
-
-done
+outDir=$2
 
+mkdir -p $outDir
 
-# gmm
 for FILE in $(ls $dir/*.csv);
 do
-	FILE=$(basename "$FILE")
-	vrtFilename=${FILE%.*}.vrt
 	echo $FILE
-	gen_vrt_gmm $dir $FILE $vrtGmmDir $vrtFilename
-	gdal_grid -outsize $resolutionX $resolutionY -a_srs EPSG:$epsg -a invdist -of GTiff -ot Float64 -l ${vrtFilename%.*} $vrtGmmDir/$vrtFilename $rasterGmmDir/$FILE.tiff --config GDAL_NUM_THREADS ALL_CPUS
-
+	FILENAME=$(basename -s .csv $FILE)
+	python3 src/generate_raster_map.py $FILE $outDir/$FILENAME.tif $FILENAME
 done
 
 

script/gen_table_figure.bash

@@ -5,7 +5,7 @@ if [[ $# -ne 1 ]]; then
 	exit 1
 fi
 
-if [[ "$(basename $(pwd))" -ne "BenthicClassifier" ]]; then
+if [[ "$(basename $(pwd))" != "BenthicClassifier" ]]; then
 	echo "please execute script from project root"
 	exit 1
 fi

script/pipeline.bash

@@ -29,7 +29,7 @@ train_models(){
 
 classify_data(){
 	hackelDirectory=$1
-	bash script/classify-mbes-data_in_separate_files.bash hackelDirectory
+	bash script/classify-mbes-data_in_separate_files.bash $hackelDirectory
 }
 
 if [[ $# -ne 3 ]]; then
@@ -38,7 +38,7 @@ if [[ $# -ne 3 ]]; then
 fi
 
 # all paths are relative to project root
-if [[ "$(basename $(pwd))" -ne "BenthicClassifier" ]]; then
+if [[ "$(basename $(pwd))" != "BenthicClassifier" ]]; then
 	echo "please execute script from project root"
 	exit 1
 fi
@@ -51,10 +51,10 @@ hackelRadius=$3
 
 ############ SCRIPT ###################
 
-gen_hackel $mbesDir $hackelDir $hackelRadius
-combine_file $hackelDir $HOME/allSurfaces.hackel
+#gen_hackel $mbesDir $hackelDir $hackelRadius
+#combine_file $hackelDir $HOME/allSurfaces.hackel
 
-train_models mbes-training-data.txt $HOME/allSurfaces.hackel
+#train_models mbes-training-data.txt $HOME/allSurfaces.hackel
 classify_data $hackelDir
 
 

src/compare-model.py

@@ -29,7 +29,7 @@ def read_classification_file(filename):
 	sys.exit(1)
 
 
-#sys.stderr.write(["0 : 19", "1 : blocks", "2 : cobble", "3 : gravels", "4 : rocky", "5 : sands", "6 : sandy mud"]) 
+#sys.stderr.write(["0 : blocks", "1 : cobble", "2 : gravels", "3 : rocky", "4 : sands", "5 : sandy mud"]) 
 
 filePath = sys.argv[1]
 classifiedData = read_classification_file(filePath)

src/generate_raster_map.py

@@ -0,0 +1,62 @@
+import pygmt
+import sys, csv
+import numpy as np
+import pandas as pd
+
+colormap = {0:"red", 1:"grey", 2:"blue", 3:"pink", 4:"orange", 5:"black", 6:"green", 7:"purple", 8:"yellow", 9:"brown"}
+gbxClassMap = {0:"blocks", 1:"cobble", 2:"gravels", 3:"rocky", 4:"sand", 5:"sandy mud"}
+
+if len(sys.argv) != 4:
+	sys.stderr.write("Usage: python3 generate_raster_map.py inputData Map_title\n")
+	sys.exit(1)
+
+inputFilePath = sys.argv[1]
+outputFilename = sys.argv[2]
+mapTitle = sys.argv[3]
+
+
+data = pd.read_csv(inputFilePath, delimiter=",", dtype=float, header=0)
+print("[+] Read {} lines".format(len(data.index)))
+
+
+minLat = data.min(axis=0)[0]
+minLon = data.min(axis=0)[1]
+maxLat = data.max(axis=0)[0]
+maxLon = data.max(axis=0)[1]
+
+gbxGroups = data.groupby("boosting class")
+
+
+fig = pygmt.Figure()
+fig.basemap(region=[minLon - 0.005, maxLon + 0.005, minLat - 0.005, maxLat + 0.005], 
+			frame=["a", "+t{}".format(mapTitle)],
+			projection="M6i"
+)
+
+fig.coast(land="lightblue", water="lightblue")
+
+for classID in gbxGroups.groups:
+	df = gbxGroups.get_group(classID)
+	fig.plot(x=df["longitude"], y=df["latitude"], style="c0.03c", fill=colormap[classID], label=gbxClassMap[classID]+"+S0.25c")
+
+
+fig.legend(box=True)
+fig.image(
+    imagefile="https://cidco.ca/themes/cidco/css/img/logo_cidco.png",
+    position="jBR+w2+o0.2c",
+    box=True,
+)
+
+with fig.inset(
+	position="jBL",
+	region=[round(minLon-0.3,3), round(maxLon+0.3,3), round(minLat-0.3,3), round(maxLat+0.3,3)],
+	box="+gwhite+p1p",
+	projection="M2c"
+):
+	fig.coast(shorelines="thinnest", land="green", water="lightblue")
+	fig.plot(data=[[minLon, minLat, maxLon, maxLat]], style="r+s", pen="1p,red")
+
+#fig.show()
+fig.savefig("{}".format(outputFilename), crop=True, dpi=1200)
+
+