Commit 790ae792 authored by luroth's avatar luroth
Browse files

.

parent 284f8315
......@@ -102,7 +102,7 @@ def create_GeoTiff(file_path, array, data_type):
driver = gdal.GetDriverByName('GTiff')
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.WriteArray(array)
......
......@@ -21,6 +21,7 @@ import csv
from scipy import ndimage
from Common import create_GeoTiff
import gdal
import colorsys
def main(argv):
camera_position_file = ''
......@@ -69,15 +70,17 @@ def main(argv):
image[image_mask_all==0,:] = 0
# Colors
R = image[:,0]/255.0
G = image[:,1]/255.0
B = image[:,2]/255.0
# ExG
R = image[:,0]
G = image[:,1]
B = image[:,2]
r= np.divide(R, (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
ExG = 2.0* g - r - b
ExR = 1.4 * r - b
misc.imsave(os.path.join(output_directory, PhotoID + '_ExG.tif'), ExG.reshape(origin_shape[0], origin_shape[1]))
......@@ -85,10 +88,11 @@ def main(argv):
ExG_ExR = ExG - ExR
plants = np.array(ExG_ExR>0)
plants = np.array(ExG_ExR, dtype=float)
plants = plants.reshape(origin_shape[0], origin_shape[1])
#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))
create_GeoTiff(os.path.join(output_directory, PhotoID + '_segmented_ExG_ExR'), plants, gdal.GDT_Float32)
#ExG_soil = ExG[image_mask_soil==1]
#print(ExG_soil.size)
......
......@@ -12,6 +12,7 @@ import numpy as np
from scipy.cluster.vq import kmeans,vq
from scipy import misc
import sklearn.cluster
import multiprocessing
import matplotlib.pyplot as plt
from osgeo import gdalnumeric
import os
......@@ -24,6 +25,86 @@ import json
from Common import *
from osgeo import ogr, gdal
import pandas as pd
from multiprocessing import Process, Queue, Manager
import time
from progressbar import ProgressBar
def zonal_worker(work_queue, result_queue, metadata_folder, sub_plots):
for job in iter(work_queue.get, 'STOP'):
PhotoID = job['PhotoID']
result = zonal_stat(metadata_folder, sub_plots, PhotoID)
result_queue.put(result)
def zonal_stat(metadata_folder, sub_plots, PhotoID):
samples = []
subplots_ds = gdal.Open(os.path.join(metadata_folder, PhotoID + '_mask_subsamples.tif'))
subplots = np.array(subplots_ds.GetRasterBand(1).ReadAsArray(), dtype=int)
origin_shape = subplots.shape
subplots = subplots.reshape(-1)
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 += DELTA_x
photo_meta_y_ds = gdal.Open(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_xy = np.stack((photo_meta_x, photo_meta_y), axis=2)
photo_meta_xy = photo_meta_xy.reshape(-1,2)
photo_meta_x = photo_meta_x.reshape(-1)
photo_meta_y = photo_meta_y.reshape(-1)
segmented = misc.imread(os.path.join(metadata_folder, PhotoID + '_segmented_ExG.tif'))
segmented = segmented.reshape(-1)
photo_meta_z = misc.imread(os.path.join(metadata_folder, PhotoID + '_metadata_coord_z.tif'))
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 = 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 = photo_meta_angle_xy.reshape(-1)
for sub_plot in sub_plots:
plot_id = sub_plot['properties']['plot_number']
subsample = sub_plot['properties']['subsample']
#print("calculating plot id", plot_id, ', subsample', subsample)
mask_subplot = subplots == (plot_id*100+subsample)
subplot_segmented = segmented[mask_subplot]
if subplot_segmented.size > 0:
subplot_mean_cc = np.mean(subplot_segmented)
subplot_mean_x = np.mean(photo_meta_x[mask_subplot])
subplot_mean_y = np.mean(photo_meta_y[mask_subplot])
subplot_mean_z = np.mean(photo_meta_z[mask_subplot])
subplot_mean_angle_z = np.mean(photo_meta_angle_z[mask_subplot])
subplot_mean_angle_xy = np.mean(photo_meta_angle_xy[mask_subplot])
image_x, image_y = np.where(mask_subplot.reshape(origin_shape))
image_x = np.mean(image_x)
image_y = np.mean(image_y)
#print(subplot_mean_cc, end="", flush=True)
samples.append({"plot_id":plot_id, 'subsample':subsample, 'CC':subplot_mean_cc,
'glob_X':subplot_mean_x,
'glob_Y':subplot_mean_y,
'glob_Z':subplot_mean_z,
'angle_z':subplot_mean_angle_z,
'angle_xy':subplot_mean_angle_xy,
'image':PhotoID,
'image_x':image_x,
'image_y':image_y})
else:
pass
return samples
def main(argv):
camera_position_file = ''
......@@ -57,82 +138,57 @@ def main(argv):
with open(sub_plots_file, 'r') as f:
sub_plots_tmp = geojson.load(f)
sub_plots = sub_plots_tmp['features']
with open(camera_position_file, 'r') as csvfile:
spamreader = csv.reader(csvfile, delimiter='\t', quotechar='|')
next(spamreader, None)
next(spamreader, None)
m = Manager()
jobs = m.Queue()
results = m.Queue()
processes = []
max_jobs = 0
count = 0
for row in spamreader:
PhotoID = os.path.splitext(row[0])[0]
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'))):
subplots_ds = gdal.Open(os.path.join(metadata_folder, PhotoID + '_mask_subsamples.tif'))
subplots = np.array(subplots_ds.GetRasterBand(1).ReadAsArray(), dtype=int)
subplots = subplots.reshape(-1)
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 += DELTA_x
photo_meta_y_ds = gdal.Open(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
job = dict()
job['PhotoID'] = PhotoID
jobs.put(job)
max_jobs += 1
photo_meta_xy = np.stack((photo_meta_x, photo_meta_y), axis=2)
photo_meta_xy = photo_meta_xy.reshape(-1,2)
photo_meta_x = photo_meta_x.reshape(-1)
photo_meta_y = photo_meta_y.reshape(-1)
segmented = misc.imread(os.path.join(metadata_folder, PhotoID + '_segmented_ExG.tif'))
segmented = segmented.reshape(-1)
photo_meta_z = misc.imread(os.path.join(metadata_folder, PhotoID + '_metadata_coord_z.tif'))
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 = 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 = photo_meta_angle_xy.reshape(-1)
progress = ProgressBar(min_value=0, max_value=max_jobs)
for sub_plot in sub_plots:
plot_id = sub_plot['properties']['plot_number']
subsample = sub_plot['properties']['subsample']
#print("calculating plot id", plot_id, ', subsample', subsample)
mask_subplot = subplots == (plot_id*100+subsample)
subplot_segmented = segmented[mask_subplot]
if subplot_segmented.size > 0:
subplot_mean_cc = np.mean(subplot_segmented)
subplot_mean_x = np.mean(photo_meta_x[mask_subplot])
subplot_mean_y = np.mean(photo_meta_y[mask_subplot])
subplot_mean_z = np.mean(photo_meta_z[mask_subplot])
subplot_mean_angle_z = np.mean(photo_meta_angle_z[mask_subplot])
subplot_mean_angle_xy = np.mean(photo_meta_angle_xy[mask_subplot])
print(subplot_mean_cc, end="", flush=True)
all_cc.append({"plot_id":plot_id, 'subsample':subsample, 'CC':subplot_mean_cc,
'glob_X':subplot_mean_x,
'glob_Y':subplot_mean_y,
'glob_Z':subplot_mean_z,
'angle_z':subplot_mean_angle_z,
'angle_xy':subplot_mean_angle_xy,
'image':PhotoID})
else:
print("-", end="", flush=True)
with open(os.path.join(metadata_folder, "CC.json"), 'w') as f:
json.dump(all_cc, f)
for w in range(multiprocessing.cpu_count()-2):
p = Process(target=zonal_worker, args=(jobs, results, metadata_folder, sub_plots))
p.daemon = True
p.start()
processes.append(p)
jobs.put('STOP')
print("jobs all started")
samples = []
progress.update(0)
while count < (max_jobs):
samples_ = results.get()
samples.extend(samples_)
progress.update(count)
count+=1
progress.finish()
for p in processes:
p.join()
pd_data = pd.read_json(json.dump(all_cc), orient='records')
pd_data.to_csv(os.path.join(metadata_folder, "CC.vsc"))
pd_data = pd.DataFrame(samples)
pd_data.to_csv(os.path.join(metadata_folder, "CC.csv"))
if __name__ == "__main__":
main(sys.argv[1:])
......
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