From 6b226489baf09bb557c5ceee74883984c53081cb Mon Sep 17 00:00:00 2001
From: Jan Mandel <Jan.Mandel@gmail.com>
Date: Fri, 9 Jan 2015 08:01:51 -0700
Subject: [PATCH] added scripts to read RSAC geotiff files

---
 other/Matlab/detection/geotiff2mat.py   | 44 +++++++++++++++++++
 other/Matlab/detection/showmod14.m      | 75 +++++++++++++++++++++++++++++++++
 other/Matlab/detection/showmod14files.m | 16 +++++++
 3 files changed, 135 insertions(+)
 create mode 100644 other/Matlab/detection/geotiff2mat.py
 create mode 100644 other/Matlab/detection/showmod14.m
 create mode 100644 other/Matlab/detection/showmod14files.m

diff --git a/other/Matlab/detection/geotiff2mat.py b/other/Matlab/detection/geotiff2mat.py
new file mode 100644
index 00000000..a4b9bd02
--- /dev/null
+++ b/other/Matlab/detection/geotiff2mat.py
@@ -0,0 +1,44 @@
+# convert geotiff file to matlab file & print basic info on screen too
+# usage: python geotiff2mat.py filename.tif
+# limitation: for now, one band only
+
+# tested on Python 3.4
+
+# arguments
+import sys
+
+file=sys.argv[1]
+infile=file
+outfile = file + '.mat'
+
+# read the geotiff file
+# see http://geoinformaticstutorial.blogspot.com/2012/09/reading-raster-data-with-python-and-gdal.html
+
+import gdal
+from gdalconst import *
+print('Reading GeoTiff file',infile)
+dataset = gdal.Open(infile, GA_ReadOnly)
+cols = dataset.RasterXSize
+rows = dataset.RasterYSize
+bands = dataset.RasterCount
+print('rows                ',rows)
+print('columns             ',cols)
+print('bands               ',bands)
+geotransform = dataset.GetGeoTransform()
+print('top left X          ',geotransform[0])
+print('W-E pixel resolution',geotransform[1])
+print('rotation, 0=North up',geotransform[2])
+print('top left Y          ',geotransform[3])
+print('rotation, 0=North up',geotransform[4])
+print('N-S pixel resolution',geotransform[5])
+band = dataset.GetRasterBand(1)
+data = band.ReadAsArray(0, 0, cols, rows)
+
+# write matlab file
+# using http://docs.scipy.org/doc/scipy-0.14.0/reference/tutorial/io.html
+import scipy.io as sio
+print('writing Matlab file ',outfile)
+sio.savemat(outfile,{'data' : data,'geotransform' : geotransform})
+
+exit()
+
diff --git a/other/Matlab/detection/showmod14.m b/other/Matlab/detection/showmod14.m
new file mode 100644
index 00000000..33b49469
--- /dev/null
+++ b/other/Matlab/detection/showmod14.m
@@ -0,0 +1,75 @@
+function showmod14(data,geo,s)
+% display MODIS14 data
+% from loading file converted by geotiff2mat.py
+% input 
+%   data     MODIS 14 data
+%   geo      geotransform from geotiff
+%   optional title string
+%
+% black = notprocessed
+% blue  = water
+% white = cloud
+% green = land
+% red   = fire (saturation is confidence level)
+
+
+cmap= [ ...
+    0   0   0    %0 not processed (missing input data), black
+    0   0   0    %1 not used, black
+    0   0   0    %2 not processed (other reason)
+    0   0   0.2    %3 water, dark blue
+    0.1 0.2 0.2    %4 cloud, purple gray
+    0   0.2   0    %5 non-fire clear land, green
+    0   0   0    %6 unknown
+    0.6 0.4 0.4    %7 low-confidence fire
+    0.8 0.5  0.5    %8 nominal-confidence fire
+    1.0 0.6  0.6    %9 high-confidence fire
+];
+% cmap= zeros(10,3);
+% cmap(8:10,:)=1;
+if any(data(:)<0 | data(:)>9), 
+    warning('Value out of range 0 to 9 for MODIS14 data')
+end
+
+[rows,cols]=size(data);
+fprintf('rows                  %i\n',rows)
+fprintf('cols                  %i\n',cols)
+fprintf('top left X            %19.15f\n',geo(1))
+fprintf('W-E pixel resolution  %19.15f\n',geo(2))
+fprintf('rotation, 0=North up  %19.15f\n',geo(3))
+fprintf('top left Y            %19.15f\n',geo(4))
+fprintf('rotation, 0=North up  %19.15f\n',geo(5))
+fprintf('N-S pixel resolution  %19.15f\n',geo(6))
+if geo(3)~=0 | geo(5)~=0,
+    error('rotation not supported')
+end
+for i=1:10,
+    count(i+1)=sum(data(:)==i);
+end
+unknown= count(1)+count(2)+count(3)+count(7);
+water  = count(4);
+cloud  = count(5);
+land   = count(6);
+fire   = count(8:10);
+fprintf('unprocessed/unknown   %i\n',unknown)
+fprintf('water                 %i\n',water)
+fprintf('land                  %i\n',land)
+fprintf('cloud                 %i\n',cloud)
+fprintf('low-confidence fire   %i\n',fire(1))
+fprintf('nominal-confid fire   %i\n',fire(2))
+fprintf('high-confidence fire  %i\n',fire(3))
+
+lon = geo(1)+[0:rows-1]*geo(2);
+lat = geo(4)+[cols-1:-1:0]*geo(6);
+image(lon,lat,data);
+colormap(cmap);
+xlabel('Longitude (deg)')
+ylabel('Latitude (deg)')
+if ~exist('s','var'), s='MOD14'; end
+t = sprintf('%s fire pixels %i %i %i',s,fire);
+title(t,'Interpreter','none')
+grid on
+% show fire pixels
+% hold on, [i,j,c]=find(data > 6); plot(lat(j),lon(i),'+r'), hold off
+end
+
diff --git a/other/Matlab/detection/showmod14files.m b/other/Matlab/detection/showmod14files.m
new file mode 100644
index 00000000..c24a099e
--- /dev/null
+++ b/other/Matlab/detection/showmod14files.m
@@ -0,0 +1,16 @@
+function showmod14files(search)
+% display one or more MOD14 Matlab files
+% input:
+%   search  file name, or search string, default '*.mat'
+if ~exist('search','var'),
+    search='*.mat';
+end
+files=dir(search);
+for i=1:length(files),
+    file=files(i).name;
+    load(file)
+    s=file;
+    showmod14(data,geotransform,s)
+    pause(2)
+end
+    
-- 
2.11.4.GIT