To receive notifications about scheduled maintenance, please subscribe to the mailing-list gitlab-operations@sympa.ethz.ch. You can subscribe to the mailing-list at https://sympa.ethz.ch

Commit acaafbd7 authored by luroth's avatar luroth
Browse files

.

parent 95daf32d
......@@ -97,17 +97,23 @@ def point_in_poly(x,y,poly):
return inside
def create_GeoTiff(file_path, array, data_type):
def create_1band_GeoTiff(file_path, array, data_type):
driver = gdal.GetDriverByName('GTiff')
file = file_path + '.tif'
ds = driver.Create(file, array.shape[1], array.shape[0], 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)
band_1.SetNoDataValue(0)
if data_type is gdal.GDT_Byte:
band_1.SetNoDataValue(0)
ds = None
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
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'):
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)
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 = 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):
photo_meta_xy = np.stack((photo_meta_x, photo_meta_y), axis=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_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)
......@@ -67,8 +76,8 @@ def create_mask(metadata_folder, output_directory, sub_plots, PhotoID):
if found:
mask_subplots = mask_subplots.reshape(photo_meta_y.shape)
gdalnumeric.SaveArray(mask_subplots, os.path.join(output_directory, PhotoID + '_mask_subsamples.tif'), format='GTiff')
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')
def main(argv):
camera_position_file = ''
......@@ -96,6 +105,10 @@ def main(argv):
sub_plots_file = arg
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:
sub_plots_tmp = geojson.load(f)
......@@ -126,7 +139,7 @@ def main(argv):
progress = ProgressBar(min_value=0, max_value=max_jobs)
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.start()
processes.append(p)
......
......@@ -5,23 +5,21 @@ Created on 16.03.2017
'''
#from pylab import plot,show
import colorsys
from numpy import vstack,array
from numpy.random import rand
import numpy.ma as ma
import numpy as np
from scipy.cluster.vq import kmeans,vq
from scipy import misc
#import sklearn.cluster
import matplotlib.pyplot as plt
from osgeo import gdalnumeric
import os
import sys
import getopt
import csv
from scipy import ndimage
from Common import create_GeoTiff
from Common import create_GeoTiff, read_1band_file
import gdal
import colorsys
def main(argv):
camera_position_file = ''
......@@ -55,59 +53,40 @@ def main(argv):
for row in spamreader:
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')):
image_mask_all = misc.imread(os.path.join(metadata_folder, PhotoID + '_mask_lot.tif'))
image_mask_soil = misc.imread(os.path.join(metadata_folder, PhotoID + '_mask_soil.tif'))
image_mask_all = read_1band_file(os.path.join(metadata_folder, PhotoID + '_mask_lot.tif'))
origin_shape = image.shape
image = image.reshape(-1,3)
image_mask_all = image_mask_all.reshape(-1)
image_mask_soil = image_mask_soil.reshape(-1)
image[image_mask_all==0,:] = 0
# Colors
R = image[:,0]/255.0
G = image[:,1]/255.0
B = image[:,2]/255.0
origin_shape = R.shape
# ExG
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.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]))
misc.imsave(os.path.join(output_directory, PhotoID + '_ExR.tif'), ExR.reshape(origin_shape[0], origin_shape[1]))
ExG_ExR = ExG - ExR
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'), plants, gdal.GDT_Float32)
#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")
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
ExR = 1.4 * r - g
create_GeoTiff(os.path.join(output_directory, PhotoID + '_ExG.tif'), ExG.reshape(origin_shape), gdal.GDT_Float32)
create_GeoTiff(os.path.join(output_directory, PhotoID + '_ExR.tif'), ExR.reshape(origin_shape), gdal.GDT_Float32)
ExG_ExR = ExG - ExR
plants = np.array(ExG_ExR>=0.0)
create_GeoTiff(os.path.join(output_directory, PhotoID + '_segmented_ExG_ExR'), plants.reshape(origin_shape[0], origin_shape[1]), gdal.GDT_Byte)
print(PhotoID, "done")
if __name__ == "__main__":
main(sys.argv[1:])
......
......@@ -11,8 +11,6 @@ import numpy.ma as ma
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
......@@ -20,91 +18,10 @@ import sys
import getopt
import csv
from scipy import ndimage
import geojson
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 = ''
......@@ -138,57 +55,79 @@ 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
os.path.isfile(os.path.join(metadata_folder, PhotoID + '_segmented_ExG.tif'))):
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_ExR.tif'))):
job = dict()
job['PhotoID'] = PhotoID
jobs.put(job)
max_jobs += 1
progress = ProgressBar(min_value=0, max_value=max_jobs)
subplots = read_1band_file(os.path.join(metadata_folder, PhotoID + '_mask_subsamples.tif'))
subplots = subplots.reshape(-1)
photo_meta_x = read_1band_file(os.path.join(metadata_folder, PhotoID + '_metadata_coord_x.tif'))
photo_meta_x += DELTA_x
photo_meta_y = read_1band_file(os.path.join(metadata_folder, PhotoID + '_metadata_coord_y.tif'))
photo_meta_y += DELTA_y
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
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 = read_1band_file(os.path.join(metadata_folder, PhotoID + '_segmented_ExG_ExR.tif'))
segmented = segmented.reshape(-1)
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_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_xy = read_1band_file(os.path.join(metadata_folder, PhotoID + '_metadata_angle_xy.tif'))
photo_meta_angle_xy = photo_meta_angle_xy.reshape(-1)
progress.finish()
for p in processes:
p.join()
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)
pd_data = pd.DataFrame(samples)
pd_data.to_csv(os.path.join(metadata_folder, "CC.csv"))
pd_data = pd.read_json(json.dump(all_cc), orient='records')
pd_data.to_csv(os.path.join(metadata_folder, "CC.vsc"))
if __name__ == "__main__":
main(sys.argv[1:])
......
Markdown is supported
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