Merge pull request #75 from GenericMappingTools/classifications

Add functions to do assisted classification.

Project.toml

@@ -9,6 +9,7 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Dates  = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
+DecisionTree = "7806a523-6efd-50cb-b5f6-3fa6f1930dbb"
 
 [weakdeps]
 SatelliteToolboxTle = "7ff27aeb-5fff-4337-a9ee-a9fe6b7ed35e"
@@ -34,3 +35,4 @@ PrecompileTools = "1.0"
 SatelliteToolboxTle = "1"
 SatelliteToolboxPropagators = "0.3"
 SatelliteToolboxTransformations = "0.1"
+DecisionTree = "0.12"

src/RemoteS.jl

@@ -1,7 +1,8 @@
 module RemoteS
 
-using GMT, Printf, Statistics, Dates#, Requires
+using GMT, Printf, Statistics, Dates
 using PrecompileTools
+using DecisionTree
 
 const SCENE_HALFW = Dict("AQUA" => 1163479, "TERRA" => 1163479, "LANDSAT8" => 92500)	# half widths
 
@@ -11,8 +12,9 @@ end
 
 export
 	cutcube, subcube, dn2temperature, dn2radiance, dn2reflectance, reflectance_surf, grid_at_sensor, truecolor,
-	clg, clre, evi, evi2, gndvi, mndwi, mtci, mcari, msavi, mbri, ndvi, ndwi, ndwi2, ndrei1,
+	clg, clre, evi, evi2, gndvi, mndwi, mtci, mcari, msavi, nbri, ndvi, ndwi, ndwi2, ndrei1,
 	ndrei2, satvi, savi, slavi,
+	classify, train_raster,
 	clip_orbits, findscenes, sat_scenes, sat_tracks, reportbands
 
 include("grid_at_sensor.jl")

src/spectral_indices.jl

@@ -208,6 +208,9 @@ Normalised Burn Ratio Index. Garcia 1991
 NBRI = (nir - swir2) / (nir + swir2)
 """
 nbri(nir, swir2; kw...) = sp_indices(swir2, nir; index="NBRI", kw...)
+nbri(cube::String;
+     bands::Vector{Int}=Int[], layers::Vector{Int}=Int[], bandnames::Vector{String}=String[], kw...) =
+	helper_si_method(cube, "NBRI"; bands=bands, layers=layers, bandnames=bandnames, defbandnames=["nir", "swir2"], kw...)
 nbri(cube::Union{GMT.GMTimage{UInt16, 3}, AbstractArray{<:AbstractFloat, 3}};
      bands::Vector{Int}=Int[], layers::Vector{Int}=Int[], bandnames::Vector{String}=String[], kw...) =
 	helper_si_method(cube, "NBRI"; bands=bands, layers=layers, bandnames=bandnames, defbandnames=["nir", "swir2"], kw...)

src/utils.jl

@@ -752,3 +752,79 @@ function reflectance_surf(fname::String; band::Int=0, mtl::String="", save::Stri
 	(save != "") && gdaltranslate(G, dest=save)
 	return (save != "") ? nothing : G
 end
+
+# ----------------------------------------------------------------------------------------------------------
+"""
+    I = classify(cube::GItype, train::Union{Vector{<:GMTdataset}, String}) -> GMTimage
+
+- `cube`: The cube wtih band data to classify.
+- `train`: A vector of GMTdatasets or a file name of one containing the polygons used to train the model.
+   NOTE: The individual datasets MUST have associated an attribute called "class" containing the class name as a string.
+   This can be achieved for text data in the form of a GMT multi-segment file (one where segments are separated by the '>'
+   symbol) if the multi-segment separator line contains the text ``Attrib(class=name)``
+
+Returns an image with the classification results where each class name was assigned a different integer number.
+That colorized image can plotted with ``viz(I, colorbar=true)``.
+"""
+function classify(cube::GItype, train::Union{Vector{<:GMTdataset}, String})
+	model, classes = train_raster(cube, train)
+	I = classify(cube, model)
+	cpt = makecpt(cmap=:categorical, range=classes);
+	image_cpt!(I, cpt)
+	return I
+end
+
+# ----------------------------------------------------------------------------------------------------------
+"""
+    I = classify(cube::GItype, model; class_names::Union{String, Vector{String}}="") -> GMTimage
+
+- `cube`: The cube wtih band data to classify.
+- `model`: The trained model obtained from the `train_raster` function.
+- `class_names`: A vector of strings with the class names to be used in the categorical colorbar or a
+   comma separated single with those class names. The number of class names must match the number used
+   when training the model with `train_raster`.
+"""
+function classify(cube::GItype, model; class_names::Union{String, Vector{String}}="")
+	nr, nc = size(cube)[1:2];
+	mat = Vector{UInt8}(undef, nr*nc);
+	t = permutedims(cube.z, (1,3,2));
+	i1 = 1;	i2 = nr;
+	for k = 1:nc			# Loop over columns
+		mat[i1:i2] = DecisionTree.predict(model, Float64.(t[:,:,k]))
+		i1 = i2 + 1
+		i2 = i1 + nr - 1
+	end
+	I = mat2img(reshape(mat, nc,nr), cube);
+	(class_names == "") && return I			# No class names, no CPT
+	classes = isa(class_names, Vector) ? join(class_names, ",") : class_names
+	cpt = makecpt(cmap=:categorical, range=classes);
+	image_cpt!(I, cpt)
+	return I
+end
+
+# ----------------------------------------------------------------------------------------------------------
+"""
+    model, classes = train_raster(cube::GItype, train::Union{Vector{<:GMTdataset}, String}; np::Int=0, density=0.1)
+
+- `cube`: The cube wtih band data to classify.
+- `train`: A vector of GMTdatasets or a file name of one containing the polygons used to train the model.
+   NOTE: The individual datasets MUST have associated an attribute called "class" containing the class name as a string.
+   This can be achieved for text data in the form of a GMT multi-segment file (one where segments are separated by the '>'
+   symbol) if the multi-segment separator line contains the text ``Attrib(class=name)``
+- `np`: Number of points per polygon to be determined by ``randinpolygon``
+- `density`: Alternative to `np`. See also the help of the ``randinpolygon`` function.
+
+Returns the trained model and the class names.
+"""
+function train_raster(cube::GItype, train::Union{Vector{<:GMTdataset}, String}; np::Int=0, density=0.1)
+	samples = isa(train, String) ? gmtread(train) : train
+	pts = randinpolygon(samples, np=np, density=density);
+	LCsamp = grdinterpolate(cube, S=pts, nocoords=true);
+	features = GMT.ds2ds(LCsamp);
+	labels = parse.(UInt8, vcat([fill(LCsamp[k].attrib["id"], size(LCsamp[k],1)) for k=1:length(LCsamp)]...));
+
+	model = DecisionTree.DecisionTreeClassifier(max_depth=3);
+	DecisionTree.fit!(model, features, labels)
+	classes = join(unique(GMT.make_attrtbl(samples, false)[1][:,1]), ",")
+	return model, classes
+end