Commit 0b29cccc authored by luroth's avatar luroth
Browse files

.

parent 284f8315
...@@ -97,17 +97,23 @@ def point_in_poly(x,y,poly): ...@@ -97,17 +97,23 @@ def point_in_poly(x,y,poly):
return inside return inside
def create_GeoTiff(file_path, array, data_type): def create_1band_GeoTiff(file_path, array, data_type):
driver = gdal.GetDriverByName('GTiff') driver = gdal.GetDriverByName('GTiff')
file = file_path + '.tif' file = file_path + '.tif'
ds = driver.Create( file, array.shape[0], array.shape[1], 1, data_type) ds = driver.Create( file, array.shape[1], array.shape[0], 1, data_type)
band_1 = ds.GetRasterBand(1) band_1 = ds.GetRasterBand(1)
band_1.WriteArray(array) band_1.WriteArray(array)
band_1.SetNoDataValue(0) if data_type is gdal.GDT_Byte or data_type is gdal.GDT_Int16:
band_1.SetNoDataValue(0)
ds = None ds = None
return file_path return file_path
def read_1band_file(file_path):
src_ds = gdal.Open(file_path)
band = src_ds.GetRasterBand(1).ReadAsArray()
return band
...@@ -25,13 +25,13 @@ from progressbar import ProgressBar ...@@ -25,13 +25,13 @@ from progressbar import ProgressBar
from Common import * from Common import *
def mask_worker(work_queue, result_queue, metadata_folder, output_directory, sub_plots): def mask_worker(work_queue, result_queue, metadata_folder, output_directory, lot, sub_plots):
for job in iter(work_queue.get, 'STOP'): for job in iter(work_queue.get, 'STOP'):
PhotoID = job['PhotoID'] PhotoID = job['PhotoID']
create_mask(metadata_folder, output_directory, sub_plots, PhotoID) create_mask(metadata_folder, output_directory, lot, sub_plots, PhotoID)
result_queue.put(PhotoID) result_queue.put(PhotoID)
def create_mask(metadata_folder, output_directory, sub_plots, PhotoID): def create_mask(metadata_folder, output_directory, lot, sub_plots, PhotoID):
photo_meta_x_ds = gdal.Open(os.path.join(metadata_folder, PhotoID + '_metadata_coord_x.tif')) photo_meta_x_ds = gdal.Open(os.path.join(metadata_folder, PhotoID + '_metadata_coord_x.tif'))
photo_meta_x = np.array(photo_meta_x_ds.GetRasterBand(1).ReadAsArray(), dtype=np.float64) photo_meta_x = np.array(photo_meta_x_ds.GetRasterBand(1).ReadAsArray(), dtype=np.float64)
...@@ -47,6 +47,15 @@ def create_mask(metadata_folder, output_directory, sub_plots, PhotoID): ...@@ -47,6 +47,15 @@ def create_mask(metadata_folder, output_directory, sub_plots, PhotoID):
photo_meta_xy = np.stack((photo_meta_x, photo_meta_y), axis=2) photo_meta_xy = np.stack((photo_meta_x, photo_meta_y), axis=2)
photo_meta_xy = photo_meta_xy.reshape(-1,2) photo_meta_xy = photo_meta_xy.reshape(-1,2)
lot_path = path.Path(lot)
mask_lot = lot_path.contains_points(photo_meta_xy, radius=0.02)
if not np.any(mask_lot):
return
mask_lot = mask_lot.reshape(photo_meta_y.shape)
#gdalnumeric.SaveArray(1*mask_lot, os.path.join(output_directory, PhotoID + '_mask_lot.tif'), format='GTiff')
create_1band_GeoTiff(os.path.join(output_directory, PhotoID + '_mask_lot.tif'), 1*mask_lot, gdal.GDT_Byte)
mask_subplots = np.empty(photo_meta_xy.shape[0:1], dtype=int) mask_subplots = np.empty(photo_meta_xy.shape[0:1], dtype=int)
...@@ -67,8 +76,8 @@ def create_mask(metadata_folder, output_directory, sub_plots, PhotoID): ...@@ -67,8 +76,8 @@ def create_mask(metadata_folder, output_directory, sub_plots, PhotoID):
if found: if found:
mask_subplots = mask_subplots.reshape(photo_meta_y.shape) mask_subplots = mask_subplots.reshape(photo_meta_y.shape)
create_1band_GeoTiff(os.path.join(output_directory, PhotoID + '_mask_subsamples.tif'), mask_subplots, gdal.GDT_Int16)
gdalnumeric.SaveArray(mask_subplots, os.path.join(output_directory, PhotoID + '_mask_subsamples.tif'), format='GTiff') #gdalnumeric.SaveArray(mask_subplots, os.path.join(output_directory, PhotoID + '_mask_subsamples.tif'), format='GTiff')
def main(argv): def main(argv):
camera_position_file = '' camera_position_file = ''
...@@ -96,6 +105,10 @@ def main(argv): ...@@ -96,6 +105,10 @@ def main(argv):
sub_plots_file = arg sub_plots_file = arg
print('Metdadata folder is', metadata_folder) print('Metdadata folder is', metadata_folder)
with open(lot_dim, 'r') as f:
lot_polygons = geojson.load(f)
lot = lot_polygons['features'][0]['geometry']['coordinates'][0]
with open(sub_plots_file, 'r') as f: with open(sub_plots_file, 'r') as f:
sub_plots_tmp = geojson.load(f) sub_plots_tmp = geojson.load(f)
...@@ -126,7 +139,7 @@ def main(argv): ...@@ -126,7 +139,7 @@ def main(argv):
progress = ProgressBar(min_value=0, max_value=max_jobs) progress = ProgressBar(min_value=0, max_value=max_jobs)
for w in range(multiprocessing.cpu_count()-2): for w in range(multiprocessing.cpu_count()-2):
p = Process(target=mask_worker, args=(jobs, results, metadata_folder, output_directory, sub_plots)) p = Process(target=mask_worker, args=(jobs, results, metadata_folder, output_directory, lot, sub_plots))
p.daemon = True p.daemon = True
p.start() p.start()
processes.append(p) processes.append(p)
......
...@@ -5,21 +5,20 @@ Created on 16.03.2017 ...@@ -5,21 +5,20 @@ Created on 16.03.2017
''' '''
#from pylab import plot,show #from pylab import plot,show
import colorsys
from numpy import vstack,array from numpy import vstack,array
from numpy.random import rand from numpy.random import rand
import numpy.ma as ma import numpy.ma as ma
import numpy as np import numpy as np
from scipy.cluster.vq import kmeans,vq from scipy.cluster.vq import kmeans,vq
from scipy import misc from scipy import misc
#import sklearn.cluster
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from osgeo import gdalnumeric
import os import os
import sys import sys
import getopt import getopt
import csv import csv
from scipy import ndimage from scipy import ndimage
from Common import create_GeoTiff from Common import create_1band_GeoTiff, read_1band_file
import gdal import gdal
def main(argv): def main(argv):
...@@ -54,56 +53,39 @@ def main(argv): ...@@ -54,56 +53,39 @@ def main(argv):
for row in spamreader: for row in spamreader:
PhotoID = os.path.splitext(row[0])[0] PhotoID = os.path.splitext(row[0])[0]
image = misc.imread(os.path.join(metadata_folder, '../../RGB', PhotoID + '.JPG')) src_ds = gdal.Open(os.path.join(metadata_folder, '../../RGB', PhotoID + '.JPG'))
R = src_ds.GetRasterBand(1).ReadAsArray()
G = src_ds.GetRasterBand(2).ReadAsArray()
B = src_ds.GetRasterBand(3).ReadAsArray()
if os.path.isfile(os.path.join(metadata_folder, PhotoID + '_mask_lot.tif')): if os.path.isfile(os.path.join(metadata_folder, PhotoID + '_mask_lot.tif')):
image_mask_all = misc.imread(os.path.join(metadata_folder, PhotoID + '_mask_lot.tif')) image_mask_all = read_1band_file(os.path.join(metadata_folder, PhotoID + '_mask_lot.tif'))
image_mask_soil = misc.imread(os.path.join(metadata_folder, PhotoID + '_mask_soil.tif'))
origin_shape = image.shape
image = image.reshape(-1,3)
image_mask_all = image_mask_all.reshape(-1) image_mask_all = image_mask_all.reshape(-1)
image_mask_soil = image_mask_soil.reshape(-1) R[image_mask_all==0] = 0
G[image_mask_all==0] = 0
image[image_mask_all==0,:] = 0 B[image_mask_all==0] = 0
# ExG origin_shape = R.shape
R = image[:,0]
G = image[:,1] r= np.divide(R, (np.mean(R)+np.mean(G)+np.mean(B)), dtype=float)
B = image[:,2] g= np.divide(G, (np.mean(R)+np.mean(G)+np.mean(B)), dtype=float)
r= np.divide(R, (np.mean(R)+np.mean(G)+np.mean(B)), dtype=float) b= np.divide(B, (np.mean(R)+np.mean(G)+np.mean(B)), dtype=float)
g= np.divide(G, (np.mean(R)+np.mean(G)+np.mean(B)), dtype=float)
b= np.divide(B, (np.mean(R)+np.mean(G)+np.mean(B)), dtype=float) ExG = 2* g - r - b
ExR = 1.4 * r - g
ExG = 2* g - r - b
ExR = 1.4 * r - b create_1band_GeoTiff(os.path.join(output_directory, PhotoID + '_ExG.tif'), ExG.reshape(origin_shape), gdal.GDT_Float32)
create_1band_GeoTiff(os.path.join(output_directory, PhotoID + '_ExR.tif'), ExR.reshape(origin_shape), gdal.GDT_Float32)
misc.imsave(os.path.join(output_directory, PhotoID + '_ExG.tif'), ExG.reshape(origin_shape[0], origin_shape[1]))
misc.imsave(os.path.join(output_directory, PhotoID + '_ExR.tif'), ExR.reshape(origin_shape[0], origin_shape[1])) ExG_ExR = ExG - ExR
ExG_ExR = ExG - ExR plants = np.array(ExG_ExR>=0.0)
plants = np.array(ExG_ExR>0) create_GeoTiff(os.path.join(output_directory, PhotoID + '_segmented_ExG_ExR'), plants.reshape(origin_shape[0], origin_shape[1]), gdal.GDT_Byte)
#gdalnumeric.SaveArray(1*(plants.reshape(origin_shape[0], origin_shape[1])), os.path.join(output_directory, PhotoID + '_segmented_ExG_ExR.tif'), format='GTiff')
create_GeoTiff(os.path.join(output_directory, PhotoID + '_segmented_ExG_ExR', ExG_ExR, gdal.GDT_Byte)) print(PhotoID, "done")
#ExG_soil = ExG[image_mask_soil==1]
#print(ExG_soil.size)
#if ExG_soil.size ==0:
# continue
#ExG_threshold = np.percentile(ExG_soil, 99.95)
#soil = np.array(ExG > ExG_threshold)
#gdalnumeric.SaveArray(1*(soil.reshape(origin_shape[0], origin_shape[1])), os.path.join(output_directory, PhotoID + '_segmented_ExG.tif'), format='GTiff')
#misc.imsave(os.path.join(output_directory, PhotoID + '_segmented_ExG.tif'), soil.reshape(origin_shape[0], origin_shape[1]))
print(PhotoID, "done")
if __name__ == "__main__": if __name__ == "__main__":
main(sys.argv[1:]) main(sys.argv[1:])
......
...@@ -11,7 +11,6 @@ import numpy.ma as ma ...@@ -11,7 +11,6 @@ import numpy.ma as ma
import numpy as np import numpy as np
from scipy.cluster.vq import kmeans,vq from scipy.cluster.vq import kmeans,vq
from scipy import misc from scipy import misc
import sklearn.cluster
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from osgeo import gdalnumeric from osgeo import gdalnumeric
import os import os
...@@ -68,18 +67,15 @@ def main(argv): ...@@ -68,18 +67,15 @@ def main(argv):
print("\n", PhotoID) print("\n", PhotoID)
if (os.path.isfile(os.path.join(metadata_folder, PhotoID + '_mask_subsamples.tif')) and if (os.path.isfile(os.path.join(metadata_folder, PhotoID + '_mask_subsamples.tif')) and
os.path.isfile(os.path.join(metadata_folder, PhotoID + '_segmented_ExG.tif'))): os.path.isfile(os.path.join(metadata_folder, PhotoID + '_segmented_ExG_ExR.tif'))):
subplots_ds = gdal.Open(os.path.join(metadata_folder, PhotoID + '_mask_subsamples.tif')) subplots = read_1band_file(os.path.join(metadata_folder, PhotoID + '_mask_subsamples.tif'))
subplots = np.array(subplots_ds.GetRasterBand(1).ReadAsArray(), dtype=int)
subplots = subplots.reshape(-1) subplots = subplots.reshape(-1)
photo_meta_x_ds = gdal.Open(os.path.join(metadata_folder, PhotoID + '_metadata_coord_x.tif')) photo_meta_x = read_1band_file(os.path.join(metadata_folder, PhotoID + '_metadata_coord_x.tif'))
photo_meta_x = np.array(photo_meta_x_ds.GetRasterBand(1).ReadAsArray(), dtype=np.float64)
photo_meta_x += DELTA_x photo_meta_x += DELTA_x
photo_meta_y_ds = gdal.Open(os.path.join(metadata_folder, PhotoID + '_metadata_coord_y.tif')) photo_meta_y = read_1band_file(os.path.join(metadata_folder, PhotoID + '_metadata_coord_y.tif'))
photo_meta_y = np.array(photo_meta_y_ds.GetRasterBand(1).ReadAsArray(), dtype=np.float64)
photo_meta_y += DELTA_y photo_meta_y += DELTA_y
photo_meta_xy = np.stack((photo_meta_x, photo_meta_y), axis=2) photo_meta_xy = np.stack((photo_meta_x, photo_meta_y), axis=2)
...@@ -88,15 +84,15 @@ def main(argv): ...@@ -88,15 +84,15 @@ def main(argv):
photo_meta_x = photo_meta_x.reshape(-1) photo_meta_x = photo_meta_x.reshape(-1)
photo_meta_y = photo_meta_y.reshape(-1) photo_meta_y = photo_meta_y.reshape(-1)
segmented = misc.imread(os.path.join(metadata_folder, PhotoID + '_segmented_ExG.tif')) segmented = read_1band_file(os.path.join(metadata_folder, PhotoID + '_segmented_ExG_ExR.tif'))
segmented = segmented.reshape(-1) segmented = segmented.reshape(-1)
photo_meta_z = misc.imread(os.path.join(metadata_folder, PhotoID + '_metadata_coord_z.tif')) photo_meta_z = read_1band_file(os.path.join(metadata_folder, PhotoID + '_metadata_coord_z.tif'))
photo_meta_z = photo_meta_z.reshape(-1) photo_meta_z = photo_meta_z.reshape(-1)
photo_meta_angle_z = misc.imread(os.path.join(metadata_folder, PhotoID + '_metadata_angle_z.tif')) photo_meta_angle_z = read_1band_file(os.path.join(metadata_folder, PhotoID + '_metadata_angle_z.tif'))
photo_meta_angle_z = photo_meta_angle_z.reshape(-1) photo_meta_angle_z = photo_meta_angle_z.reshape(-1)
photo_meta_angle_xy = misc.imread(os.path.join(metadata_folder, PhotoID + '_metadata_angle_xy.tif')) photo_meta_angle_xy = read_1band_file(os.path.join(metadata_folder, PhotoID + '_metadata_angle_xy.tif'))
photo_meta_angle_xy = photo_meta_angle_xy.reshape(-1) photo_meta_angle_xy = photo_meta_angle_xy.reshape(-1)
for sub_plot in sub_plots: for sub_plot in sub_plots:
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment