diff --git a/tlib/tgeo/__init__.py b/tlib/tgeo/__init__.py index f6965fcd2ec75c62907829caec52bb38f04ac8bf..915e56fbcf0e2faa7088e2791242a3974e58d810 100644 --- a/tlib/tgeo/__init__.py +++ b/tlib/tgeo/__init__.py @@ -1,2 +1,2 @@ from distutils.command.config import config -from . import geotif, kml, landcover, plotter, tiling, utils +from . import geotif, kml, landcover, s2, tiling, utils diff --git a/tlib/tgeo/plotter.py b/tlib/tgeo/plotter.py deleted file mode 100644 index 4de29662c8b75878deec7e888114efe421b63cb7..0000000000000000000000000000000000000000 --- a/tlib/tgeo/plotter.py +++ /dev/null @@ -1,712 +0,0 @@ -import os - -import numpy as np -import matplotlib.pyplot as plt -from tqdm.auto import tqdm - -from . import geotif, landcover, utils -from . import kml as tgeo_kml -from tlib import tutils - - -class Plotter: - """ - This class helps to plot images, that were exported using the tgeo.gee.s2.S2Exporter() class. - However, you can also plot other Sentinel-2 images, that were not exported with this class. - - :param folder_s2: folder which contains Sentinel-2 geotifs; if given, it is enough to use file names for plotting; - if None, you need the whole file paths - :param folder_scl: folder which contains geotifs with the SCL band of Sentinel-2; if None, no SCL can be plotted - :param folder_lcs: folder which contains geotifs with landcover dataset; if None, no landcover data can be plotted - :param folder_masks: folder which contains geotifs with masks; if None, no masks can be plotted - :param val_range_s2: value range that is plotted if rgb images are plotted of Sentinel-2 data - :param channel_indices: if your geotif channels have names (e.g. 'B2', 'B3') etc., you do not need this; - if your geotif channels do not have names, you need to define a dictionary here indicating the indices - of the channels, e.g. {'B2': 0, 'B3: 1, ...} - - Note: You can replace the following function variables to define more compplicated directories: - self.get_s2_dir - self.get_scl_dir - self.get_lcs_dir - self.get_masks_dir - """ - - fig_height = 4.6 # for figures plotted with matplotlib - output_folder = None # is determined in to_tif() function - - def __init__(self, folder, channel_indices: dict = None): - - self.folder = folder - self.channel_indices = channel_indices - - def get_file_dirs(self, *files): - """ - Return file paths according to given self folders for s2, scl, lcs and masks. - - :param files: list of files; if only file name (without directory) is given, paths for all self folders are - returned, if self folders are given; else, only path to files_s2 is returned and the others are None - :return: (list of files_s2, list of files_scl, list of files_lcs, list of files_masks) - """ - - # if only file name and not directory ... add folders to file name - if files[0].split('/')[0] != self.folder.split('/')[0] or files[0][0] == '/': - files = [os.path.join(self.folder, file) for file in files] if self.folder is not None else None - - else: # assume that given file is s2 file and set others None - files = files - - return files - - def plot(self): - pass - - def plot_on_map(self): - pass - - def to_tif(self): - pass - - def tif_to_kml(self, main_folder=None, markings='s2', save=True, kml=None): - - if main_folder == None: - main_folder = self.output_folder - - print(main_folder) - kml = tgeo_kml.subfolders_to_kml( - main_folder=main_folder, - overlay_folder='legends', - points_for=markings, - polygons_for=markings, - save=save, - kml=kml - ) - - return kml - - -class Plotters(Plotter): - - def __init__(self, plotters: list): - - self.plotters = plotters - - def plot(self, file, **kwargs): - - for plotter in self.plotters: - these_kwargs = {key: value for (key, value) in kwargs.items() if key in plotter.plot.__code__.co_varnames} - plotter.plot(file, **these_kwargs) - - def plot_on_map(self, *files, m=None, **kwargs): - - for plotter in self.plotters: - these_kwargs = { - key: value for (key, value) in kwargs.items() - if key in plotter.plot_on_map.__code__.co_varnames or key in utils.get_map.__code__.co_varnames - } - m = plotter.plot_on_map(*files, m=m, **these_kwargs) - - return m - - def to_tif(self, *files, output_folder, **kwargs): - - self.output_folder = output_folder - - for plotter in self.plotters: - these_kwargs = {key: value for (key, value) in kwargs.items() if key in plotter.to_tif.__code__.co_varnames} - plotter.to_tif(*files, output_folder=output_folder, **these_kwargs) - - -class LCPlotter(Plotter): - - def __init__(self, folder, channel_indices: dict = None): - - self.folder = folder - self.channel_indices = channel_indices - - - def plot(self, file, plot_legends: bool = False): - """ - Plots using matplotlib and plt.show(). No figure or axis is returned. - - :param file: file name or file dir - :return: - """ - - print(f'file: {file}') - file = self.get_file_dirs(file)[0] - - # plot lcs - fig, axs = plt.subplots(1, 4, figsize=(self.fig_height * 4, self.fig_height)) - axs[0] = geotif.plot_landcover(file=file, channel='CORINE', plot_legend=plot_legends, ax=axs[0]) - axs[1] = geotif.plot_landcover(file=file, channel='MODIS_1', plot_legend=plot_legends, ax=axs[1]) - axs[2] = geotif.plot_landcover(file=file, channel='CGLS', plot_legend=plot_legends, ax=axs[2]) - axs[3] = geotif.plot_landcover(file=file, channel='GlobCover', plot_legend=plot_legends, ax=axs[3]) - - axs[0].set_title('CORINE') - axs[1].set_title('MODIS_1') - axs[2].set_title('CGLS') - axs[3].set_title('GlobCover') - fig.tight_layout() - plt.show() - - def plot_on_map( - self, - *files, - - plot_lcs: list = None, # CORINE, MODIS_1, CGLS, GlobCover - plot_legends: bool = False, - plot_locations: bool = False, - - m=None, - **map_kwargs, - ): - """ - Plots on ipyleaflet map. - - :param files: list of file names or file dirs - :param lcs: list of land cover data that shall be plotted, e.g. ['CORINE', 'MODIS_1', 'CGLS', 'GlobCover'] - :param plot_legend: bool if legends of land cover data shall be plotted with matplotlib - :param plot_locations: bool if polygons shall be plotted - :param m: ipyleaflet map; if None a new one is created - :param map_kwargs: ipyleaflet map kwargs - :return: - """ - - files = self.get_file_dirs(*files) - - if m is None: - m = utils.get_map(**map_kwargs) - - # plot locations - if plot_locations: - m = geotif.location_on_map(*files, m=m) - - # plot landcovers - if plot_lcs is not None: - for lc in plot_lcs: - m = geotif.landcover_on_map(*files, channel=lc, plot_legend=plot_legends, m=m) - - return m - - def to_tif(self, *files, output_folder, plot_lcs: list = None): - """ - Writes readable tif files (rgb or grey-channel) into given output_folder. - - :param files: list of file names or file dirs - :param output_folder: output folder in which tifs are stored - :param lcs: list of land cover data that shall be plotted, e.g. ['CORINE', 'MODIS_1', 'CGLS', 'GlobCover'] - :return: - """ - - files = self.get_file_dirs(*files) - self.output_folder = output_folder - - # plot lcs - if plot_lcs is not None: - for lc in plot_lcs: - this_output_folder = tutils.files.join_paths(output_folder, f'lcs_{lc}') - - geotif.landcover_to_tif(*files, output_folder=this_output_folder, channel=lc) - - # legend - fig = landcover.plot_legend(lc) - this_output_folder = tutils.files.join_paths(output_folder, 'legends') - fig.savefig(os.path.join(this_output_folder, f'legend_{lc}.png'), facecolor='white', bbox_inches='tight') - plt.close() - - -class MaskPlotter(Plotter): - - def __init__(self, folder, channel_indices: dict = None): - - self.folder = folder - self.channel_indices = channel_indices - - self.cmap_mask = 'Reds'# 'gray' - - def plot(self, file): - """ - Plots using matplotlib and plt.show(). No figure or axis is returned. - - :param file: file name or file dir - :return: - """ - - print(f'file: {file}') - file = self.get_file_dirs(file)[0] - - # plot masks - channel_names = geotif.get_channel_names(file) - if len(channel_names) == 1: # TODO: because of using plt.subplots, a single plot is not possible with this coding and the mask is doubled in this case - channel_names = [channel_names[0]] * 2 - - fig, axs = plt.subplots(1, len(channel_names), figsize=(self.fig_height * len(channel_names), self.fig_height)) - for i in range(len(channel_names)): - array = geotif.get_array(file, channels=channel_names[i], transpose=True) - array = array.astype(bool) - - axs[i].imshow(array, cmap=self.cmap_mask, vmin=0, vmax=1) - axs[i].set_xticks([]), axs[i].set_yticks([]) - axs[i].set_title(channel_names[i]) - fig.tight_layout() - plt.show() - - def plot_on_map( - self, - *files, - - plot_masks: list = None, # ROI, Valid_Area - plot_locations: bool = False, - - m=None, - **map_kwargs, - ): - """ - Plots on ipyleaflet map. - - :param masks: list of mask channels that shall be plotted, e.g. ['ROI', 'Valid_Area'] - :param plot_locations: bool if polygons shall be plotted - :param m: ipyleaflet map; if None a new one is created - :param map_kwargs: ipyleaflet map kwargs - :return: - """ - - files = self.get_file_dirs(*files) - - if m is None: - m = utils.get_map(**map_kwargs) - - # plot locations - if plot_locations: - m = geotif.location_on_map(*files, m=m) - - # plot masks - if plot_masks is not None: - for mask in plot_masks: - m = geotif.image_on_map(*files, channels=mask, val_range=(0, 1), m=m) - - return m - - def to_tif(self, *files, output_folder, plot_masks: list=None): - - """ - Writes readable tif files (rgb or grey-channel) into given output_folder. - - :param files: list of file names or file dirs - :param output_folder: output folder in which tifs are stored - :param masks: list of mask channels that shall be plotted, e.g. ['ROI', 'Valid_Area'] - :return: - """ - - files = self.get_file_dirs(*files) - self.output_folder = output_folder - - # plot masks - if plot_masks is not None and files is not None: - for mask in plot_masks: - this_output_folder = tutils.files.join_paths(output_folder, f'masks_{mask}') - - geotif.image_to_tif(*files, output_folder=this_output_folder, channels=mask, val_range=(0, 1)) - - -class S2Plotter(Plotter): - - def __init__(self, folder, channel_indices: dict = None, val_range_s2: tuple = (0, 2**10)): - - self.folder = folder - self.channel_indices = channel_indices - - self.val_range_s2 = val_range_s2 - self.cmap_ndvi = 'RdYlGn' - - def plot(self, file, plot_s2_channels: bool = False): - """ - Plots using matplotlib and plt.show(). No figure or axis is returned. - - :param file: file name or file dir - :param s2_channels: bool if all channels of Sentinel-2 image are plotted individually - :return: - """ - - print(f'file: {file}') - file = self.get_file_dirs(file)[0] - - # plot s2 - fig, axs = plt.subplots(1,4, figsize=(self.fig_height * 4, self.fig_height)) - - # rgb image - if self.channel_indices is None: - axs[0] = geotif.plot_image( - file=file, channels=['B4', 'B3', 'B2'], val_range=self.val_range_s2, ax=axs[0]) - else: - channels = [self.channel_indices['B4'], self.channel_indices['B3'], self.channel_indices['B2']] - axs[0] = geotif.plot_image(file=file, channels=channels, val_range=self.val_range_s2, ax=axs[0]) - - axs[0].set_title('RGB') - - # false color image - if self.channel_indices is None: - axs[1] = geotif.plot_image( - file=file, channels=['B8', 'B4', 'B3'], val_range=self.val_range_s2, ax=axs[1]) - else: - channels = [self.channel_indices['B8'], self.channel_indices['B4'], self.channel_indices['B3']] - axs[1] = geotif.plot_image(file=file, channels=channels, val_range=self.val_range_s2, ax=axs[1]) - - axs[1].set_title('NIR-R-G') - - # ndvi - if self.channel_indices is None: - image = geotif.get_array(file=file, channels=('B8', 'B4')) - else: - channels = [self.channel_indices['B8'], self.channel_indices['B4']] - image = geotif.get_array(file=file, channels=channels) - - ndvi = utils.get_ndvi(nir=image[0], red=image[1]) - - axs[2] = tutils.plots.plot_original_image_size( - array=ndvi, plot_cbar=True, cmap=self.cmap_ndvi, vmin=-1, vmax=1, fig=fig, ax=axs[2]) - axs[2].set_title('NDVI') - - # evi - if self.channel_indices is None: - image = geotif.get_array(file=file, channels=('B8', 'B4', 'B2')) - else: - channels = [self.channel_indices['B8'], self.channel_indices['B4'], self.channel_indices['B2']] - image = geotif.get_array(file=file, channels=channels) - - evi = utils.get_evi(nir=image[0], red=image[1], blue=image[2]) - - q = 0.02 - - axs[3] = tutils.plots.plot_original_image_size( - array=evi, - plot_cbar=True, - cmap=self.cmap_ndvi, - vmin=np.quantile(evi, q=q), - vmax=np.quantile(evi, q=1-q), - fig=fig, - ax=axs[3]) - axs[3].set_title('EVI') - - plt.show() - - # plot s2 channels - if plot_s2_channels: - geotif.plot_channels(file, ncols=5, figsize=(self.fig_height * 5, self.fig_height * 2)) - plt.show() - - def plot_on_map( - self, - *files, - - plot_rgb: bool=False, - plot_false_colors: bool=False, - plot_channels: list = None, - plot_ndvi: bool=False, - plot_evi: bool = False, - plot_locations: bool = False, - - m=None, - **map_kwargs, - ): - """ - Plots on ipyleaflet map. - - :param files: list of file names or file dirs - :param plot_rgb: bool if Sentinel-2 rgbs are plotted - :param plot_false_colors: bool if false color images of Sentinel-2 are plotted - :param plot_channels: list of channels that shall be plotted from Sentinel-2, e.g. ['B2', 'B3'] - :param plot_ndvi: bool if ndvi of Sentinel-2 images shall be plotted - :param plot_evi: bool if evi of Sentinel-2 images shall be plotted - :param plot_locations: bool if polygons shall be plotted - :param m: ipyleaflet map; if None a new one is created - :param map_kwargs: ipyleaflet map kwargs - :return: - """ - - files = self.get_file_dirs(*files) - - if m is None: - m = utils.get_map(**map_kwargs) - - # plot locations - if plot_locations: - m = geotif.location_on_map(*files, m=m) - - # plot rgb - if plot_rgb: - if self.channel_indices is None: - m = geotif.image_on_map( - *files, channels=['B4', 'B3', 'B2'], val_range=self.val_range_s2, m=m) # RGB - else: - channels = [self.channel_indices['B4'], self.channel_indices['B3'], self.channel_indices['B2']] - m = geotif.image_on_map(*files, channels=channels, val_range=self.val_range_s2, m=m) # RGB - - # plot false_colors - if plot_false_colors: - if self.channel_indices is None: - m = geotif.image_on_map( - *files, channels=['B8', 'B4', 'B3'], val_range=self.val_range_s2, m=m) # NIR-R-G - else: - channels = [self.channel_indices['B8'], self.channel_indices['B4'], self.channel_indices['B3']] - m = geotif.image_on_map(*files, channels=channels, val_range=self.val_range_s2, m=m) # NIR-R-G - - # plot channels - if plot_channels is not None: - for channel in plot_channels: - m = geotif.image_on_map(*files, channels=channel, val_range='auto', m=m) - - # plot ndvi - if plot_ndvi: - if self.channel_indices is None: - images = [geotif.get_array(file=file, channels=('B8', 'B4')) for file in files] - else: - channels = [self.channel_indices['B8'], self.channel_indices['B4']] - images = [geotif.get_array(file=file, channels=channels) for file in files] - - ndvis = [utils.get_ndvi(nir=image[0], red=image[1]) for image in images] - - m = geotif.arrays_on_map( - files=files, arrays=ndvis, description='NDVI', cmap=self.cmap_ndvi, vmin=-1, vmax=1, m=m) - - # plot evi - if plot_evi: - if self.channel_indices is None: - images = [geotif.get_array(file=file, channels=('B8', 'B4', 'B2'))for file in files] - else: - channels = [self.channel_indices['B8'], self.channel_indices['B4'], self.channel_indices['B2']] - images = [geotif.get_array(file=file, channels=channels) for file in files] - - evis = [utils.get_evi(nir=image[0], red=image[1], blue=image[2]) for image in images] - - q = 0.02 # here, the quantile is applied to all EVIs (images) at once - evis_values = np.hstack([evi.reshape(-1) for evi in evis]) - - m = geotif.arrays_on_map( - files=files, - arrays=evis, - description='EVI', - cmap=self.cmap_ndvi, - vmin=np.quantile(evis_values, q=q), - vmax=np.quantile(evis_values, q=1-q), - m=m - ) - - return m - - def to_tif( - self, - *files, - output_folder, - - plot_rgb: bool=False, - plot_false_colors: bool=False, - plot_channels:list = None, - plot_ndvi: bool=False, - plot_evi: bool = False, - ): - - """ - Writes readable tif files (rgb or grey-channel) into given output_folder. - - :param files: list of file names or file dirs - :param output_folder: output folder in which tifs are stored - :param plot_rgb: bool if Sentinel-2 rgbs are plotted - :param plot_false_colors: bool if false color images of Sentinel-2 are plotted - :param plot_channels: list of channels that shall be plotted from Sentinel-2, e.g. ['B2', 'B3'] - :param plot_ndvi: bool if ndvi of Sentinel-2 images shall be plotted - :param plot_evi: bool if evi of Sentinel-2 images shall be plotted - :return: - """ - - files = self.get_file_dirs(*files) - self.output_folder = output_folder - - # plot rgb - if plot_rgb: - this_output_folder = tutils.files.join_paths(output_folder, 's2') - - if self.channel_indices is None: - geotif.image_to_tif( - *files, output_folder=this_output_folder, channels=['B4', 'B3', 'B2'], val_range=self.val_range_s2) - else: - channels = [self.channel_indices['B4'], self.channel_indices['B3'], self.channel_indices['B2']] - geotif.image_to_tif( - *files, output_folder=this_output_folder, channels=channels, val_range=self.val_range_s2) - - # plot false colors - if plot_false_colors: - this_output_folder = tutils.files.join_paths(output_folder, 's2_nir-r-g') - - if self.channel_indices is None: - geotif.image_to_tif( - *files, output_folder=this_output_folder, channels=['B8', 'B4', 'B3'], val_range=self.val_range_s2) - else: - channels = [self.channel_indices['B8'], self.channel_indices['B4'], self.channel_indices['B3']] - geotif.image_to_tif( - *files, output_folder=this_output_folder, channels=channels, val_range=self.val_range_s2) - - # plot channels - if plot_channels is not None: - for channel in plot_channels: - this_output_folder = tutils.files.join_paths(output_folder, f's2_{channel}') - - geotif.image_to_tif(*files, output_folder=this_output_folder, channels=channel, val_range='auto') - - # plot ndvi - if plot_ndvi: - this_output_folder = tutils.files.join_paths(output_folder, 's2_ndvi') - - for file in tqdm(files, desc='ndvi', leave=False): - if self.channel_indices is None: - image = geotif.get_array(file=file, channels=('B8', 'B4')) - else: - channels = [self.channel_indices['B8'], self.channel_indices['B4']] - image = geotif.get_array(file=file, channels=channels) - - ndvi = utils.get_ndvi(nir=image[0], red=image[1]) - - geotif.arrays_to_rgb_tif( - files=file, - arrays=ndvi, - output_folder=this_output_folder, - use_plot=True, - cmap=self.cmap_ndvi, - vmin=-1, - vmax=1, - disable_tqdm=True, - ) - - if plot_evi: - this_output_folder = tutils.files.join_paths(output_folder, 's2_evi') - - for file in tqdm(files, desc='evi', leave=False): - if self.channel_indices is None: - image = geotif.get_array(file=file, channels=('B8', 'B4', 'B2')) - else: - channels = [self.channel_indices['B8'], self.channel_indices['B4'], self.channel_indices['B2']] - image = geotif.get_array(file=file, channels=channels) - - evi = utils.get_evi(nir=image[0], red=image[1], blue=image[2]) - - q = 0.02 - geotif.arrays_to_rgb_tif( - files=file, - arrays=evi, - output_folder=this_output_folder, - use_plot=True, - cmap=self.cmap_ndvi, - vmin=np.quantile(evi, q=q), - vmax=np.quantile(evi, q=1-q), - disable_tqdm=True, - ) - -class S2SCLPlotter(Plotter): - - def __init__(self, folder, channel_indices: dict = None): - - self.folder = folder - self.channel_indices = channel_indices - - def plot(self, file, plot_scl_without_artifacts: bool = False, plot_legends: bool = False): - """ - Plots using matplotlib and plt.show(). No figure or axis is returned. - - :param file: file name or file dir - :param plot_scl_without_artifacts: bool if scl is plotted again, without artifacts - :return: - """ - - print(f'file: {file}') - file = self.get_file_dirs(file)[0] - - fig, ax = plt.subplots(1, 1, figsize=(self.fig_height, self.fig_height)) - - if file is not None: - ax = geotif.plot_landcover(file=file, channel='SCL', plot_legend=plot_legends, ax=ax) - ax.set_title('SCL') - else: - ax.set_title('SCL file not available') - ax.axis(False) - - plt.show() - - # plot scl (without artifacts) - if plot_scl_without_artifacts and file is not None: - if self.channel_indices is None: - geotif.plot_landcover(file=file, channel='SCL', plot_legend=False) - else: - channel = self.channel_indices['SCL'] - geotif.plot_landcover(file=file, channel=channel, plot_legend=False) - - plt.show() - - def plot_on_map( - self, - *files, - - plot_scl: bool = False, - plot_legends: bool = False, - plot_locations: bool = False, - - m=None, - **map_kwargs, - ): - """ - Plots on ipyleaflet map. - - :param files: list of file names or file dirs - :param plot_legend: bool if legends of land cover data shall be plotted with matplotlib - :param plot_locations: bool if polygons shall be plotted - :param m: ipyleaflet map; if None a new one is created - :param map_kwargs: ipyleaflet map kwargs - :return: - """ - - files = self.get_file_dirs(*files) - - if m is None: - m = utils.get_map(**map_kwargs) - - # plot locations - if plot_locations: - m = geotif.location_on_map(*files, m=m) - - # plot scl - if plot_scl: - if self.channel_indices is None: - m = geotif.landcover_on_map(*files, channel='SCL', plot_legend=plot_legends, m=m) - else: - channel = self.channel_indices['SCL'] - m = geotif.landcover_on_map(*files, channel=channel, plot_legend=plot_legends, m=m) - - return m - - def to_tif(self, *files, output_folder, plot_scl: bool = False): - - """ - Writes readable tif files (rgb or grey-channel) into given output_folder. - - :param files: list of file names or file dirs - :param output_folder: output folder in which tifs are stored - :return: - """ - - files = self.get_file_dirs(*files) - self.output_folder = output_folder - - if plot_scl: - - this_output_folder = tutils.files.join_paths(output_folder, 's2_scl') - - if self.channel_indices is None: - geotif.landcover_to_tif(*files, output_folder=this_output_folder, channel='SCL') - else: - channel = self.channel_indices['SCL'] - geotif.landcover_to_tif(*files, output_folder=this_output_folder, channel=channel) - - # legend - fig = landcover.plot_legend('SCL') - this_output_folder = tutils.files.join_paths(output_folder, 'legends') - fig.savefig(os.path.join(this_output_folder, f'legend_SCL.png'), facecolor='white', bbox_inches='tight') - plt.close() diff --git a/tlib/tgeo/s2.py b/tlib/tgeo/s2.py new file mode 100644 index 0000000000000000000000000000000000000000..63bd3d497dbe5962abf2e5172e3e39a0a27a649b --- /dev/null +++ b/tlib/tgeo/s2.py @@ -0,0 +1,520 @@ +import os + +import numpy as np +import matplotlib.pyplot as plt +from tqdm.auto import tqdm + +from . import geotif, landcover, utils +from . import kml as tgeo_kml +from tlib import tutils + + +class Plotter: + """ + This class helps to plot images, that were exported using the tgeo.gee.s2.AdvancedS2Exporter() class. + However, you can also plot other Sentinel-2 images, that were not exported with this class. + + :param folder_s2: folder which contains Sentinel-2 geotifs; if given, it is enough to use file names for plotting; + if None, you need the whole file paths + :param folder_scl: folder which contains geotifs with the SCL band of Sentinel-2; if None, no SCL can be plotted + :param folder_lcs: folder which contains geotifs with landcover dataset; if None, no landcover data can be plotted + :param folder_masks: folder which contains geotifs with masks; if None, no masks can be plotted + :param val_range_s2: value range that is plotted if rgb images are plotted of Sentinel-2 data + :param channel_indices: if your geotif channels have names (e.g. 'B2', 'B3') etc., you do not need this; + if your geotif channels do not have names, you need to define a dictionary here indicating the indices + of the channels, e.g. {'B2': 0, 'B3: 1, ...} + + Note: You can replace the following function variables to define more compplicated directories: + self.get_s2_dir + self.get_scl_dir + self.get_lcs_dir + self.get_masks_dir + """ + + def __init__( + self, + folder_s2=None, + folder_scl=None, + folder_lcs=None, + folder_masks=None, + val_range_s2: tuple = (0, 2**10), + channel_indices: dict = None + ): + + # functions to get directories + self.folder_s2 = folder_s2 + self.get_s2_dir = None if self.folder_s2 is None else lambda file: os.path.join(self.folder_s2, file) + self.get_scl_dir = None if folder_scl is None else lambda file: os.path.join(folder_scl, file) + self.get_lcs_dir = None if folder_lcs is None else lambda file: os.path.join(folder_lcs, file) + self.get_masks_dir = None if folder_masks is None else lambda file: os.path.join(folder_masks, file) + + # plot parameters + self.val_range_s2 = val_range_s2 + self.channel_indices = channel_indices + + self.cmap_mask = 'gray' + self.cmap_ndvi = 'RdYlGn' + self.output_folder = None # is determined in to_tif() function + + def get_file_dirs(self, *files): + """ + Return file paths according to given self folders for s2, scl, lcs and masks. + + :param files: list of files; if only file name (without directory) is given, paths for all self folders are + returned, if self folders are given; else, only path to files_s2 is returned and the others are None + :return: (list of files_s2, list of files_scl, list of files_lcs, list of files_masks) + """ + + # if only file name and not directory ... add folders to file name + + if files[0].split('/')[0] != self.folder_s2.split('/')[0] or files[0][0] == '/': + + files_s2 = [self.get_s2_dir(file) for file in files] if self.get_s2_dir is not None else None + files_scl = [self.get_scl_dir(file) for file in files] if self.get_scl_dir is not None else None + files_lcs = [self.get_lcs_dir(file) for file in files] if self.get_lcs_dir is not None else None + files_masks = [self.get_masks_dir(file) for file in files] if self.get_masks_dir is not None else None + + else: # assume that given file is s2 file and set others None + + files_s2 = files + files_scl = None + files_lcs = None + files_masks = None + + return files_s2, files_scl, files_lcs, files_masks + + def plot( + self, + file, + + plot_s2: bool = False, + plot_s2_channels: bool = False, + plot_scl_without_artifacts: bool = False, + plot_lcs: bool = False, + plot_legends: bool = False, + plot_masks: bool = False, + + fig_height: float = 4.6, + ): + """ + Plots using matplotlib and plt.show(). No figure or axis is returned. + + :param file: file name or file dir + :param plot_s2: bool if the following is plotted from Sentinel-2 image: rgb, false color, scl, ndvi, evi; + note, that the scl plot might have artifacts + :param s2_channels: bool if all channels of Sentinel-2 image are plotted individually + :param plot_scl_without_artifacts: bool if scl is plotted again, without artifacts + :param plot_lcs: bool if land cover data is plotted + :param plot_masks: bool if masks are plotted + :param fig_height: height of the plotted figures + :return: + """ + + # get file dirs + files_s2, files_scl, files_lcs, files_masks = self.get_file_dirs(file) + file_s2 = files_s2[0] if files_s2 is not None else None + file_scl = files_scl[0] if files_scl is not None else None + file_lcs = files_lcs[0] if files_lcs is not None else None + file_masks = files_masks[0] if files_masks is not None else None + + # print file name and path + print(f'file: {file}') + print(f'path: {file_s2}') + + # plot masks + if plot_masks and file_masks is not None: + channel_names = geotif.get_channel_names(file_masks) + if len(channel_names) == 1: # TODO: because of using plt.subplots, a single plot is not possible with this coding and the mask is doubled in this case + channel_names = [channel_names[0]] * 2 + + fig, axs = plt.subplots(1, len(channel_names), figsize=(fig_height * len(channel_names), fig_height)) + for i in range(len(channel_names)): + array = geotif.get_array(file_masks, channels=channel_names[i], transpose=True) + array = array.astype(bool) + + axs[i].imshow(array, cmap=self.cmap_mask, vmin=0, vmax=1) + axs[i].set_xticks([]), axs[i].set_yticks([]) + axs[i].set_title(channel_names[i]) + fig.tight_layout() + plt.show() + + # plot lcs + if plot_lcs and file_lcs is not None: + fig, axs = plt.subplots(1, 4, figsize=(fig_height * 4, fig_height)) + axs[0] = geotif.plot_landcover(file=file_lcs, channel='CORINE', plot_legend=plot_legends, ax=axs[0]) + axs[1] = geotif.plot_landcover(file=file_lcs, channel='MODIS_1', plot_legend=plot_legends, ax=axs[1]) + axs[2] = geotif.plot_landcover(file=file_lcs, channel='CGLS', plot_legend=plot_legends, ax=axs[2]) + axs[3] = geotif.plot_landcover(file=file_lcs, channel='GlobCover', plot_legend=plot_legends, ax=axs[3]) + + axs[0].set_title('CORINE') + axs[1].set_title('MODIS_1') + axs[2].set_title('CGLS') + axs[3].set_title('GlobCover') + fig.tight_layout() + plt.show() + + # plot s2 channels + if plot_s2_channels: + geotif.plot_channels(file_s2, ncols=5, figsize=(fig_height * 5, fig_height * 2)) + plt.show() + + # plot s2 + if plot_s2: + fig, axs = plt.subplots(1, 5, figsize=(fig_height * 5, fig_height)) + + # rgb image + try: # with channel names + axs[0] = geotif.plot_image( + file=file_s2, channels=['B4', 'B3', 'B2'], val_range=self.val_range_s2, ax=axs[0]) + except: # use channel indices + channels = [self.channel_indices['B4'], self.channel_indices['B3'], self.channel_indices['B2']] + axs[0] = geotif.plot_image(file=file_s2, channels=channels, val_range=self.val_range_s2, ax=axs[0]) + + axs[0].set_title('RGB') + + # false color image + try: # with channel names + axs[1] = geotif.plot_image( + file=file_s2, channels=['B8', 'B4', 'B3'], val_range=self.val_range_s2, ax=axs[1]) + except: # use channel indices + channels = [self.channel_indices['B8'], self.channel_indices['B4'], self.channel_indices['B3']] + axs[1] = geotif.plot_image(file=file_s2, channels=channels, val_range=self.val_range_s2, ax=axs[1]) + + axs[1].set_title('NIR-R-G') + + # ndvi + try: # with channel names + image = geotif.get_array(file=file_s2, channels=('B8', 'B4')) + except: # use channel indices + channels = [self.channel_indices['B8'], self.channel_indices['B4']] + image = geotif.get_array(file=file_s2, channels=channels) + + ndvi = utils.get_ndvi(nir=image[0], red=image[1]) + + axs[2] = tutils.plots.plot_original_image_size( + array=ndvi, plot_cbar=True, cmap=self.cmap_ndvi, vmin=-1, vmax=1, fig=fig, ax=axs[2]) + axs[2].set_title('NDVI') + + # evi + try: # with channel names + image = geotif.get_array(file=file_s2, channels=('B8', 'B4', 'B2')) + except: # use channel indices + channels = [self.channel_indices['B8'], self.channel_indices['B4'], self.channel_indices['B2']] + image = geotif.get_array(file=file_s2, channels=channels) + + evi = utils.get_evi(nir=image[0], red=image[1], blue=image[2]) + + q = 0.02 + + axs[3] = tutils.plots.plot_original_image_size( + array=evi, + plot_cbar=True, + cmap=self.cmap_ndvi, + vmin=np.quantile(evi, q=q), + vmax=np.quantile(evi, q=1-q), + fig=fig, + ax=axs[3]) + axs[3].set_title('EVI') + + # scl + if file_scl is not None: + axs[4] = geotif.plot_landcover(file=file_scl, channel='SCL', plot_legend=plot_legends, ax=axs[4]) + axs[4].set_title('SCL') + else: + axs[4].set_title('SCL file not available') + axs[4].axis(False) + + plt.show() + + # plot scl (without artifacts) + if plot_scl_without_artifacts and file_scl is not None: + try: # with channel name + geotif.plot_landcover(file=file_scl, channel='SCL', plot_legend=plot_legends) + except: # use channel index + channel = self.channel_indices['SCL'] + geotif.plot_landcover(file=file_scl, channel=channel, plot_legend=plot_legends) + + plt.show() + + def plot_on_map( + self, + *files, + + plot_rgb: bool=False, + plot_false_colors: bool=False, + plot_channels:list = None, + plot_ndvi: bool=False, + plot_evi: bool = False, + plot_scl: bool = False, + plot_lcs: list = None, # CORINE, MODIS_1, CGLS, GlobCover + plot_legends: bool = False, + plot_masks: list = None, # ROI, Valid_Area + plot_locations: bool = False, + + m=None, + **map_kwargs, + ): + """ + Plots on ipyleaflet map. + + :param files: list of file names or file dirs + :param plot_rgb: bool if Sentinel-2 rgbs are plotted + :param plot_false_colors: bool if false color images of Sentinel-2 are plotted + :param plot_channels: list of channels that shall be plotted from Sentinel-2, e.g. ['B2', 'B3'] + :param plot_ndvi: bool if ndvi of Sentinel-2 images shall be plotted + :param plot_evi: bool if evi of Sentinel-2 images shall be plotted + :param plot_scl: bool if scl band of Sentinel-2 shall be plotted + :param plot_lcs: list of land cover data that shall be plotted, e.g. ['CORINE', 'MODIS_1', 'CGLS', 'GlobCover'] + :param plot_legend: bool if legends of land cover data shall be plotted with matplotlib + :param plot_masks: list of mask channels that shall be plotted, e.g. ['ROI', 'Valid_Area'] + :param plot_locations: bool if polygons shall be plotted + :param m: ipyleaflet map; if None a new one is created + :param map_kwargs: ipyleaflet map kwargs + :return: + """ + + files_s2, files_scl, files_lcs, files_masks = self.get_file_dirs(*files) + + if m is None: + m = utils.get_map(**map_kwargs) + + # plot locations + if plot_locations: + m = geotif.location_on_map(*files_s2, m=m) + + # plot masks + if plot_masks is not None and files_masks is not None: + for mask in plot_masks: + m = geotif.image_on_map(*files_masks, channels=mask, val_range=(0, 1), m=m) + + # plot scl + if plot_scl and files_scl is not None: + try: # with channel names + m = geotif.landcover_on_map(*files_scl, channel='SCL', plot_legend=plot_legends, m=m) + except: # use channel indices + channel = self.channel_indices['SCL'] + m = geotif.landcover_on_map(*files_scl, channel=channel, plot_legend=plot_legends, m=m) + + # plot landcovers + if plot_lcs is not None and files_lcs is not None: + for lc in plot_lcs: + m = geotif.landcover_on_map(*files_lcs, channel=lc, plot_legend=plot_legends, m=m) + + # plot rgb + if plot_rgb: + try: # with channel names + m = geotif.image_on_map( + *files_s2, channels=['B4', 'B3', 'B2'], val_range=self.val_range_s2, m=m) # RGB + except: # use channel indices + channels = [self.channel_indices['B4'], self.channel_indices['B3'], self.channel_indices['B2']] + m = geotif.image_on_map(*files_s2, channels=channels, val_range=self.val_range_s2, m=m) # RGB + + # plot false_colors + if plot_false_colors: + try: # with channel names + m = geotif.image_on_map( + *files_s2, channels=['B8', 'B4', 'B3'], val_range=self.val_range_s2, m=m) # NIR-R-G + except: # use channel indices + channels = [self.channel_indices['B8'], self.channel_indices['B4'], self.channel_indices['B3']] + m = geotif.image_on_map(*files_s2, channels=channels, val_range=self.val_range_s2, m=m) # NIR-R-G + + # plot channels + if plot_channels is not None: + for channel in plot_channels: + m = geotif.image_on_map(*files_s2, channels=channel, val_range='auto', m=m) + + # plot ndvi + if plot_ndvi: + try: # with channel names + images = [geotif.get_array(file=file, channels=('B8', 'B4')) for file in files_s2] + except: # use channel indices + channels = [self.channel_indices['B8'], self.channel_indices['B4']] + images = [geotif.get_array(file=file, channels=channels) for file in files_s2] + + ndvis = [utils.get_ndvi(nir=image[0], red=image[1]) for image in images] + + m = geotif.arrays_on_map( + files=files_s2, arrays=ndvis, description='NDVI', cmap=self.cmap_ndvi, vmin=-1, vmax=1, m=m) + + # plot evi + if plot_evi: + try: # with channel names + images = [geotif.get_array(file=file, channels=('B8', 'B4', 'B2'))for file in files_s2] + except: # use channel indices + channels = [self.channel_indices['B8'], self.channel_indices['B4'], self.channel_indices['B2']] + images = [geotif.get_array(file=file, channels=channels) for file in files_s2] + + evis = [utils.get_evi(nir=image[0], red=image[1], blue=image[2]) for image in images] + + q = 0.02 # here, the quantile is applied to all EVIs (images) at once + evis_values = np.hstack([evi.reshape(-1) for evi in evis]) + + m = geotif.arrays_on_map( + files=files_s2, + arrays=evis, + description='EVI', + cmap=self.cmap_ndvi, + vmin=np.quantile(evis_values, q=q), + vmax=np.quantile(evis_values, q=1-q), + m=m + ) + + return m + + def to_tif( + self, + *files, + output_folder, + + plot_rgb: bool=False, + plot_false_colors: bool=False, + plot_channels:list = None, + plot_ndvi: bool=False, + plot_evi: bool = False, + plot_scl: bool = False, + plot_lcs: list = None, # CORINE, MODIS_1, CGLS, GlobCover + plot_masks: list=None, # ROI, Valid_Area + ): + + """ + Writes readable tif files (rgb or grey-channel) into given output_folder. + + :param files: list of file names or file dirs + :param output_folder: output folder in which tifs are stored + :param plot_rgb: bool if Sentinel-2 rgbs are plotted + :param plot_false_colors: bool if false color images of Sentinel-2 are plotted + :param plot_channels: list of channels that shall be plotted from Sentinel-2, e.g. ['B2', 'B3'] + :param plot_ndvi: bool if ndvi of Sentinel-2 images shall be plotted + :param plot_evi: bool if evi of Sentinel-2 images shall be plotted + :param plot_scl: bool if scl band of Sentinel-2 shall be plotted + :param plot_lcs: list of land cover data that shall be plotted, e.g. ['CORINE', 'MODIS_1', 'CGLS', 'GlobCover'] + :param plot_masks: list of mask channels that shall be plotted, e.g. ['ROI', 'Valid_Area'] + :return: + """ + + files_s2, files_scl, files_lcs, files_masks = self.get_file_dirs(*files) + + # plot rgb + if plot_rgb: + this_output_folder = tutils.files.join_paths(output_folder, 's2') + + try: # with channel names + geotif.image_to_tif( + *files_s2, output_folder=this_output_folder, channels=['B4', 'B3', 'B2'], val_range=self.val_range_s2) + except: # use channel indices + channels = [self.channel_indices['B4'], self.channel_indices['B3'], self.channel_indices['B2']] + geotif.image_to_tif( + *files_s2, output_folder=this_output_folder, channels=channels, val_range=self.val_range_s2) + + # plot false colors + if plot_false_colors: + this_output_folder = tutils.files.join_paths(output_folder, 's2_nir-r-g') + + try: # with channel names + geotif.image_to_tif( + *files_s2, output_folder=this_output_folder, channels=['B8', 'B4', 'B3'], val_range=self.val_range_s2) + except: # use channel indices + channels = [self.channel_indices['B8'], self.channel_indices['B4'], self.channel_indices['B3']] + geotif.image_to_tif( + *files_s2, output_folder=this_output_folder, channels=channels, val_range=self.val_range_s2) + + # plot channels + if plot_channels is not None: + for channel in plot_channels: + this_output_folder = tutils.files.join_paths(output_folder, f's2_{channel}') + + geotif.image_to_tif(*files_s2, output_folder=this_output_folder, channels=channel, val_range='auto') + + # plot ndvi + if plot_ndvi: + this_output_folder = tutils.files.join_paths(output_folder, 's2_ndvi') + + for file in tqdm(files_s2, desc='ndvi', leave=False): + try: # with channel names + image = geotif.get_array(file=file, channels=('B8', 'B4')) + except: # use channel indices + channels = [self.channel_indices['B8'], self.channel_indices['B4']] + image = geotif.get_array(file=file, channels=channels) + + ndvi = utils.get_ndvi(nir=image[0], red=image[1]) + + geotif.arrays_to_rgb_tif( + files=file, + arrays=ndvi, + output_folder=this_output_folder, + use_plot=True, + cmap=self.cmap_ndvi, + vmin=-1, + vmax=1, + disable_tqdm=True, + ) + + if plot_evi: + this_output_folder = tutils.files.join_paths(output_folder, 's2_evi') + + for file in tqdm(files_s2, desc='evi', leave=False): + try: # with channel names + image = geotif.get_array(file=file, channels=('B8', 'B4', 'B2')) + except: # use channel indices + channels = [self.channel_indices['B8'], self.channel_indices['B4'], self.channel_indices['B2']] + image = geotif.get_array(file=file, channels=channels) + + evi = utils.get_evi(nir=image[0], red=image[1], blue=image[2]) + + q = 0.02 + geotif.arrays_to_rgb_tif( + files=file, + arrays=evi, + output_folder=this_output_folder, + use_plot=True, + cmap=self.cmap_ndvi, + vmin=np.quantile(evi, q=q), + vmax=np.quantile(evi, q=1-q), + disable_tqdm=True, + ) + + # plot scl + if plot_scl and files_scl is not None: + this_output_folder = tutils.files.join_paths(output_folder, 's2_scl') + + try: # with channel names + geotif.landcover_to_tif(*files_scl, output_folder=this_output_folder, channel='SCL') + except: # use channel indices + channel = self.channel_indices['SCL'] + geotif.landcover_to_tif(*files_scl, output_folder=this_output_folder, channel=channel) + + # legend + fig = landcover.plot_legend('SCL') + this_output_folder = tutils.files.join_paths(output_folder, 'legends') + fig.savefig(os.path.join(this_output_folder, f'legend_SCL.png'), facecolor='white', bbox_inches='tight') + plt.close() + + # plot lcs + if plot_lcs is not None and files_lcs is not None: + for lc in plot_lcs: + this_output_folder = tutils.files.join_paths(output_folder, f'lcs_{lc}') + + geotif.landcover_to_tif(*files_lcs, output_folder=this_output_folder, channel=lc) + + # legend + fig = landcover.plot_legend(lc) + this_output_folder = tutils.files.join_paths(output_folder, 'legends') + fig.savefig(os.path.join(this_output_folder, f'legend_{lc}.png'), facecolor='white', bbox_inches='tight') + plt.close() + + # plot masks + if plot_masks is not None and files_masks is not None: + for mask in plot_masks: + this_output_folder = tutils.files.join_paths(output_folder, f'masks_{mask}') + + geotif.image_to_tif(*files_masks, output_folder=this_output_folder, channels=mask, val_range=(0, 1)) + + def tif_to_kml(self, main_folder=None, save=True, kml=None): + + if main_folder == None: + main_folder = self.output_folder + + kml = tgeo_kml.subfolders_to_kml( + main_folder=main_folder, overlay_folder='legends', points_for='s2', polygons_for='s2', save=save, kml=kml) + + return kml