diff --git a/idepix-aatsr/.gitignore b/idepix-aatsr/.gitignore new file mode 100644 index 00000000..eba1e305 --- /dev/null +++ b/idepix-aatsr/.gitignore @@ -0,0 +1,15 @@ +target/ +.idea +*.iml +Technical Note_cloudshadow_draft_v0.2.docx +pom.xml.tag +pom.xml.releaseBackup +pom.xml.versionsBackup +pom.xml.next +release.properties +dependency-reduced-pom.xml +buildNumber.properties +.mvn/timing.properties + +# Avoid ignoring Maven wrapper jar file (.jar files are usually ignored) +!/.mvn/wrapper/maven-wrapper.jar diff --git a/idepix-aatsr/pom.xml b/idepix-aatsr/pom.xml new file mode 100644 index 00000000..c06fcfcb --- /dev/null +++ b/idepix-aatsr/pom.xml @@ -0,0 +1,95 @@ + + + + + 4.0.0 + + org.esa.snap + snap-idepix + 8.0.1 + + + idepix-aatsr + 8.0.3 + + nbm + + IdePix AATSR + Classification of pixels (cloud, cloud-shadow and land) originating from AATSR 4th reprocessing data. + + + 8.0.1 + 8.0.6 + + + + + org.esa.snap + idepix-core + ${idepix.version} + + + + org.esa.snap + ceres-core + + + org.esa.snap + ceres-jai + + + org.esa.snap + ceres-glayer + + + org.esa.snap + snap-core + + + org.esa.snap + snap-collocation + + + org.esa.snap + snap-gpf + + + + org.esa.s3tbx + s3tbx-sentinel3-reader + ${s3tbx.version} + + + + junit + junit + + + + + + + org.apache.netbeans.utilities + nbm-maven-plugin + + + + + + + diff --git a/idepix-aatsr/resource/aatsr_cloud_shadow_dev.py b/idepix-aatsr/resource/aatsr_cloud_shadow_dev.py new file mode 100644 index 00000000..6550edd9 --- /dev/null +++ b/idepix-aatsr/resource/aatsr_cloud_shadow_dev.py @@ -0,0 +1,763 @@ +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import sys +import os +import time +from scipy import interpolate +from scipy.signal import fftconvolve + +sys.path.append("C:\\Users\Dagmar\.snap\snap-python") +import snappy as snp +from snappy import Product +from snappy import ProductData +from snappy import ProductIO +from snappy import jpy +from snappy import PixelPos #org.esa.snap.core.datamodel.PixelPos + +GeoPos = jpy.get_type('org.esa.snap.core.datamodel.GeoPos') +PixelPos = jpy.get_type('org.esa.snap.core.datamodel.PixelPos') +GPF = jpy.get_type('org.esa.snap.core.gpf.GPF') +File = jpy.get_type('java.io.File') +ProgressMonitor = jpy.get_type('com.bc.ceres.core.ProgressMonitor') + +def getRefinedHeightFromCtp(ctp, slp, temperatures): #ctp: cloud top pressure, slp: sea level pressure, temperature-profile + height = 0.0 + prsLevels = np.array((1000., 950., 925., 900., 850., 800., 700., 600., 500., 400., 300., 250., 200., 150., 100., + 70., 50., 30., 20., 10., 7., 5., 3., 2., 1.)) + + if ctp >= prsLevels[-1]: + for i in range(len(prsLevels)- 1): + if ctp > prsLevels[0] or (ctp < prsLevels[i] and ctp > prsLevels[i + 1]): + t1 = temperatures[i] + t2 = temperatures[i + 1] + ts = (t2 - t1) / (prsLevels[i + 1] - prsLevels[i]) * (ctp - prsLevels[i]) + t1 + height = getHeightFromCtp(ctp, slp, ts) + break + + else: + # CTP < 1 hPa? This should never happen... + t1 = temperatures[-2] + t2 = temperatures[-1] + ts = (t2 - t1) / (prsLevels[-2] - prsLevels[-1]) * (ctp - prsLevels[-1]) + t1 + height = getHeightFromCtp(ctp, slp, ts) + + return height + +def getRefinedHeightFromBT(BT, slp, temperatures): #BT: brightness temperature, slp: sea level pressure, temperature-profile + #this function is not necessarily eindeutig! inversion in T(p) can lead to two (or more) solutions for p. + height1, height2, height1b, height2b = None, None, None, None + prsLevels = np.array((1000., 950., 925., 900., 850., 800., 700., 600., 500., 400., 300., 250., 200., 150., 100., + 70., 50., 30., 20., 10., 7., 5., 3., 2., 1.)) + + if BT >= np.min(temperatures) and BT <= np.max(temperatures): + for i in range(len(temperatures)- 1): + if BT > temperatures[i] and BT < temperatures[i + 1]: + ctp1 = prsLevels[i] + ctp2 = prsLevels[i + 1] + ctp = (ctp2 - ctp1) / (temperatures[i + 1] - temperatures[i]) * (BT - temperatures[i]) + ctp1 + height1 = getHeightFromCtp(ctp, slp, BT) + height1b = computeHeightFromPressure(ctp) + if BT < temperatures[i] and BT > temperatures[i + 1]: #inversion! + ctp1 = prsLevels[i] + ctp2 = prsLevels[i + 1] + ctp = (ctp2 - ctp1) / (temperatures[i + 1] - temperatures[i]) * (BT - temperatures[i]) + ctp1 + height2 = getHeightFromCtp(ctp, slp, BT) + height2b = computeHeightFromPressure(ctp) + + return height1, height2, height1b, height2b + +def getHeightFromCtp(ctp, p0, ts): + return -ts * (np.power(ctp / p0, 1. / 5.255) - 1) / 0.0065 + +def computeHeightFromPressure(pressure): + return -8000 * np.log(pressure / 1013.0) + +def investigate_AATSR4th_transect(plot_BT=False, plot_Temp=False): + path ="E:\Documents\projects\QA4EO\data\\transect\\" + # fname = "subset_0_of_ENV_AT_1_RBT____20021129T235200_20021130T013735_20210315T024827_6334_011_359______DSI_R_NT_004_TRANSECT.txt" + # todo: unfortunately, no radiance or BT is extracted! + fname = "subset_WATER_of_ENV_AT_1_RBT____20020810T083508_20020810T102042_20210303T040313_6334_008_264______DSI_R_NT_004_TRANSECT.txt" + # fname = "subset_LAND_of_ENV_AT_1_RBT____20020810T083508_20020810T102042_20210303T040313_6334_008_264______DSI_R_NT_004_TRANSECT.txt" + ##data from the snap-transect plot! Combine both for analysis. + fname_BT = "subset_WATER_S8_S9_of_ENV_AT_1_RBT____20020810T083508_TRANSECT.txt" + + d = pd.read_csv(path + fname, header=0, comment='#', sep='\t') + d_BT = pd.read_csv(path + fname_BT, header=0, comment='#', sep='\t') + print(d.shape) + print(d_BT.shape) + print(d.columns.values) + referencePressureLevels = np.array((1000., 950., 925., 900., 850., 800., 700., 600., 500., 400., 300., 250., 200., + 150., 100., 70., 50., 30., 20., 10., 7., 5., 3., 2., 1.)) + + col_BT = [a for a in d_BT.columns.values if 'BT' in a] + col_BT = [a for a in col_BT if not 'sigma' in a] + col_T = [a for a in d.columns.values if 'temperature_profile_tx_pressure' in a] + col_slp = [a for a in d.columns.values if 'surface_pressure' in a] + print(col_slp) + if plot_BT: + for v in col_BT: + plt.plot(d_BT[v].values, '-', label=v) + plt.legend() + plt.show() + + if plot_Temp: + ID = np.array(d_BT[col_BT[0]].values < 270) #todo: cloud flags? + print(np.sum(ID)) + d_BT = d_BT.loc[ID,:] + d_BT = d_BT.reset_index(drop=True) + d = d.loc[ID,:] + d = d.reset_index(drop=True) + number = np.random.random(10)*d.shape[0] + for i in number: + fig, ax = plt.subplots(1, 1, figsize=(5,5)) + ax.plot(d[col_T].loc[int(i),:].values, referencePressureLevels, '-') + for v in col_BT: + BT = d_BT[v].values[int(i)] + ax.plot((BT, BT), (np.min(referencePressureLevels), np.max(referencePressureLevels)), 'r-') + height1, height2, h1, h2 = getRefinedHeightFromBT(BT, d[col_slp].values[int(i)], d[col_T].loc[int(i),:].values) + print(i, BT, height1, height2, h1, h2) + ax.invert_yaxis() + ax.set_xlabel('Temperature T [K]') + ax.set_ylabel('pressure [hPa]') + fig.tight_layout() + plt.show() + +def get_band_or_tiePointGrid(product, name, dtype='float32', reshape=True, subset=None): + ## + # This function reads a band or tie-points, identified by its name , from SNAP product + # The fuction returns a numpy array of shape (height, width) + ## + if subset is None: + height = product.getSceneRasterHeight() + width = product.getSceneRasterWidth() + sline, eline, scol, ecol = 0, height-1, 0, width -1 + else: + sline,eline,scol,ecol = subset + height = eline - sline + 1 + width = ecol - scol + 1 + + var = np.zeros(width * height, dtype=dtype) + if name in list(product.getBandNames()): + product.getBand(name).readPixels(scol, sline, width, height, var) + elif name in list(product.getTiePointGridNames()): + var.shape = (height, width) + for i,iglob in enumerate(range(sline,eline+1)): + for j,jglob in enumerate(range(scol,ecol+1)): + var[i, j] = product.getTiePointGrid(name).getPixelDouble(jglob, iglob) + var.shape = (height*width) + else: + raise Exception('{}: neither a band nor a tie point grid'.format(name)) + + if reshape: + var.shape = (height, width) + + return var + +def set_cloud_mask_from_Bayesian(cloudB): + single_moderate_mask = np.bitwise_and(cloudB, 2 ** 1) == 2 ** 1 + no_bayes_available_mask = np.bitwise_and(cloudB, 2 ** 7) == 2 ** 7 + out = np.zeros(cloudB.shape) + np.nan + if np.sum(no_bayes_available_mask)==0: + out = single_moderate_mask + return out + + +def set_cloud_mask_from_Confidence(cloudB): + summary_cloud_mask = np.bitwise_and(cloudB, 2 ** 14) == 2 ** 14 + return summary_cloud_mask + + +def set_cloud_mask_from_cloudBand(cloudB): + visible_mask = np.bitwise_and(cloudB, 2 ** 0) == 2 ** 0 + gross_cloud_mask = np.bitwise_and(cloudB, 2 ** 7) == 2 ** 7 + thin_cirrus_mask = np.bitwise_and(cloudB, 2 ** 8) == 2 ** 8 + medium_high_mask = np.bitwise_and(cloudB, 2 ** 9) == 2 ** 9 + return np.array(visible_mask + gross_cloud_mask + thin_cirrus_mask + medium_high_mask > 0) + +def set_land_mask(product, subset=None): + confid_in = get_band_or_tiePointGrid(product, 'confidence_in', 'int32', subset=subset) + #LAND : coastline, tidal, land, inland_water + coastline_mask = np.bitwise_and(confid_in, 2 ** 0) == 2 ** 0 + tidal_mask = np.bitwise_and(confid_in, 2 ** 2) == 2 ** 2 + land_mask = np.bitwise_and(confid_in, 2 ** 3) == 2 ** 3 + inland_water_mask = np.bitwise_and(confid_in, 2 ** 4) == 2 ** 4 + return np.array(coastline_mask + tidal_mask + land_mask + inland_water_mask > 0) + +def calculate_view_azimuth_interpolation_singleline( vaa_line, width, plot=False): + + nx0 = 0 + nx1 = 200 + nx2 = 340 + nx3 = width + nx_change = 272 # nadir line between 271 + 272 + + vaa_out = np.copy(vaa_line) + + ## + # interpolation for row ny: + y1 = vaa_line[nx0:nx1] + y2 = vaa_line[nx2:] + x1 = np.arange(nx0, nx1, 1.) + x2 = np.arange(nx2, nx3, 1.) + + p = np.polyfit(x1, y1, 2) + print(p) + delta = vaa_line[nx1] - (p[0] * (nx1) ** 2 + p[1] * nx1 + p[2]) + print(vaa_line[nx1]) + print(p[0] * (nx1) ** 2 + p[1] * nx1 + p[2]) + x_out = np.arange(nx1, nx_change, 1.) + y_out1 = p[0] * (x_out) ** 2 + p[1] * x_out + p[2] + delta + vaa_out[nx1:nx_change] = y_out1 + + p = np.polyfit(x2, y2, 2) + delta = vaa_line[nx2] - (p[0] * (nx2) ** 2 + p[1] * nx2 + p[2]) + x_out = np.arange(nx_change, nx2, 1.) + y_out2 = p[0] * (x_out) ** 2 + p[1] * x_out + p[2] + delta + + vaa_out[nx_change:nx2] = y_out2 + + print(vaa_out[(nx_change-1):(nx_change+1)]) + if plot: + plt.plot(vaa_line, '-', label='from TPG') + plt.plot(vaa_out, 'r-', label='Interpolation') + plt.xlabel('Pixel X') + plt.ylabel('View Azimuth Angle') + plt.legend() + plt.show() + + return vaa_out + + +def calculate_view_zenith_interpolation_singleline_OLCI(vza, ny): + + nx0 = 3500 + nx1 = 3580 + nx2 = 3649 + nx3 = nx2 + 80 + nx_change = 3616 + + vza_out = np.copy(vza) + + ## + # interpolation for row ny: + + y1 = vza[ny, nx0:nx1] + y2 = vza[ny, nx2:nx3] + x1 = np.arange(nx0, nx1, 1.) + x2 = np.arange(nx2, nx3, 1.) + + + p = np.polyfit(x1, y1, 2) + x_out = np.arange(nx1, nx_change, 1.) + y_out1 = p[0] * (x_out) ** 2 + p[1] * x_out + p[2] + vza_out[ny, nx1:nx_change] = y_out1 + + p = np.polyfit(x2, y2, 2) + x_out = np.arange(nx_change, nx2, 1.) + y_out2 = p[0] * (x_out) ** 2 + p[1] * x_out + p[2] + + vza_out[ny, nx_change:nx2] = y_out2 + + return vza_out + +def setRelativePathIndex_and_TheoreticalHeight(sza, saa, oza, x_tx, orientation, spatialResolution, maxObjectAltitude, minSurfaceAltitude): #oaa replaced by x_tx + shadow_angle_PixelCoord = (saa - orientation) + 180. + cosSaa = np.cos(shadow_angle_PixelCoord*np.pi/180. - np.pi / 2.) + sinSaa = np.sin(shadow_angle_PixelCoord*np.pi/180. - np.pi / 2.) + + #modified sza for influence of oza: + if saa - orientation < 180.: + if x_tx < 0: + sza = np.arctan(np.tan(sza * np.pi / 180.) - np.tan(oza * np.pi / 180.)) * 180. / np.pi #negative?? + else: + sza = np.arctan(np.tan(sza * np.pi / 180.) + np.tan(oza * np.pi / 180.)) * 180. / np.pi + else: + if x_tx < 0: + sza = np.arctan(np.tan(sza * np.pi / 180.) + np.tan(oza * np.pi / 180.)) * 180. / np.pi #negative + else: + sza = np.arctan(np.tan(sza * np.pi / 180.) - np.tan(oza * np.pi / 180.)) * 180. / np.pi #positive + + + deltaProjX = ((maxObjectAltitude - minSurfaceAltitude) * np.tan(sza*np.pi/180.) * cosSaa) / spatialResolution + deltaProjY = ((maxObjectAltitude - minSurfaceAltitude) * np.tan(sza*np.pi/180.) * sinSaa) / spatialResolution + + x0 = 0 + y0 = 0 + + x1 = x0 + deltaProjX + np.sign(deltaProjX) * 1.5 + y1 = y0 + deltaProjY + np.sign(deltaProjY) * 1.5 + + #create index steps + # Path touches which pixels? + #setup all pixel centers from x0/y0 to x1/y1. + # calculate distance between pixel center and line (X0, X1) + # all distances below/equal sqrt(2*0.5^2): the pixel is touched by the line and a potential shadow pixel. + + if x0 0) > 0: + eline = np.arange(0, height-1, 1)[diffSZATest>0] + eline = eline[eline > sline] + eline = np.min(eline) + else: + eline = height + sline = int(sline) + eline = int(eline) + # eline = np.max(np.arange(0, height, 1)[SZATest]) + sza = sza[sline:eline,:] + subset = (int(sline), int(eline-1), 0, width-1) + saa = get_band_or_tiePointGrid(product, 'solar_azimuth_tn', subset=subset) + oza = get_band_or_tiePointGrid(product, 'sat_zenith_tn', subset=subset) + x_tx = get_band_or_tiePointGrid(product, 'x_tx', subset=subset) # at nadir line x_tx changes sign. x<272 : x_tx<0; x>=272 : x_tx >0 + + ## OAA is not needed! use fixed position of nadir line between x = 271, 272 + # oaa = get_band_or_tiePointGrid(product, 'sat_azimuth_tn') #needs adjustment! + # test sat azimuth interpolation. + # a = calculate_view_azimuth_interpolation_singleline(oaa[0,:], width, plot=True) + # oaa_correct = np.zeros((height, width)) + # for i in range(height): + # oaa_correct[i,:] = a + + + ## derive orientation at coarse raster positions. + ## interpolate/extrapolate spatially to entire grid + GridStep = 100 + X, Y = np.mgrid[1:(width-1):GridStep, 1:(eline-sline-1):GridStep] + print("start Position") + orientationS = np.zeros(X.shape) + for iy in range(orientationS.shape[0]): + for ix in range(orientationS.shape[1]): + x = X[iy, ix] + y = Y[iy, ix] + pixelPos1 = PixelPos(x-1 + 0.5, sline + y + 0.5) + geoPos = product.getSceneGeoCoding().getGeoPos(pixelPos1, None) + lat1 = geoPos.getLat() + lon1 = geoPos.getLon() + pixelPos2 = PixelPos(x + 1 + 0.5, sline + y + 0.5) + geoPos = product.getSceneGeoCoding().getGeoPos(pixelPos2, None) + lat2 = geoPos.getLat() + lon2 = geoPos.getLon() + orientationS[iy, ix] = calculate_orientation(lat1, lon1, lat2, lon2) * 180. / np.pi + + f = interpolate.RectBivariateSpline(y=np.arange(1,(eline-sline-1), GridStep), x=np.arange(1, (width-1), GridStep), + z=orientationS, bbox=[0, width, 0, (eline-sline)] , kx=3, ky=3, s=0) + xx = np.arange(0, width, 1) + yy = np.arange(0, (eline-sline), 1) + orientation = np.transpose(f(xx, yy)) # orientationFull, + + ### by Pixel: Lat, Lon, Orientation + # deriving North direction on pixel grid, orientation + # lat = np.zeros((eline - sline , width)) + # lon = np.zeros((eline - sline , width)) + # print("start Position") + # for ix, x in enumerate(range(0, width)): + # for iy, y in enumerate(range(sline, eline)): + # pixelPos = PixelPos(x + 0.5, y + 0.5) + # geoPos = product.getSceneGeoCoding().getGeoPos(pixelPos, None) + # lat[iy, ix] = geoPos.getLat() + # lon[iy, ix] = geoPos.getLon() + # + # #direction of North on the grid. Used to adjust the sun azimuth angle to grid directions. + # #in degree! + # print("start orientation") + # orientation = np.zeros(lat.shape) + # for i in range(orientation.shape[0]): + # for j in range(width - 2): + # orientation[i, j + 1] = calculate_orientation(lat[i, j], lon[i, j], lat[i, j + 2], lon[i, j + 2]) * 180./np.pi + # orientation[:, 0] = orientation[:, 1] + # orientation[:, width-1] = orientation[:, width-2] + + ### plot orientation comparison + # fig, ax = plt.subplots(1, 2, figsize=(8,5)) + # im = ax[0].imshow(orientationFull) + # fig.colorbar(im, ax=ax[0], shrink=0.5) + # ax[0].set_title('orient. Interpolated') + # im=ax[1].imshow(orientation) + # fig.colorbar(im, ax=ax[1], shrink=0.5) + # ax[1].set_title('all pixl orientation') + # plt.show() + # + # fig,ax = plt.subplots(1, 2, figsize=(8,5)) + # im = ax[0].imshow((orientationFull- orientation)/orientation*100.) + # fig.colorbar(im, ax=ax[0], shrink=0.5) + # ax[0].set_title('rel.Diff orient. Interpolated %') + # im = ax[1].plot(orientation[:, 1], 'r-', label='all') + # im = ax[1].plot(orientationFull[:, 1], 'b--', label='interpol.') + # im = ax[1].plot( Y[0,:], orientationS[0, :], 'k+', label='points') + # ax[1].legend() + # plt.show() + fig,ax = plt.subplots(1, 2, figsize=(8,5)) + im = ax[0].imshow(orientation) + fig.colorbar(im, ax=ax[0], shrink=0.5) + ax[0].set_title('orient. Interpolated ') + im = ax[1].plot(orientation[:, 1], 'r-', label='all') + im = ax[1].plot(Y[0,:], orientationS[0, :], 'k+', label='points') + ax[1].legend() + plt.show() + + ### cloud mask + + if cloudflagType == 'new': + bayes_in = get_band_or_tiePointGrid(product, 'bayes_in', 'int32', subset=subset) + bayesMask = set_cloud_mask_from_Bayesian(bayes_in) + + confid_in = get_band_or_tiePointGrid(product, 'confidence_in', 'int32', subset=subset) + confidMask = set_cloud_mask_from_Confidence(confid_in) + + if np.sum(np.isnan(bayesMask)) > 0: + cloudMask = confidMask + else: + #todo: what to do with partial information in bayesian cloud mask? + print('cloud mask to be done') + cloudMask = np.logical_or(bayesMask, confidMask) + + else: + cloud_in = get_band_or_tiePointGrid(product, 'cloud_in', 'int32', subset=subset) + cloudMask = set_cloud_mask_from_cloudBand(cloud_in) + + ## landmask + landMask = set_land_mask(product, subset=subset) + + ### convolution cloudmask and search radius, convolution landmask and search radius + # every 1000 or 2000 pixels (y-direction), the convolution is done with the current, mean search radius defined by SZA and CTH. + + ConvolStep = 2000 + upperLimits = np.arange(0, cloudMask.shape[0], ConvolStep) + if len(upperLimits)> 1: + upperLimits[-1] = cloudMask.shape[0] + else: + upperLimits = np.array((0, cloudMask.shape[0])) + + startSearchMask = np.zeros((cloudMask.shape)) + landConvolveMask = np.zeros((landMask.shape)) + for i, up in enumerate(upperLimits[1:]): + #setup kernel + radius = cth * np.tan(np.median(sza[upperLimits[i]:up,:])*np.pi/180.) + print(i, radius) + kernel = setup_round_kernel(radius=radius, spacing=(1000.,1000.)) #todo: use an elongated shape in direction of illumination! + # kNy, kNx = kernel.shape + # convolveMatrixSubset = cloudMask[upperLimits[i]:up,:] + # convolveMatrix = np.zeros((convolveMatrixSubset.shape[0]+2*kernel.shape[0], convolveMatrixSubset.shape[1]+2*kernel.shape[1])) + # convolveMatrix[kNy:(-kNy), kNx:(-kNx)] = convolveMatrixSubset + # for iy in range(kNy): + # convolveMatrix[iy,kNx:(-kNx)] = convolveMatrixSubset[0,:] + # convolveMatrix[-(iy+1), kNx:(-kNx)] = convolveMatrixSubset[-1,:] + # for ix in range(kNx): + # convolveMatrix[:, ix] = convolveMatrix[:, kNx] + # convolveMatrix[:, -(ix + 1)] = convolveMatrix[:, -(kNx+1)] + # + # perc_circle = fftconvolve(convolveMatrix, kernel, mode='same') + startSearchMask[upperLimits[i]:up, :] = convolve_mask_kernel(cloudMask[upperLimits[i]:up, :], kernel) + landConvolveMask[upperLimits[i]:up, :] = convolve_mask_kernel(landMask[upperLimits[i]:up, :], kernel) + + startSearchMask = np.logical_and(startSearchMask >0.001, startSearchMask < 0.998) + startSearchMask = np.logical_and(startSearchMask, landConvolveMask > 0.001) + print(np.sum(startSearchMask), np.sum(cloudMask)) + print(np.sum(np.logical_and(startSearchMask, cloudMask==1))) + + fig, ax = plt.subplots(1, 2, figsize=(5, 8)) + ax[0].imshow(cloudMask) + ax[0].set_title('AATSR cloud mask') + ax[1].imshow(startSearchMask) + ax[1].set_title('Start Points') + plt.show() + + elevation = get_band_or_tiePointGrid(product, 'elevation_in', subset=subset) + spatialResolution = 1000. # todo: read from product + # maxObjectAltitude = np.nanmax(elevation) + minSurfaceAltitude = np.nanmin(elevation) + ShadowArray = np.zeros(elevation.shape) + + ### + # print(filename) + # print('SZA', sza[yline, :]) + # print('orientation', orientation[yline, :]) + # print('minSurfaceAltitude', minSurfaceAltitude) + # print('startSearchMask', startSearchMask[yline, :]) + + this_height = sza.shape[0] + this_width = sza.shape[1] + count = 0 + startTime = time.time() + for i in range(this_height): + for j in range(this_width): + if cloudMask[i, j]==1 and startSearchMask[i, j]: + count += 1 + if count % 10000 == 0: + print(count) + #search for cloud shadow. + # calculate theoretical height at search path position. + illuPathSteps, illuPathHeight, threshHeight = setRelativePathIndex_and_TheoreticalHeight(sza=sza[i, j], saa=saa[i, j], + oza=oza[i, j], x_tx = x_tx[i, j],#oaa = oaa_correct[i, j], + orientation = orientation[i, j], + spatialResolution=spatialResolution, + maxObjectAltitude=cth, + minSurfaceAltitude=minSurfaceAltitude) + X = np.sqrt(illuPathSteps[:, 0] ** 2 + illuPathSteps[:, 1] ** 2) + + IndexArray = np.copy(illuPathSteps) + IndexArray[:, 0] += j + IndexArray[:, 1] += i + # find cloud free positions along the search path: + ID = np.logical_and(np.logical_and(IndexArray[:, 0] >= 0, IndexArray[:, 0] < this_width), + np.logical_and(IndexArray[:, 1] >= 0, IndexArray[:, 1] < this_height)) + + if np.sum(ID) > 3: # path positions + IndexArray = IndexArray[ID, :] + BaseHeightArray = illuPathHeight[ID] + # Xx = X[ID] + elevPath = elevation[IndexArray[:, 1], IndexArray[:, 0]] + + # plt.plot(Xx, elevPath, 'g+') + # plt.plot(Xx, BaseHeightArray, 'rx') + # plt.show() + ID2 = np.logical_and(np.abs(BaseHeightArray - elevPath) < threshHeight , + np.logical_not(cloudMask[IndexArray[:, 1], IndexArray[:, 0]])) + + ShadowArray[IndexArray[ID2, 1], IndexArray[ID2, 0]] = 1 + + # plt.imshow(ShadowArray) + # plt.show() + # + # out = cloudMask + ShadowArray *2 + # plt.imshow(out) + # plt.show() + endTime = time.time() + + print("Time", startTime, endTime) + + outpath = scene_path + filename[:42] + "_testCloudShadow_Cloud_" +cloudflagType+ ".dim" + outProduct = Product('AASTR_cloudShadow', 'AASTR_cloudShadow', width, this_height) #height + outProduct.setFileLocation(File(outpath)) + + ProductSubsetDef = jpy.get_type('org.esa.snap.core.dataio.ProductSubsetDef') + subset_def = ProductSubsetDef() + subset_def.setRegion(0, sline, width, eline) # (0,0,width, height) + product.transferGeoCodingTo(outProduct, subset_def) + + rad = get_band_or_tiePointGrid(product, 'S2_radiance_in', subset=subset) + band = outProduct.addBand("S2_radiance_in", ProductData.TYPE_FLOAT32) + band.setNoDataValue(np.nan) + band.setNoDataValueUsed(True) + sourceData = rad.reshape(cloudMask.shape).astype('float32') + band.setRasterData(ProductData.createInstance(sourceData)) + + band = outProduct.addBand("cloudMask", ProductData.TYPE_INT16) + band.setNoDataValue(np.nan) + band.setNoDataValueUsed(True) + sourceData = cloudMask.reshape(cloudMask.shape).astype('int16') + band.setRasterData(ProductData.createInstance(sourceData)) + + band = outProduct.addBand("shadowMask", ProductData.TYPE_INT16) + band.setNoDataValue(np.nan) + band.setNoDataValueUsed(True) + sourceData = ShadowArray.reshape(cloudMask.shape).astype('int16') + band.setRasterData(ProductData.createInstance(sourceData)) + + ProductIO.writeProduct(outProduct, outpath, 'BEAM-DIMAP') + + product.closeIO() + outProduct.closeIO() + + +def analyse_AATSR4th_transect(varname='sat_azimuth_tn', start = 10950, end = 12950, step=50): + scene_path = "E:\Documents\projects\QA4EO\AATSR4th\\" + filename = "ENV_AT_1_RBT____20020810T083508_20020810T102042_20210303T040313_6334_008_264______DSI_R_NT_004.dim" + # filename = "subset_0_of_ENV_AT_1_RBT____20020810T083508_20020810T102042_20210303T040313_6334_008_264______DSI_R_NT_004.dim" + # filename = "subset_SouthernHem_of_ENV_AT_1_RBT____20021129T235200_20021130T013735_20210315T024827_6334_011_359______DSI_R_NT_004.dim" + # filename = "subset_LowSunNorth_of_ENV_AT_1_RBT____20020810T083508_20020810T102042_20210303T040313_6334_008_264______DSI_R_NT_004.dim" + product = snp.ProductIO.readProduct(os.path.join(scene_path, filename)) + height = product.getSceneRasterHeight() + width = product.getSceneRasterWidth() + + # sza = get_band_or_tiePointGrid(product, 'solar_zenith_tn') # sline,eline,scol,ecol = subset + subset = (start, end, 0, width - 1) + vaa = get_band_or_tiePointGrid(product, varname, subset=subset) + + + yList = np.arange(start, end, step, dtype='int16') + latPart = np.zeros((len(yList), width)) + lonPart = np.zeros((len(yList), width)) + orientationPart = np.zeros((len(yList), width)) + + for j, y in enumerate(yList): + for ix, x in enumerate(range(0, width)): + pixelPos = PixelPos(x + 0.5, y + 0.5) + geoPos = product.getSceneGeoCoding().getGeoPos(pixelPos, None) + latPart[j, ix] = geoPos.getLat() + lonPart[j, ix] = geoPos.getLon() + + for i in range(orientationPart.shape[0]): + for j in range(width - 2): + orientationPart[i, j + 1] = calculate_orientation(latPart[i, j], lonPart[i, j], latPart[i, j + 2], + lonPart[i, j + 2]) * 180. / np.pi + orientationPart[:, 0] = orientationPart[:, 1] + orientationPart[:, width - 1] = orientationPart[:, width - 2] + + fig, ax = plt.subplots(1, 2, figsize=(8,5)) + for j,y in enumerate(yList): + ax[0].plot(vaa[y-start, :], '-', label=y) + ax[1].plot(vaa[y-start, :]-orientationPart[j,:], '-', label=y) + ax[0].set_xlabel('pixel No X') + ax[0].set_ylabel('view azimuth angle') + ax[1].set_xlabel('pixel No X') + ax[1].set_ylabel('view azimuth angle - orientation') + plt.show() + + # yList = np.arange(end, end+2000, 50, dtype='int16') + # for y in yList: + # plt.plot(vaa[y,:], '-', label=y) + # plt.legend() + # plt.show() + + +# investigate_AATSR4th_transect(plot_BT=True) +# investigate_AATSR4th_transect(plot_Temp=True) + +### cloud shadow processor +cloud_shadow_processor_AASTR(cloudflagType='old') + +# analyse_AATSR4th_transect() +# analyse_AATSR4th_transect(start=11600, end=30860, step=1000) \ No newline at end of file diff --git a/idepix-aatsr/src/main/java/org/esa/snap/idepix/aatsr/IdepixAatsrOp.java b/idepix-aatsr/src/main/java/org/esa/snap/idepix/aatsr/IdepixAatsrOp.java new file mode 100644 index 00000000..b7f382af --- /dev/null +++ b/idepix-aatsr/src/main/java/org/esa/snap/idepix/aatsr/IdepixAatsrOp.java @@ -0,0 +1,734 @@ +/* + * Copyright (C) 2022 Brockmann Consult GmbH (info@brockmann-consult.de) + * + * This program is free software; you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation; either version 3 of the License, or (at your option) + * any later version. + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for + * more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, see http://www.gnu.org/licenses/ + */ +package org.esa.snap.idepix.aatsr; + +import com.bc.ceres.core.ProgressMonitor; +import com.bc.ceres.core.SubProgressMonitor; +import org.esa.snap.core.dataio.geocoding.ComponentGeoCoding; +import org.esa.snap.core.datamodel.Band; +import org.esa.snap.core.datamodel.FlagCoding; +import org.esa.snap.core.datamodel.GeoCoding; +import org.esa.snap.core.datamodel.Mask; +import org.esa.snap.core.datamodel.Product; +import org.esa.snap.core.datamodel.ProductData; +import org.esa.snap.core.datamodel.RasterDataNode; +import org.esa.snap.core.gpf.Operator; +import org.esa.snap.core.gpf.OperatorException; +import org.esa.snap.core.gpf.OperatorSpi; +import org.esa.snap.core.gpf.Tile; +import org.esa.snap.core.gpf.annotations.OperatorMetadata; +import org.esa.snap.core.gpf.annotations.Parameter; +import org.esa.snap.core.gpf.annotations.SourceProduct; +import org.esa.snap.core.gpf.annotations.TargetProduct; +import org.esa.snap.core.image.ImageManager; +import org.esa.snap.core.util.BitSetter; +import org.esa.snap.core.util.ImageUtils; +import org.esa.snap.core.util.ProductUtils; +import org.esa.snap.core.util.SystemUtils; +import org.esa.snap.core.util.math.MathUtils; +import org.esa.snap.idepix.core.IdepixConstants; +import org.esa.snap.idepix.core.IdepixFlagCoding; + +import javax.imageio.ImageIO; +import javax.media.jai.BorderExtender; +import javax.media.jai.ImageLayout; +import javax.media.jai.Interpolation; +import javax.media.jai.JAI; +import javax.media.jai.KernelJAI; +import javax.media.jai.OpImage; +import javax.media.jai.RenderedOp; +import javax.media.jai.operator.ClampDescriptor; +import javax.media.jai.operator.ConvolveDescriptor; +import javax.media.jai.operator.CropDescriptor; +import javax.media.jai.operator.DivideByConstDescriptor; +import javax.media.jai.operator.ExtremaDescriptor; +import javax.media.jai.operator.FormatDescriptor; +import javax.media.jai.operator.MosaicDescriptor; +import javax.media.jai.operator.MultiplyConstDescriptor; +import javax.media.jai.operator.ScaleDescriptor; +import javax.media.jai.operator.SubtractConstDescriptor; +import java.awt.Color; +import java.awt.Dimension; +import java.awt.Rectangle; +import java.awt.RenderingHints; +import java.awt.image.DataBuffer; +import java.awt.image.Raster; +import java.awt.image.RenderedImage; +import java.io.File; +import java.io.IOException; +import java.nio.file.Files; +import java.util.ArrayList; +import java.util.Arrays; +import java.util.List; +import java.util.Locale; +import java.util.concurrent.Callable; +import java.util.concurrent.ExecutionException; +import java.util.concurrent.ExecutorService; +import java.util.concurrent.Executors; +import java.util.concurrent.Future; +import java.util.logging.Level; +import java.util.stream.IntStream; + +/** + * The IdePix pixel classification operator for AATSR products (4th repro). + */ +@OperatorMetadata(alias = "Idepix.Aatsr", + category = "Optical/Preprocessing/Masking", + version = "1.0", + authors = "Dagmar Mueller, Marco Peters", + copyright = "(c) 2022 by Brockmann Consult", + description = "Pixel identification and classification for AATSR 4th repro data.") +public class IdepixAatsrOp extends Operator { + + private static final boolean DEBUG = false; + private final static int SPATIAL_RESOLUTION = 1000; // in meter // better to get it from product + + @SourceProduct(label = "AATSR L1b product", + description = "The AATSR L1b source product.") + private Product sourceProduct; + + @TargetProduct(description = "The target product.") + private Product targetProduct; + + @Parameter(label = "Copy source bands", defaultValue = "false") + private boolean copySourceBands; + + @Parameter(label = "Assumed cloud top height", defaultValue = "6000") + private int cloudTopHeight; + private Mask startSearchMask; + private Rectangle dayTimeROI; + private RenderedOp orientationImage; + private double minSurfaceAltitude; + private int maxShadowDistance; + private Band idepixFlagBand; + + // overall parameters + + @Override + public void initialize() throws OperatorException { + // 1) + validate(sourceProduct); + + // 2) create TargetProduct + final String targetProductName = sourceProduct.getName() + "_idepix"; + targetProduct = createCompatibleProduct(sourceProduct, targetProductName); + + // init ComponentGeoCoding is necessary for SNAP8.x, can be removed for SNAP9 + final GeoCoding sceneGeoCoding = sourceProduct.getSceneGeoCoding(); + if (sceneGeoCoding instanceof ComponentGeoCoding) { + ComponentGeoCoding compGC = (ComponentGeoCoding) sceneGeoCoding; + compGC.initialize(); + } + + if (copySourceBands) { + ProductUtils.copyProductNodes(sourceProduct, targetProduct); + for (Band band : sourceProduct.getBands()) { + if (!targetProduct.containsBand(band.getName())) { + ProductUtils.copyBand(band.getName(), sourceProduct, targetProduct, true); + } + } + } else { + ProductUtils.copyGeoCoding(sourceProduct, targetProduct); + targetProduct.setStartTime(sourceProduct.getStartTime()); + targetProduct.setEndTime(sourceProduct.getEndTime()); + } + + // 2.1) copy source bands (todo - which source bands to include?) + // 2.2) create flag band compatible with other IdePix processors but only + idepixFlagBand = targetProduct.addBand(IdepixConstants.CLASSIF_BAND_NAME, ProductData.TYPE_INT16); + FlagCoding flagCoding = IdepixFlagCoding.createDefaultFlagCoding(IdepixConstants.CLASSIF_BAND_NAME); + targetProduct.getFlagCodingGroup().add(flagCoding); + idepixFlagBand.setSampleCoding(flagCoding); + IdepixFlagCoding.setupDefaultClassifBitmask(targetProduct); + } + + @Override + public void doExecute(ProgressMonitor pm) throws OperatorException { + pm.beginTask("Executing cloud shadow detection...", 10); + try { + final int sceneWidth = sourceProduct.getSceneRasterWidth(); + final int sceneHeight = sourceProduct.getSceneRasterHeight(); + + // 1) detect day time area where cloud shadow can occur + dayTimeROI = getDayTimeArea(sourceProduct); + + // 2) create north-corrected orientation image + orientationImage = computeOrientationImage(sourceProduct); +// writeDebugImage(orientationImage, "orientationImage.png"); + + // 3) create cloudMaskImage and landMaskImage + // as alternative the bayesian_in and confidence_in could be used. See TechNote. + // But currently the bayes_in.no_bayesian_probabilities_available is always set. so it makes no sense to use it. + final Mask cloudMask = Mask.BandMathsType.create("__cloud_mask", "", sceneWidth, sceneHeight, + "cloud_in.visible or cloud_in.12_gross_cloud or cloud_in.11_12_thin_cirrus or cloud_in.3_7_12_medium_high", + Color.white, 0.5f); + cloudMask.setOwner(sourceProduct); + final Mask landMask = Mask.BandMathsType.create("__land_mask", "", sceneWidth, sceneHeight, + "confidence_in.coastline or confidence_in.tidal or confidence_in.land or confidence_in.inland_water", + Color.green, 0.5f); + landMask.setOwner(sourceProduct); + +// writeDebugImage(cloudMask.getSourceImage(), "cloud_mask.png"); +// writeDebugImage(landMask.getSourceImage(), "land_mask.png"); + + // 4) create startSearchMask using cloudMaskImage, landMaskImage and search radius + // splitting cloudMask image into 2000 y-pixel slices but only in dayTimeArea + final List convolvedCloudSlices = new ArrayList<>(); + final List convolvedLandSlices = new ArrayList<>(); + final List slices = sliceRect(dayTimeROI, 2000); + // merge slices +// final Rectangle remove1 = slices.remove(slices.size() - 2); +// final Rectangle remove2 = slices.remove(slices.size() - 1); +// slices.add(remove1.union(remove2)); + + convolveLandAndCloudSlices(slices, dayTimeROI, cloudMask, convolvedCloudSlices, landMask, convolvedLandSlices); + + // mosaic the generated slices + final Rectangle mosaicBounds = new Rectangle(sceneWidth, sceneHeight); + final Rectangle mosaicTileSize = new Rectangle(cloudMask.getSourceImage().getTileWidth(), cloudMask.getSourceImage().getTileHeight()); + final RenderedImage mosaicCloudImage = createMosaic(slices, convolvedCloudSlices, mosaicBounds, mosaicTileSize); + final RenderedImage mosaicLandImage = createMosaic(slices, convolvedLandSlices, mosaicBounds, mosaicTileSize); + writeDebugImage(mosaicCloudImage, "mosaicCloudImage.png"); + writeDebugImage(mosaicLandImage, "mosaicLandImage.png"); + + + // create temporary product to compute the start-search-mask + final Product tempMaskProduct = new Product("temp", "tempType", sceneWidth, sceneHeight); + final Band convCloudMaskBand = new Band("__convCloudMask", ProductData.TYPE_FLOAT32, sceneWidth, sceneHeight); + convCloudMaskBand.setSourceImage(mosaicCloudImage); + tempMaskProduct.addBand(convCloudMaskBand); + final Band convLandMaskBand = new Band("__convLandMask", ProductData.TYPE_FLOAT32, sceneWidth, sceneHeight); + convLandMaskBand.setSourceImage(mosaicLandImage); + tempMaskProduct.addBand(convLandMaskBand); + + // use band math to combine mosaics to startSearchMask + startSearchMask = Mask.BandMathsType.create("__startSearch", "", sceneWidth, sceneHeight, + "__convCloudMask > 0.001 && __convCloudMask < 0.998 && __convLandMask > 0.001", + Color.BLUE, 0.0); + tempMaskProduct.addBand(startSearchMask); + writeDebugImage(startSearchMask.getSourceImage(), "startSearchMask.png"); + + final RasterDataNode elevationRaster = sourceProduct.getRasterDataNode("elevation_in"); + final RenderedOp croppedElev = CropDescriptor.create(elevationRaster.getSourceImage(), (float) dayTimeROI.x, (float) dayTimeROI.y, (float) dayTimeROI.width, (float) dayTimeROI.height, null); + final RenderedOp clampedElev = ClampDescriptor.create(croppedElev, new double[]{0}, new double[]{Short.MAX_VALUE}, null); + final RenderedOp extrema = ExtremaDescriptor.create(clampedElev, null, 10, 10, Boolean.FALSE, 1, null); + minSurfaceAltitude = ((double[]) extrema.getProperty("minimum"))[0]; + + pm.worked(1); + + doShadowDetectionPerSlice(cloudMask, landMask, slices, new SubProgressMonitor(pm, 9)); + // force creation of source image to prevent creation of new image by GPF + idepixFlagBand.getSourceImage(); + } catch (IOException e) { + throw new OperatorException("Could not read source data", e); + } finally { + pm.done(); + } + } + + private void doShadowDetectionPerSlice(Mask cloudMask, Mask landMask, List slices, ProgressMonitor pm) { + final int shadowValue = BitSetter.setFlag(0, IdepixConstants.IDEPIX_CLOUD_SHADOW); + idepixFlagBand.setRasterData(idepixFlagBand.createCompatibleRasterData()); + + final ExecutorService executorService = Executors.newFixedThreadPool((int) (Runtime.getRuntime().availableProcessors() * 0.8)); // Use 80% use cores; is this good? + + final Rectangle extendedDayTimeRoi = (Rectangle) dayTimeROI.clone(); + extendedDayTimeRoi.grow(0, maxShadowDistance); + final Tile elevation = getSourceTile(sourceProduct.getRasterDataNode("elevation_in"), extendedDayTimeRoi); + final Tile cloudMaskData = getSourceTile(cloudMask, extendedDayTimeRoi); + + final ArrayList> tasks = new ArrayList<>(); + for (Rectangle slice : slices) { + tasks.add(() -> { + final Tile sza = getSourceTile(sourceProduct.getRasterDataNode("solar_zenith_tn"), slice); + final Tile saa = getSourceTile(sourceProduct.getRasterDataNode("solar_azimuth_tn"), slice); + final Tile oza = getSourceTile(sourceProduct.getRasterDataNode("sat_zenith_tn"), slice); + final Tile x_tx = getSourceTile(sourceProduct.getRasterDataNode("x_tx"), slice); + final Raster orientation = orientationImage.getData(slice); + final Tile landMaskData = getSourceTile(landMask, slice); + final Raster startData = ((RenderedImage) startSearchMask.getSourceImage()).getData(slice); + + for (int i = slice.y; i < slice.y + slice.height; ++i) { + for (int j = slice.x; j < slice.x + slice.width; ++j) { + if (startData.getSample(j, i, 0) > 0 && cloudMaskData.getSampleInt(j, i) > 0) { + final PathAndHeightInfo pathAndHeightInfo = calcPathAndTheoreticalHeight(sza.getSampleFloat(j, i), saa.getSampleFloat(j, i), + oza.getSampleFloat(j, i), x_tx.getSampleFloat(j, i), + orientation.getSampleFloat(j, i, 0), + SPATIAL_RESOLUTION, + cloudTopHeight, + minSurfaceAltitude); + + final int[][] indexArray = pathAndHeightInfo.illuPathSteps.clone(); + final int shiftJ = j; + final int shiftI = i; + IntStream.range(0, indexArray.length).parallel().forEach(n -> {indexArray[n][0]+=shiftJ;indexArray[n][1]+=shiftI;}); + + // find cloud free positions along the search path: + final Boolean[] id = Arrays.stream(indexArray).parallel().map( + ints -> ints[0] >= 0 && ints[0] < slice.x + dayTimeROI.width && + ints[1] >= 0 && ints[1] < slice.y + dayTimeROI.height).toArray(Boolean[]::new); + + final int sum = Arrays.stream(id).parallel().mapToInt(aBoolean -> aBoolean ? 1 : 0).sum(); + if (sum > 3) { // path positions + + final int[][] curIndexArray = IntStream.range(0, indexArray.length).parallel().filter(n -> id[n]).mapToObj(n -> indexArray[n]).toArray(int[][]::new); + final double[] illuPathHeight = pathAndHeightInfo.illuPathHeight; + final double[] baseHeightArray = IntStream.range(0, curIndexArray.length).parallel().filter(n -> id[n]).mapToDouble(n -> illuPathHeight[n]).toArray(); + final double[] elevPath = Arrays.stream(curIndexArray).parallel().mapToDouble(index -> elevation.getSampleFloat(index[0], index[1])).toArray(); + final Boolean[] cloudPath = Arrays.stream(curIndexArray).parallel().map(index -> cloudMaskData.getSampleInt(index[0], index[1]) > 0).toArray(Boolean[]::new); + + final Boolean[] id2 = IntStream.range(0, baseHeightArray.length).parallel().mapToObj(value -> (Math.abs(baseHeightArray[value] - elevPath[value]) < pathAndHeightInfo.threshHeight) && !cloudPath[value]).toArray(Boolean[]::new); + + IntStream.range(0, baseHeightArray.length).parallel().filter(n -> id2[n]).forEach( + n -> idepixFlagBand.setPixelInt(curIndexArray[n][0], curIndexArray[n][1], shadowValue)); + } + } + + int flagValue = idepixFlagBand.getPixelInt(j, i); + if (cloudMaskData.getSampleInt(j, i) > 0) { + flagValue = BitSetter.setFlag(flagValue, IdepixConstants.IDEPIX_CLOUD); + } + if (landMaskData.getSampleInt(j, i) > 0) { + flagValue = BitSetter.setFlag(flagValue, IdepixConstants.IDEPIX_LAND); + } + + idepixFlagBand.setPixelInt(j, i, flagValue); + } + } + return slice; + }); + } + + try { + final List> futures = new ArrayList<>(); + for (Callable task : tasks) { + futures.add(executorService.submit(task)); + } + + pm.beginTask("Detecting clouds shadows...", futures.size()); + try { + executorService.shutdown(); + int stillRunning = futures.size(); + while (!executorService.isTerminated()) { + synchronized (cloudMask) { + cloudMask.wait(100); + } + + int numRunning = numRunningFutures(futures); + pm.worked(stillRunning - numRunning); + stillRunning = numRunning; + } + } finally { + pm.done(); + } + } catch (InterruptedException e) { + throw new OperatorException(e); + } + + } + + private int numRunningFutures(List> futures) throws InterruptedException { + int numRunning = 0; + for (int i = 0; i < futures.size(); i++) { + Future future = futures.get(i); + numRunning = numRunning + (future.isDone() ? 0 : 1); + if (future.isDone()) { + final Future doneFuture = futures.remove(i); + try { + doneFuture.get(); + } catch (ExecutionException e) { + throw new OperatorException("Error during calculation of cloud shadow", e); + } + } + } + return numRunning; + } + + @Override + public void computeTile(Band targetBand, Tile targetTile, ProgressMonitor pm) throws OperatorException { + // We need to place the data into the targetTile even though it is already in the raster. + // Otherwise, the data is not written to the product by GPF. + final Rectangle rect = targetTile.getRectangle(); + final int[] pixels = new int[rect.width * rect.height]; + final int[] samplesInt = idepixFlagBand.getPixels(rect.x, rect.y, rect.width, rect.height, pixels); + targetTile.setSamples(samplesInt); + } + + @SuppressWarnings("SameParameterValue") + static PathAndHeightInfo calcPathAndTheoreticalHeight(double sza, double saa, double oza, double xtx, double orientation, int spatialRes, int maxObjectAltitude, double minSurfaceAltitude) { + sza = correctSzaForOzaInfluence(sza, saa, oza, xtx, orientation); + + final double shadowAngle = ((saa - orientation) + 180) * Math.PI / 180 - Math.PI / 2; + double cosSaa = Math.cos(shadowAngle); + double sinSaa = Math.sin(shadowAngle); + + double deltaProjX = ((maxObjectAltitude - minSurfaceAltitude) * Math.tan(sza * Math.PI / 180) * cosSaa) / spatialRes; + double deltaProjY = ((maxObjectAltitude - minSurfaceAltitude) * Math.tan(sza * Math.PI / 180) * sinSaa) / spatialRes; + + double x0 = 0; + double y0 = 0; + double x1 = x0 + deltaProjX + Math.signum(deltaProjX) * 1.5; + double y1 = y0 + deltaProjY + Math.signum(deltaProjY) * 1.5; + + // create index steps + // Path touches which pixels? + // setup all pixel centers from x0/y0 to x1/y1. + // calculate distance between pixel center and line (X0, X1) + // all distances below/equal sqrt(2*0.5^2): the pixel is touched by the line and a potential shadow pixel. + double[] xCenters; + double[] yCenters; + if (x0 < x1) { + xCenters = arange(x0 + 0.5, Math.round(x1) + 0.5, 1.0); + } else { + xCenters = arange(Math.round(x1) + 0.5, x0 + 0.5, 1.0); + } + if (y0 < y1) { + yCenters = arange(y0 + 0.5, Math.round(y1) + 0.5, 1.0); + } else { + yCenters = arange(Math.round(y1) + 0.5, y0 + 0.5, 1.0); + } + + double[][] distance = new double[xCenters.length][yCenters.length]; + int[][] xIndex = new int[xCenters.length][yCenters.length]; + int[][] yIndex = new int[xCenters.length][yCenters.length]; + int nxCenter = xCenters.length; + int nyCenter = yCenters.length; + + final double divider = Math.pow(x0 - x1, 2) + Math.pow(y0 - y1, 2); + for (int i = 0; i < nxCenter; i++) { + for (int j = 0; j < nyCenter; j++) { + double r = -((x0 - xCenters[i]) * (x0 - x1) + (y0 - yCenters[j]) * (y0 - y1)) / divider; + double d = Math.sqrt(Math.pow((x0 - xCenters[i]) + r * (x0 - x1), 2) + Math.pow((y0 - yCenters[j]) + r * (y0 - y1), 2)); + distance[i][j] = d; + xIndex[i][j] = deltaProjX < 0 ? i - (nxCenter - 1) : i; + yIndex[i][j] = deltaProjY < 0 ? j - (nyCenter - 1) : j; + } + } + double halfPixelDistance = 0.5 * Math.sqrt(2); + int[][] id = new int[xCenters.length][yCenters.length]; + for (int x = 0; x < id.length; x++) { + for (int y = 0; y < id[x].length; y++) { + id[x][y] = distance[x][y] <= halfPixelDistance ? 1 : 0; + } + } + final int numCoords = Arrays.stream(id).flatMapToInt(Arrays::stream).sum(); + int[][] stepIndex = new int[numCoords][2]; + int s = 0; + for (int i = 0; i < id.length; i++) { + for (int j = 0; j < id[i].length; j++) { + if (id[i][j] == 1) { + stepIndex[s++] = new int[]{xIndex[i][j], yIndex[i][j]}; + } + } + } + + double[] theoretHeight = new double[stepIndex.length]; + final double tanSza = Math.tan(sza * Math.PI / 180.); + for (int i = 0; i < theoretHeight.length; i++) { + theoretHeight[i] = maxObjectAltitude - Math.sqrt(Math.pow(stepIndex[i][0], 2) + Math.pow(stepIndex[i][1], 2)) * spatialRes / tanSza; + } + double threshHeight = 1000. / tanSza; + + return new PathAndHeightInfo(stepIndex, theoretHeight, threshHeight); + } + + @SuppressWarnings("SameParameterValue") + static double[] arange(double startValue, double endValue, double step) { + final int numValues = (int) (Math.ceil((endValue - startValue) / step)); + return IntStream.range(0, numValues).mapToDouble(x -> x * step + startValue).toArray(); + } + + private static double correctSzaForOzaInfluence(double sza, double saa, double oza, double xtx, double orientation) { + if (saa - orientation < 180) { + sza = xtx < 0 ? correctSzaNegative(sza, oza) : correctSzaPositive(sza, oza); + } else { + sza = xtx < 0 ? correctSzaPositive(sza, oza) : correctSzaNegative(sza, oza); + } + return sza; + } + + private static double correctSzaPositive(double sza, double oza) { + return Math.atan(Math.tan(sza * Math.PI / 180) + Math.tan(oza * Math.PI / 180)) * 180 / Math.PI; + } + + private static double correctSzaNegative(double sza, double oza) { + return Math.atan(Math.tan(sza * Math.PI / 180) - Math.tan(oza * Math.PI / 180)) * 180 / Math.PI; + } + + void validate(Product sourceProduct) throws OperatorException { + if (!"ENV_AT_1_RBT".equals(sourceProduct.getProductType())) { + throw new OperatorException("An AATSR product from the 4th reprocessing is needed as input"); + } + String[] usedRasterNames = {"solar_zenith_tn", "solar_azimuth_tn", "sat_zenith_tn", "latitude_tx", "longitude_tx", "x_tx", "elevation_in", "confidence_in", "cloud_in"}; + for (String usedRasterName : usedRasterNames) { + if (!sourceProduct.containsRasterDataNode(usedRasterName)) { + throw new OperatorException(String.format("Missing raster '%s' in source product", usedRasterName)); + } + } + } + + /** + * Returns the daytime area of the scene and where the sun angles are usable + * + * @param scene the scene to compute the daytime area + * @return a rectangle specifying the daytime area with usable sun angles + */ + private Rectangle getDayTimeArea(Product scene) { + final int sceneWidth = scene.getSceneRasterWidth(); + final int sceneHeight = scene.getSceneRasterHeight(); + + // test sun elevation and find first and last row with SZA<85°, daylight zone + // day flag is necessary because there can be an area with SZA<85° but it is not marked as DAY. + final Mask szaMask = Mask.BandMathsType.create("__SZA_DAY_Mask", "", sceneWidth, sceneHeight, + "solar_zenith_tn * (confidence_in.day ? 1 : NaN) < 85", + Color.yellow, 0.5f); + szaMask.setOwner(scene); + final int[] szaRange0 = detectMaskedPixelRangeInColumn(szaMask, 0); + final int[] szaRange1 = detectMaskedPixelRangeInColumn(szaMask, sceneWidth - 1); + // free used memory + szaMask.setOwner(null); + szaMask.dispose(); + + int szaUpperLimit = Math.max((szaRange0[0]), szaRange1[0]); + int szaLowerLimit = Math.min((szaRange0[1]), szaRange1[1]); + + return new Rectangle(0, szaUpperLimit, sceneWidth, szaLowerLimit - szaUpperLimit); + } + + /** + * Return the min and maximum line number in the where pixels are masked. + * + * @param mask the mask + * @param column the column to check masked values + * @return an integer array containing the minimum (index=0) and maximum (index=1) + */ + static int[] detectMaskedPixelRangeInColumn(Mask mask, int column) { + final int rasterHeight = mask.getRasterHeight(); + final Raster data = mask.getSourceImage().getData(new Rectangle(column, 0, 1, rasterHeight)); + // masks have byte data type + final byte[] dataBufferArray = (byte[]) ImageUtils.createDataBufferArray(data.getTransferType(), rasterHeight); + data.getDataElements(column, 0, 1, rasterHeight, dataBufferArray); + int[] minMax = new int[2]; + // todo - can be further optimised and parallelized + for (int i = 0; i < dataBufferArray.length; i++) { + if (dataBufferArray[i] == -1) { + minMax[0] = i; + break; + } + } + for (int i = dataBufferArray.length - 1; i >= 0; i--) { + if (dataBufferArray[i] == -1) { + minMax[1] = i; + break; + } + } + return minMax; + } + + /** + * Creates a north-corrected orientation image on a sub-sampled grid. The result is scale back to the original size + * using bicubic interpolation. The computation is based on 'latitude_tx' and 'longitude_tx'. + * + * @param scene the scene to compute the north-corrected orientation image for. + * @return the north-corrected orientation image + */ + private RenderedOp computeOrientationImage(Product scene) { + final RasterDataNode latRaster = scene.getRasterDataNode("latitude_tx"); // todo - better use latitude_in, but this gives strange results + final RasterDataNode lonRaster = scene.getRasterDataNode("longitude_tx"); // todo - better use longitude_in, but this gives strange results + final Interpolation nearest = Interpolation.getInstance(Interpolation.INTERP_NEAREST); + final float downScaleFactor = 1 / 10F; + final RenderedOp lowResLats = ScaleDescriptor.create(latRaster.getSourceImage(), downScaleFactor, downScaleFactor, 0.0F, 0.0F, nearest, null); + final RenderedOp lowResLons = ScaleDescriptor.create(lonRaster.getSourceImage(), downScaleFactor, downScaleFactor, 0.0F, 0.0F, nearest, null); + final OpImage lowResOrientation = new OrientationOpImage(lowResLats, lowResLons); + final float upScaleWidthFactor = scene.getSceneRasterWidth() / (float) lowResOrientation.getWidth(); + final float upScaleHeightFactor = scene.getSceneRasterHeight() / (float) lowResOrientation.getHeight(); + return ScaleDescriptor.create(lowResOrientation, upScaleWidthFactor, upScaleHeightFactor, 0.0F, 0.0F, + Interpolation.getInstance(Interpolation.INTERP_BICUBIC), + new RenderingHints(JAI.KEY_BORDER_EXTENDER, BorderExtender.createInstance(BorderExtender.BORDER_COPY))); + } + + private RenderedImage createMosaic(List slices, List imageSlices, Rectangle mosaicBounds, Rectangle mosaicTileSize) { + double[][] sourceThresholds = new double[slices.size()][1]; + for (int i = 0; i < sourceThresholds.length; i++) { + sourceThresholds[i][0] = Double.MIN_VALUE; + } + final ImageLayout mosaicImageLayout = ImageManager.createSingleBandedImageLayout(DataBuffer.TYPE_FLOAT, + mosaicBounds.width, + mosaicBounds.height, + mosaicTileSize.width, mosaicTileSize.height); + return MosaicDescriptor.create(imageSlices.toArray(new RenderedImage[0]), + MosaicDescriptor.MOSAIC_TYPE_OVERLAY, + null, + null, + sourceThresholds, + new double[]{0.0},//backgroundValue, + new RenderingHints(JAI.KEY_IMAGE_LAYOUT, mosaicImageLayout)); + } + + private void convolveLandAndCloudSlices(List slices, Rectangle dayTimeROI, Mask cloudMask, List convolvedCloudSlices, Mask landMask, List convolvedLandSlices) throws IOException { + final RasterDataNode sza = sourceProduct.getRasterDataNode("solar_zenith_tn"); + // convolution shall take place with float data type. We need to format the source images of cloudMask and landMask. + final RenderedImage floatCloudMaskImage = FormatDescriptor.create(cloudMask.getSourceImage(), DataBuffer.TYPE_FLOAT, null); + final RenderedImage floatLandMaskImage = FormatDescriptor.create(landMask.getSourceImage(), DataBuffer.TYPE_FLOAT, null); + // convolve + for (Rectangle slice : slices) { + if (slice.intersects(dayTimeROI)) { + // only convolve slices intersecting with dayTimeROI. Areas outside are handled by default background value when creating the mosaic. + double radius = computeKernelRadiusForSlice(sza, slice); + maxShadowDistance = MathUtils.ceilInt(Math.max(radius, maxShadowDistance)); + final KernelJAI jaiKernel = createJaiKernel(radius, new Dimension(1000, 1000)); + convolveSlice(floatCloudMaskImage, slice, jaiKernel, convolvedCloudSlices, "convCloudImage"); + convolveSlice(floatLandMaskImage, slice, jaiKernel, convolvedLandSlices, "convLandImage"); + } + } + } + + @SuppressWarnings("unused") + private void convolveSlice(RenderedImage sourceImage, Rectangle slice, KernelJAI jaiKernel, List convolvedSlices, String debugId) { + final RenderedOp curSlice = CropDescriptor.create(sourceImage, (float) slice.x, (float) slice.y, (float) slice.width, (float) slice.height, null); + final RenderedImage convImage = createConvolvedImage(curSlice, jaiKernel); + convolvedSlices.add(convImage); + writeDebugImage(convImage, String.format("%s_%d_%dx%d.png", debugId, slice.y, jaiKernel.getWidth(), jaiKernel.getWidth())); + } + + private double computeKernelRadiusForSlice(RasterDataNode sza, Rectangle slice) throws IOException { + final float[] szaData = new float[slice.width * slice.height]; + sza.readPixels(slice.x, slice.y, slice.width, slice.height, szaData); + Arrays.sort(szaData); + float szaMedian; + if (szaData.length % 2 == 0) { + szaMedian = (szaData[szaData.length / 2] + szaData[szaData.length / 2 - 1]) / 2; + } else { + szaMedian = szaData[szaData.length / 2]; + } + return Math.abs(cloudTopHeight * Math.tan(szaMedian * Math.PI / 180.0)); + } + + @SuppressWarnings("SameParameterValue") + static List sliceRect(Rectangle sourceBounds, int sliceHeight) { + List slices = new ArrayList<>(); + for (int i = sourceBounds.y; i < sourceBounds.y + sourceBounds.height; i += sliceHeight) { + Rectangle curRect = new Rectangle(0, i, sourceBounds.width, sliceHeight); + curRect = sourceBounds.intersection(curRect); + if (!curRect.isEmpty()) { + slices.add(curRect); + } + } + return slices; + } + + private RenderedImage createConvolvedImage(RenderedImage sourceImage, KernelJAI jaiKernel) { + final RenderedOp convImage = ConvolveDescriptor.create(sourceImage, jaiKernel, new RenderingHints(JAI.KEY_BORDER_EXTENDER, BorderExtender.createInstance(BorderExtender.BORDER_COPY))); + // due to convolution value can be higher than 255, so we clamp + final RenderedOp clampImage = ClampDescriptor.create(convImage, new double[]{0}, new double[]{255}, null); + // normalise to [0,1.0] value range + return DivideByConstDescriptor.create(clampImage, new double[]{255}, null); + } + + private KernelJAI createJaiKernel(double radius, Dimension spacing) { + Dimension kernelHalfDim = new Dimension(MathUtils.ceilInt(radius / spacing.width), + MathUtils.ceilInt(radius / spacing.height)); + int[] xDist = createDistanceArray(kernelHalfDim.width, spacing.width); + int[] yDist = createDistanceArray(kernelHalfDim.height, spacing.height); + return createKernelData(xDist, yDist, radius); + } + + static int[] createDistanceArray(int kernelHalfDim, int spacing) { + int[] xDist = new int[kernelHalfDim * 2 + 1]; + for (int i = 0; i < xDist.length; i++) { + xDist[i] = -kernelHalfDim * spacing + i * spacing; + } + return xDist; + } + + private KernelJAI createKernelData(final int[] xDist, final int[] yDist, double radius) { + final double[][] kernel = new double[yDist.length][xDist.length]; + for (int y = 0; y < kernel.length; y++) { + double[] rowData = kernel[y]; + for (int x = 0; x < rowData.length; x++) { + rowData[x] = Math.sqrt(Math.pow(yDist[y], 2) + Math.pow(xDist[x], 2)) <= radius ? 1.0 : 0.0; + } + } + final double kernelSum = Arrays.stream(kernel).flatMapToDouble(Arrays::stream).sum(); + for (double[] rowData : kernel) { + for (int x = 0; x < rowData.length; x++) { + rowData[x] = rowData[x] / kernelSum; + } + } + final double[] oneDimKernel = Arrays.stream(kernel).flatMapToDouble(Arrays::stream).toArray(); + final float[] oneDimFloatKernel = new float[oneDimKernel.length]; + for (int i = 0; i < oneDimKernel.length; i++) { + oneDimFloatKernel[i] = (float) oneDimKernel[i]; + } + + return new KernelJAI(xDist.length, yDist.length, oneDimFloatKernel); + } + + private Product createCompatibleProduct(Product sourceProduct, String name) { + int sceneWidth = sourceProduct.getSceneRasterWidth(); + int sceneHeight = sourceProduct.getSceneRasterHeight(); + return new Product(name, "AATSR_IDEPIX", sceneWidth, sceneHeight); + } + + @SuppressWarnings("might be used for debug purpose") + private void writeDebugImage(RenderedImage image, String filename) { + if (DEBUG) { + final File outputDir = new File("target/images"); + final File output = new File(outputDir, filename); + try { + Files.createDirectories(output.toPath().getParent()); + if (!ImageIO.write(image, "PNG", output)) { + SystemUtils.LOG.log(Level.WARNING, "No writer found for image '" + filename + "', trying to reformat the image"); + final RenderedOp extrema = ExtremaDescriptor.create(image, null, 10, 10, Boolean.FALSE, 1, null); + final double[] minimum = (double[]) extrema.getProperty("minimum"); + final double[] maximum = (double[]) extrema.getProperty("maximum"); + final RenderedOp step1 = SubtractConstDescriptor.create(image, minimum, null); + final RenderedOp normImage = DivideByConstDescriptor.create(step1, new double[]{maximum[0] - minimum[0]}, null); + final RenderedOp scaledImage = MultiplyConstDescriptor.create(normImage, new double[]{255}, null); + final RenderedOp formattedImage = FormatDescriptor.create(scaledImage, DataBuffer.TYPE_BYTE, null); + if (!ImageIO.write(formattedImage, "PNG", output)) { + SystemUtils.LOG.log(Level.WARNING, "Retry of writing if image '" + filename + "'did not work too. Giving up!"); + } + } + } catch (IOException e) { + throw new RuntimeException(e); + } + } + } + + /** + * The Service Provider Interface (SPI) for the operator. + * It provides operator meta-data and is a factory for new operator instances. + */ + public static class Spi extends OperatorSpi { + + public Spi() { + super(IdepixAatsrOp.class); + } + } + + static class PathAndHeightInfo { + final int[][] illuPathSteps; + final double[] illuPathHeight; + final double threshHeight; + + public PathAndHeightInfo(int[][] stepIndex, double[] theoretHeight, double threshHeight) { + this.illuPathSteps = stepIndex; + this.illuPathHeight = theoretHeight; + this.threshHeight = threshHeight; + } + } +} diff --git a/idepix-aatsr/src/main/java/org/esa/snap/idepix/aatsr/OrientationOpImage.java b/idepix-aatsr/src/main/java/org/esa/snap/idepix/aatsr/OrientationOpImage.java new file mode 100644 index 00000000..7c0da415 --- /dev/null +++ b/idepix-aatsr/src/main/java/org/esa/snap/idepix/aatsr/OrientationOpImage.java @@ -0,0 +1,79 @@ +/* + * Copyright (c) 2022. Brockmann Consult GmbH (info@brockmann-consult.de) + * + * This program is free software; you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation; either version 3 of the License, or (at your option) + * any later version. + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for + * more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, see http://www.gnu.org/licenses/ + * + */ + +package org.esa.snap.idepix.aatsr; + +import org.esa.snap.core.image.ImageManager; +import org.esa.snap.core.util.math.MathUtils; + +import javax.media.jai.OpImage; +import javax.media.jai.PlanarImage; +import java.awt.Rectangle; +import java.awt.image.DataBuffer; +import java.awt.image.Raster; +import java.awt.image.RenderedImage; +import java.awt.image.WritableRaster; + +/** + * @author Marco Peters + */ +public class OrientationOpImage extends OpImage { + + public OrientationOpImage(RenderedImage lats, RenderedImage lons) { + super(vectorize(lats, lons), ImageManager.createSingleBandedImageLayout(DataBuffer.TYPE_DOUBLE, lats.getWidth(), lats.getHeight(), lats.getTileWidth(), lats.getTileHeight()), + null, false); + } + + @Override + protected void computeRect(PlanarImage[] planarImages, WritableRaster writableRaster, Rectangle destRect) { + PlanarImage lats = planarImages[0]; + PlanarImage lons = planarImages[1]; + int x0 = destRect.x; + int y0 = destRect.y; + final int x1 = x0 + writableRaster.getWidth()-1; + final int y1 = y0 + writableRaster.getHeight()-1; + final Raster latsData = lats.getData(destRect); + final Raster lonsData = lons.getData(destRect); + for (int y = y0; y <= y1; y++) { + for (int x = x0; x <= x1; x++) { + final int x_a = MathUtils.crop(x - 1, x0, x1); + final int x_b = MathUtils.crop(x + 1, x0, x1); + final float lat1 = latsData.getSampleFloat(x_a, y, 0); + final float lon1 = lonsData.getSampleFloat(x_a, y, 0); + final float lat2 = latsData.getSampleFloat(x_b, y, 0); + final float lon2 = lonsData.getSampleFloat(x_b, y, 0); + final double orientation = computeOrientation(lat1, lon1, lat2, lon2); + writableRaster.setSample(x, y, 0, orientation); + } + } + } + + static double computeOrientation(float lat1, float lon1, float lat2, float lon2) { + return Math.atan2(-(lat2 - lat1), (lon2 - lon1) * Math.cos(lat1 * Math.PI / 180.0)) * 180.0 / Math.PI; + } + + @Override + public Rectangle mapSourceRect(Rectangle rectangle, int i) { + return new Rectangle(rectangle); + } + + @Override + public Rectangle mapDestRect(Rectangle rectangle, int i) { + return new Rectangle(rectangle); + } + +} diff --git a/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/AatsrAlgo.html b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/AatsrAlgo.html new file mode 100644 index 00000000..97528e78 --- /dev/null +++ b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/AatsrAlgo.html @@ -0,0 +1,141 @@ + + + + + SNAP Data Processors - IdePix AATSR Cloud Shadow + + + + + + + + + + +
  + IdePix AATSR Cloud Shadow Algorithm + +
+ +

Algorithm Description

+

The cloud shadow algorithm is based on geometrical considerations alone. With the cloud top height and the + illumination direction the shadow position can be determined. As elevated objects are mapped in their apparent + position on the ground if not viewed directly from nadir, the sun and observation angles have to be transformed + accordingly. In that way, the algorithm translates the question of shadows from spherical geometry of the Earth + observation and geo-positions into the space of the regular projection grid.

+ +

North Direction

+

As all calculations of the shadow are translated to the geometry of the pixel grid, it is necessary to calculate the + north direction (also called orientation or bearing) for each pixel individually. + The orientation for a pixel (i,j) is derived from the neighboring pixel (i, j-1) and (i, j+1) from pixel-geocoded + location:

+

+

Defining a cloud mask

+The following cloud flags have been combined in a cloud mask:
+Cloud_in.visible (bit 0), cloud_in.gross_cloud (bit 7), cloud_in.thin_cirrus (bit 8) and the cloud_in.medium_high (bit +9). +

Adjustment of Sun Zenith Angle for Elevated Objects

+Under tilted view (view zenith angle > 0°) elevated objects of unknown height like clouds are projected along the line +of view on the surface, so that their apparent location differs from the actual position over ground (nadir view). +If view and sun azimuth angle are positioned in the same halfspace, both left or right of the nadir line (center line of +grid) in the projected grid, x_tx>0 (VAA*180°), SAA*>180° or x_tx<0 + (VAA*<180°), +SAA*<180° (angles are corrected by +orientation), the geometry of the apparent sun zenith angle which causes the shadow position can be described in the +following way: +From +

+follows +

+

For view and sun direction in opposite directions (x_tx>0 (VAA*>180°), SAA*<180° + or + x_tx<0 (VAA*<180°), SAA*>180°, angles corrected for orientation) + follows accordingly:

+

+ +

The sun azimuth angles are corrected by the orientation (North direction) at the current cloud pixel, so that they + represent the azimuth angle in the projected grid coordinates with 0° in upwards direction on the grid. The view + azimuth angle is replaced by the tie point grid x_tx, which gives the distance from the nadir line at the center + position of a pixel. It changes its sign from left of nadir x_tx <0 to right of nadir x_tx>0. This is the easiest + way to find the viewing direction (without interpolation and corrections), and it allows the algorithm to process + subsets of the swath width.

+ + + + + + + + + +
Geometrical correction for apparent sun zenith angle (VAA>180°, SAA<180° or VAA<180°, + SAA>180°) + Geometrical correction for apparent sun zenith angle (VAA>180°, SAA>180° or VAA<180°, + SAA<180°) +
+ +

Determining the search path in illumination direction

+

Starting from a cloud pixel, which is defined by the cloud flag expression, the illumination path is projected on the + grid and all pixels up to a maximum distance are identified which are intersected by this path.

+

With the adjusted sun zenith angle θ*S and the azimuth angles adjusted for North + direction, so that they represent + the azimuth on the grid against the Y-direction, the geometry of the illumination path on the projection grid can be + fully described.

+

Orientation (North direction) at pixel [i, j] is calculated by the positions at neighboring pixels
+ p[i, j-1] = (lon1, lat1) and p[i, j+1]= (lon2, lat2).

+

+

The theoretical maximum length of the projected path in grid coordinates is defined by the range of surface elevation + and the adjusted angles in x and y direction (spatial resolution 1km in AATSR products):

+

+

The relative grid coordinates of the start point at the cloud pixel at set to (x0, y0) = (0, 0).

+

As the spatial resolution of the grid is quite coarse with regard to the expected cloud top heights, which are + currently fixed at 6km, the maximum extent of the search path is defined as

+

+

With sign(x) defined as: for x>=0, sign(x)=1; for x<0, sign(x)=-1.
+ For all integer value combinations between (x0,y0) and (x1,y1), the center position of the + pixel is set to relative grid coordinates + (0.5, 0.5). The distance from each pixel center to the + theoretical line of illumination is calculated. If the distance is smaller than 0.5*√2, the pixel + area is intersected by the line, and the pixel is potentially a shadow pixel.

+

+

For the intersected pixels (distance < 0.5*√2, relative grid coordinates x, y) the theoretical + height of the illumination path is calculated as:

+

+

Where the theoretical line intersects the surface defined by its altitude, the shadow can be found. The discreet + nature of the grid calls for a height threshold, so that shadow pixels can be identified:

+

+

The cloud shadow flag is raised for all the pixels of the path, which are not masked as cloud and the theoretical + height intersects with the surface elevation:

+

+ +

Algorithmic Limitations

+

Cloud Top Height

+

The AATSR Level-1 product does not provide cloud top height information. Therefor this values is predefined to a + value of 6000 meter. The processor interface allows to adjust this value for each processed scene.

+ +

Limited to Land

+

The detection of cloud shadow is limited to pixels marked as land. Above water no shadow detection is performed.

+ +

Limited to Daytime

+

The cloud shadow calculation is limited to the daytime. As daytime is the part of the orgbit considered where SZA<85° + and the confidence_in flag day (bit value 10) is raised.

+
+
+ + diff --git a/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/AatsrIntro.html b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/AatsrIntro.html new file mode 100644 index 00000000..59c252be --- /dev/null +++ b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/AatsrIntro.html @@ -0,0 +1,67 @@ + + + + + SNAP Data Processors - IdePix - AATSR Cloud Shadow + + + + + + + + + + +
  + IdePix AATSR Cloud Shadow Intro + +
+ +

Overview

+ +

The nadir product of AATSR 4th reprocessing is missing a flag indicating cloud shadows. The IdePix AATSR Cloud Shadow + processor aims to provide such a flag over land. + It has been the goal to calculate this cloud shadow flag based on the information in the Level-1 product only. From + the + cloud flag and the geometry of sun illumination and observation almost all necessary information is given, which + allows a translation of the cloud shadow algorithm for OLCI to AATSR. The information which is missing in the + Level-1 product, is the cloud top height value. Thus this value can be provided as a parameter which is then used + for the whole scene. + The general idea of the shadow detection algorithm is to start from a pixel, which has been flagged as a cloud, and + follow the direction of illumination until the surface is reached. + This Idepix processor does not provide all flags which are defined by other IdePix processors. This one sets only + 3 of the general IdePix flags.
+ These are:
+ IDEPIX_CLOUD, IDEPIX_CLOUD_SHADOW and IDEPIX_LAND +

+ +

Further details are provided

+ +
+ +
+
+ + diff --git a/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/AatsrProcessor.html b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/AatsrProcessor.html new file mode 100644 index 00000000..b05ee881 --- /dev/null +++ b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/AatsrProcessor.html @@ -0,0 +1,91 @@ + + + + + SNAP Data Processors - IdePix AATSR Cloud Shadow + + + + + + + + + + +
  + IdePix AATSR Cloud Shadow Processor + +
+ +

Processor Description

+ +

I/O Parameters Tab

+ +

+ +
Source Product Group
+ +

+ Name: + Used to select the spectral source product. The source product shall + contain spectral bands providing a source spectrum at each pixel. Use the ... button to + open a data product currently not opened in the Sentinel Toolbox.
+ Supported Source Products: The processor supports AATSR L1 product from the 4th reprocessing in the Safe + format or converted to e.g. BEAM-DIMAP. +

+ +
Target Product Group
+ +

+ Name: + Used to specify the name of the target product. +

+ +

+ Save to: + Used to specify whether the target product should be saved to the file system. The + combo box presents a list of file formats. +

+ +

+ Open in SNAP: + Used to specify whether the target product should be opened in the Sentinel Toolbox. + When the target product is not saved, it is opened in the Sentinel Toolbox automatically. +

+ +

The Processing Parameters

+ +

+ +

+ Copy source bands:
+ If enabled the bands of the source product are copied to the target product. If disabled, the target will only contain the flag band. +

+ +

+ Assumed cloud top height:
+ This defines the cloud top height value used by the algorithm. The default is 6000 meter. +

+ + +
+
+ + diff --git a/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/Aatsr_IO-Params.png b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/Aatsr_IO-Params.png new file mode 100644 index 00000000..4f2857a9 Binary files /dev/null and b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/Aatsr_IO-Params.png differ diff --git a/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/Aatsr_Proc-Params.png b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/Aatsr_Proc-Params.png new file mode 100644 index 00000000..e6e1dbaa Binary files /dev/null and b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/Aatsr_Proc-Params.png differ diff --git a/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/IlluminationPath_geometry1_small.png b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/IlluminationPath_geometry1_small.png new file mode 100644 index 00000000..3933454f Binary files /dev/null and b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/IlluminationPath_geometry1_small.png differ diff --git a/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/IlluminationPath_geometry2_small.png b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/IlluminationPath_geometry2_small.png new file mode 100644 index 00000000..87613259 Binary files /dev/null and b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/IlluminationPath_geometry2_small.png differ diff --git a/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq001_orientation.png b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq001_orientation.png new file mode 100644 index 00000000..fb207ba3 Binary files /dev/null and b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq001_orientation.png differ diff --git a/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq002_apparentSza_1.png b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq002_apparentSza_1.png new file mode 100644 index 00000000..767c4948 Binary files /dev/null and b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq002_apparentSza_1.png differ diff --git a/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq002_apparentSza_2.png b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq002_apparentSza_2.png new file mode 100644 index 00000000..5ceb7349 Binary files /dev/null and b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq002_apparentSza_2.png differ diff --git a/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq002_apparentSza_3.png b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq002_apparentSza_3.png new file mode 100644 index 00000000..9bdbe8eb Binary files /dev/null and b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq002_apparentSza_3.png differ diff --git a/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq003_searchPath_1.png b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq003_searchPath_1.png new file mode 100644 index 00000000..406d4d7a Binary files /dev/null and b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq003_searchPath_1.png differ diff --git a/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq003_searchPath_2.png b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq003_searchPath_2.png new file mode 100644 index 00000000..9fe3190a Binary files /dev/null and b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq003_searchPath_2.png differ diff --git a/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq003_searchPath_3.png b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq003_searchPath_3.png new file mode 100644 index 00000000..8f00a633 Binary files /dev/null and b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq003_searchPath_3.png differ diff --git a/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq003_searchPath_4.png b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq003_searchPath_4.png new file mode 100644 index 00000000..71af0317 Binary files /dev/null and b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq003_searchPath_4.png differ diff --git a/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq003_searchPath_5.png b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq003_searchPath_5.png new file mode 100644 index 00000000..705943ad Binary files /dev/null and b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq003_searchPath_5.png differ diff --git a/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq003_searchPath_6.png b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq003_searchPath_6.png new file mode 100644 index 00000000..067918bc Binary files /dev/null and b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq003_searchPath_6.png differ diff --git a/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq003_searchPath_7.png b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq003_searchPath_7.png new file mode 100644 index 00000000..608a052d Binary files /dev/null and b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/aatsr/images/eq003_searchPath_7.png differ diff --git a/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/help.hs b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/help.hs new file mode 100644 index 00000000..6b261657 --- /dev/null +++ b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/help.hs @@ -0,0 +1,24 @@ + + + + + IdePix AATSR Help + + top + + + + TOC + + javax.help.TOCView + toc.xml + + + Search + + javax.help.SearchView + JavaHelpSearch + + diff --git a/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/images/snap_header.jpg b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/images/snap_header.jpg new file mode 100644 index 00000000..d034251e Binary files /dev/null and b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/images/snap_header.jpg differ diff --git a/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/map.jhm b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/map.jhm new file mode 100644 index 00000000..05fa8c97 --- /dev/null +++ b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/map.jhm @@ -0,0 +1,12 @@ + + + + + + + + + + diff --git a/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/style.css b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/style.css new file mode 100644 index 00000000..6f023170 --- /dev/null +++ b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/style.css @@ -0,0 +1,181 @@ +body { + background-color: #FFFFFF; + color: #000000; + font-family: Verdana, Arial, Helvetica, sans-serif; + font-size: 12; + padding: 6pt; + margin: 6pt; +} + +p { + font-family: Verdana, Arial, Helvetica, sans-serif; + margin: 6pt; +} + +.i1 { + margin-left: 30; +} + +.i2 { + margin-left: 50; +} + +.note { + font-family: Verdana, Arial, Helvetica, sans-serif; + color: red; +} + +.inote { + font-family: Verdana, Arial, Helvetica, sans-serif; + background-color: #DDDDDD; +} + +.code { + color: blue; + font-family: courier new, courier, Monospaced, cursive; + text-align: left; + font-weight: normal; + text-decoration: none; + font-style: normal; +} + +h1 { + font-family: Verdana, Arial, Helvetica, sans-serif; + font-size: 25; + font-weight: bold; + margin-top: 18pt; + margin-right: 6pt; + margin-bottom: 12pt; + margin-left: 6pt; +} + +h2 { + font-family: Verdana, Arial, Helvetica, sans-serif; + font-size: 21; + font-weight: bold; + margin-top: 18pt; + margin-right: 6pt; + margin-bottom: 12pt; + margin-left: 6pt; +} + +h3 { + font-family: Verdana, Arial, Helvetica, sans-serif; + font-size: 18; + font-weight: bold; + margin-top: 18pt; + margin-right: 6pt; + margin-bottom: 12pt; + margin-left: 6pt; +} + +h4 { + font-family: Verdana, Arial, Helvetica, sans-serif; + font-size: 15; + font-weight: bold; + margin-top: 18pt; + margin-right: 6pt; + margin-bottom: 6pt; + margin-left: 6pt; +} + +h5 { + font-family: Verdana, Arial, Helvetica, sans-serif; + font-size: 12; + font-weight: bold; + margin-top: 18pt; + margin-right: 6pt; + margin-bottom: 6pt; + margin-left: 6pt; +} + +h6 { + font-family: Verdana, Arial, Helvetica, sans-serif; + font-size: 11; + font-weight: bold; + margin-top: 18pt; + margin-right: 6pt; + margin-bottom: 6pt; + margin-left: 6pt; +} + +li { + margin-top: 4pt; + margin-bottom: 4pt; +} + +table { + border-width: 0; + border-collapse: collapse; + padding: 1px; + margin: 0; +} + +invisibletable { + border-width: 0; + border-style: solid; + border-color: #FFFFFF; + border-collapse: collapse; + padding: 1px; + margin: 0; + font-family: Verdana, Arial, Helvetica, sans-serif; + font-size: 12px; + font-style: normal; +} + +th { + font-size: 12; + text-align: left; + border-width: 1; + border-color: #BBBBBB; + border-style: solid; + border-collapse: collapse; + padding-left: 6px; + padding-right: 6px; + padding-top: 3px; + padding-bottom: 3px; +} + +td { + font-size: 11; + border-width: 1; + border-color: #BBBBBB; + border-style: solid; + border-collapse: collapse; + border-spacing: 0; + padding: 0; + padding-left: 6px; + padding-right: 6px; + vertical-align: top; +} + +table.header { + width: 100%; + height: 29; + font-size: 17; + color: #ffffff; + vertical-align: middle; + background-color: #274351; + border-width: 0; + padding: 0; + margin: 0; +} + +.header { + font-size: 17; + height: 29; + color: #ffffff; + vertical-align: middle; + background-color: #274351; + border-width: 0; + padding: 0; + white-space: nowrap; +} + +hr { + height: 4; + color: #000063; + border-color: #000063; + border-style: solid; + border-width: 1; +} diff --git a/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/toc.xml b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/toc.xml new file mode 100644 index 00000000..fd6596d2 --- /dev/null +++ b/idepix-aatsr/src/main/javahelp/org/esa/snap/idepix/aatsr/docs/toc.xml @@ -0,0 +1,22 @@ + + + + + + + + + + + + + + + + + + diff --git a/idepix-aatsr/src/main/nbm/manifest.mf b/idepix-aatsr/src/main/nbm/manifest.mf new file mode 100644 index 00000000..1d767356 --- /dev/null +++ b/idepix-aatsr/src/main/nbm/manifest.mf @@ -0,0 +1,14 @@ +Manifest-Version: 1.0 +AutoUpdate-Show-In-Client: true +AutoUpdate-Essential-Module: false +OpenIDE-Module-Java-Dependencies: Java > 1.8 +OpenIDE-Module-Display-Category: SNAP Supported Plugins +OpenIDE-Module-Long-Description:

Classification of pixels (cloud, snow, ice, land, water) originating from AATSR 4th reprocessing data

+

Vendor: Brockmann Consult GmbH

+

Release notes: Release notes on GitHub +

Contact address: Chrysanderstr. 1, 21029 Hamburg (Germany)

+

Copyright: (C) 2022 by Brockmann Consult GmbH

+

Vendor: Brockmann Consult GmbH

+

License: GPLv3

+OpenIDE-Module-Layer: layer.xml + diff --git a/idepix-aatsr/src/main/resources/META-INF/services/org.esa.snap.core.gpf.OperatorSpi b/idepix-aatsr/src/main/resources/META-INF/services/org.esa.snap.core.gpf.OperatorSpi new file mode 100644 index 00000000..aa028b2c --- /dev/null +++ b/idepix-aatsr/src/main/resources/META-INF/services/org.esa.snap.core.gpf.OperatorSpi @@ -0,0 +1 @@ +org.esa.snap.idepix.aatsr.IdepixAatsrOp$Spi diff --git a/idepix-aatsr/src/main/resources/helpset.xml b/idepix-aatsr/src/main/resources/helpset.xml new file mode 100644 index 00000000..b5e1b35c --- /dev/null +++ b/idepix-aatsr/src/main/resources/helpset.xml @@ -0,0 +1,7 @@ + + + + + + diff --git a/idepix-aatsr/src/main/resources/layer.xml b/idepix-aatsr/src/main/resources/layer.xml new file mode 100644 index 00000000..938eb6dc --- /dev/null +++ b/idepix-aatsr/src/main/resources/layer.xml @@ -0,0 +1,44 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/idepix-aatsr/src/test/java/org/esa/snap/idepix/aatsr/IdepixAatsrOpTest.java b/idepix-aatsr/src/test/java/org/esa/snap/idepix/aatsr/IdepixAatsrOpTest.java new file mode 100644 index 00000000..7d095b06 --- /dev/null +++ b/idepix-aatsr/src/test/java/org/esa/snap/idepix/aatsr/IdepixAatsrOpTest.java @@ -0,0 +1,207 @@ +/* + * Copyright (C) 2022 Brockmann Consult GmbH (info@brockmann-consult.de) + * + * This program is free software; you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation; either version 3 of the License, or (at your option) + * any later version. + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for + * more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, see http://www.gnu.org/licenses/ + */ + +package org.esa.snap.idepix.aatsr; + +import org.esa.snap.core.datamodel.Mask; +import org.esa.snap.core.datamodel.Product; +import org.esa.snap.core.datamodel.ProductData; +import org.esa.snap.core.gpf.GPF; +import org.esa.snap.core.gpf.OperatorException; +import org.esa.snap.core.gpf.OperatorSpi; +import org.esa.snap.core.gpf.OperatorSpiRegistry; +import org.esa.snap.core.util.DummyProductBuilder; +import org.esa.snap.core.util.math.Range; +import org.junit.Test; + +import java.awt.Color; +import java.awt.Rectangle; +import java.util.List; + +import static org.junit.Assert.*; + +public class IdepixAatsrOpTest { + + @Test + public void testOperatorSpiIsLoaded() { + OperatorSpiRegistry registry = GPF.getDefaultInstance().getOperatorSpiRegistry(); + OperatorSpi operatorSpi = registry.getOperatorSpi("Idepix.Aatsr"); + assertNotNull(operatorSpi); + assertEquals("Idepix.Aatsr", operatorSpi.getOperatorAlias()); + assertNotNull(operatorSpi.getOperatorDescriptor()); + assertSame(operatorSpi.getOperatorClass(), operatorSpi.getOperatorDescriptor().getOperatorClass()); + } + + @Test + public void testTargetProductSignature() { + final Product aatsr = createDummyAatsrSource(); + final IdepixAatsrOp idepixAatsrOp = new IdepixAatsrOp(); + idepixAatsrOp.setSourceProduct(aatsr); + + idepixAatsrOp.setParameterDefaultValues(); + + final Product targetProduct = idepixAatsrOp.getTargetProduct(); + assertEquals(aatsr.getName() + "_idepix", targetProduct.getName()); + assertEquals("AATSR_IDEPIX", targetProduct.getProductType()); + assertEquals(aatsr.getSceneRasterSize(), targetProduct.getSceneRasterSize()); + + } + + @Test + public void arangeTest() { + double[] expected = {-22.5, -21.5, -20.5, -19.5, -18.5, -17.5, -16.5, -15.5, -14.5, -13.5, -12.5, -11.5, -10.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5}; + final double[] arange = IdepixAatsrOp.arange(-22.5, 0.5, 1); + assertArrayEquals(expected, arange, 1.0e-6); + } + + @Test + public void validate_ValidSourceProduct() { + final Product aatsr = createDummyAatsrSource(); + + try { + final IdepixAatsrOp idepixAatsrOp = new IdepixAatsrOp(); + idepixAatsrOp.validate(aatsr); + } catch (Throwable t) { + fail("No exception expected here! Source should be valid."); + } + + } + + @Test(expected = OperatorException.class) + public void validate_InvalidSourceProduct_WrongType() { + final Product aatsr = createDummyAatsrSource(); + aatsr.setProductType("differentType"); + + final IdepixAatsrOp idepixAatsrOp = new IdepixAatsrOp(); + idepixAatsrOp.validate(aatsr); + } + + @Test(expected = OperatorException.class) + public void validate_InvalidSourceProduct_MissingBand() { + final Product aatsr = createDummyAatsrSource(); + aatsr.removeBand(aatsr.getBand("cloud_in")); + + final IdepixAatsrOp idepixAatsrOp = new IdepixAatsrOp(); + idepixAatsrOp.validate(aatsr); + } + + @Test + public void detectFirstLastMaskedPixel() { + final Product aatsr = createDummyAatsrSource(); + // Col 0 = 25-30 + // Col 560 = 10-12 + // Col 1300 = 300-400 + final Mask testMask = Mask.BandMathsType.create("TEST_MASK", "", aatsr.getSceneRasterWidth(), aatsr.getSceneRasterHeight(), + "(X==0.5 && Y>=25.5 && Y<=30.5) ||" + + "(X==560.5 && Y>=10.5 && Y<=12.5) ||" + + "(X==1300.5 && Y>=300.5 && Y<=400.5)", + Color.yellow, 0.5f); + + + aatsr.addMask(testMask); + int[] range; + range = IdepixAatsrOp.detectMaskedPixelRangeInColumn(testMask, 0); + assertEquals(25, range[0], 1.0e-6); + assertEquals(30, range[1], 1.0e-6); + + range = IdepixAatsrOp.detectMaskedPixelRangeInColumn(testMask, 560); + assertEquals(10, range[0], 1.0e-6); + assertEquals(12, range[1], 1.0e-6); + + range = IdepixAatsrOp.detectMaskedPixelRangeInColumn(testMask, 1300); + assertEquals(300, range[0], 1.0e-6); + assertEquals(400, range[1], 1.0e-6); + + } + + @Test + public void sliceRectangle_atZero() { + final List rectangles = + IdepixAatsrOp.sliceRect(new Rectangle(0, 0, 512, 43138), 2000); + + assertEquals(22, rectangles.size()); + assertEquals(new Rectangle(0, 6000, 512, 2000), rectangles.get(3)); + assertEquals(new Rectangle(0, 42000, 512, 1138), rectangles.get(21)); + } + + @Test + public void sliceRectangle_withYOffset() { + final List rectangles = + IdepixAatsrOp.sliceRect(new Rectangle(0, 12356, 512, 20222), 2000); + + assertEquals(11, rectangles.size()); + assertEquals(new Rectangle(0, 12356, 512, 2000), rectangles.get(0)); + assertEquals(new Rectangle(0, 18356, 512, 2000), rectangles.get(3)); + assertEquals(new Rectangle(0, 32356, 512, 222), rectangles.get(10)); + } + + @Test + public void createDistanceArray() { + final int[] distanceArray = IdepixAatsrOp.createDistanceArray(25, 1000); + assertEquals(51, distanceArray.length); + assertEquals(-25000, distanceArray[0]); + assertEquals(-24000, distanceArray[1]); + assertEquals(-23000, distanceArray[2]); + assertEquals(24000, distanceArray[distanceArray.length - 2]); + assertEquals(25000, distanceArray[distanceArray.length - 1]); + } + + @Test + public void calcPathAndTheoreticalHeight() { + // input values taken from Python + float maxObjAlt = 6000; + float minSurfAlt = -65; + double orientation = 154.91488; + float oza = 16.4364f; + float saa = 313.2662f; + float spatialRes = 1000; + float sza = 83.9982f; + float x_tx = 173500.0f; + + // expected values taken from Python + double[] illuPathHeight = {-40.56836, 53.990723, 148.30078, 186.63232, 281.31152, 375.73584, 469.8916, 413.7661, 508.5747, 603.1211, 697.3926, 735.772, 830.4507, 924.84717, 1018.9448, 962.89453, 1057.7168, 1152.2485, 1246.4727, 1284.9097, 1379.5894, 1473.9512, 1512.0195, 1606.8594, 1701.3726, 1795.5366, 1834.0474, 1928.7266, 2023.0447, 2061.1404, 2156.0017, 2250.4902, 2344.5764, 2288.1187, 2383.1836, 2477.8623, 2572.122, 2610.2534, 2705.1443, 2799.6, 2893.5803, 2837.186, 2932.3176, 3026.9954, 3121.175, 3159.3547, 3254.2874, 3348.6963, 3386.2214, 3481.4473, 3576.1243, 3670.1855, 3708.437, 3803.4297, 3897.7683, 3935.1997, 4030.5703, 4125.245, 4219.1143, 4161.652, 4257.481, 4352.5723, 4446.79, 4484.059, 4579.6772, 4674.348, 4767.853, 4710.1294, 4806.434, 4901.715, 4995.68, 5032.597, 5128.74, 5223.395, 5257.624, 5355.065, 5450.8574, 5543.962, 5579.5537, 5677.532, 5771.981, 5796.0537, 5898.027, 6000.0}; + int[][] illuPathSteps = {{-22, -55}, {-22, -54}, {-22, -53}, {-21, -53}, {-21, -52}, {-21, -51}, {-21, -50}, {-20, -51}, {-20, -50}, {-20, -49}, {-20, -48}, {-19, -48}, {-19, -47}, {-19, -46}, {-19, -45}, {-18, -46}, {-18, -45}, {-18, -44}, {-18, -43}, {-17, -43}, {-17, -42}, {-17, -41}, {-16, -41}, {-16, -40}, {-16, -39}, {-16, -38}, {-15, -38}, {-15, -37}, {-15, -36}, {-14, -36}, {-14, -35}, {-14, -34}, {-14, -33}, {-13, -34}, {-13, -33}, {-13, -32}, {-13, -31}, {-12, -31}, {-12, -30}, {-12, -29}, {-12, -28}, {-11, -29}, {-11, -28}, {-11, -27}, {-11, -26}, {-10, -26}, {-10, -25}, {-10, -24}, {-9, -24}, {-9, -23}, {-9, -22}, {-9, -21}, {-8, -21}, {-8, -20}, {-8, -19}, {-7, -19}, {-7, -18}, {-7, -17}, {-7, -16}, {-6, -17}, {-6, -16}, {-6, -15}, {-6, -14}, {-5, -14}, {-5, -13}, {-5, -12}, {-5, -11}, {-4, -12}, {-4, -11}, {-4, -10}, {-4, -9}, {-3, -9}, {-3, -8}, {-3, -7}, {-2, -7}, {-2, -6}, {-2, -5}, {-2, -4}, {-1, -4}, {-1, -3}, {-1, -2}, {0, -2}, {0, -1}, {0, 0}}; + double thresHeight = 101.97246511092132; + + final IdepixAatsrOp.PathAndHeightInfo pathAndHeightInfo = IdepixAatsrOp.calcPathAndTheoreticalHeight(sza, saa, oza, x_tx, orientation, (int) spatialRes, (int) maxObjAlt, minSurfAlt); + + assertArrayEquals(illuPathHeight, pathAndHeightInfo.illuPathHeight, 1.0e-2); + for (int i = 0; i < illuPathSteps.length; i++) { + int[] illuPathStep = illuPathSteps[i]; + assertArrayEquals(illuPathStep, pathAndHeightInfo.illuPathSteps[i]); + } + assertEquals(thresHeight, pathAndHeightInfo.threshHeight, 1.0e-2); + } + + private Product createDummyAatsrSource() { + final DummyProductBuilder pb = new DummyProductBuilder(); + pb.size(DummyProductBuilder.Size.MEDIUM); + pb.gc(DummyProductBuilder.GC.PER_PIXEL); + final Product product = pb.create(); + product.setProductType("ENV_AT_1_RBT"); + product.addBand("solar_zenith_tn", ProductData.TYPE_FLOAT32); + product.addBand("solar_azimuth_tn", ProductData.TYPE_FLOAT32); + product.addBand("sat_zenith_tn", ProductData.TYPE_FLOAT32); + product.addBand("confidence_in", ProductData.TYPE_FLOAT32); + product.addBand("elevation_in", ProductData.TYPE_FLOAT32); + product.addBand("cloud_in", ProductData.TYPE_FLOAT32); + product.addBand("latitude_tx", ProductData.TYPE_FLOAT32); + product.addBand("longitude_tx", ProductData.TYPE_FLOAT32); + product.addBand("x_tx", ProductData.TYPE_FLOAT32); + return product; + } + +} diff --git a/idepix-aatsr/src/test/java/org/esa/snap/idepix/aatsr/OrientationOpImageTest.java b/idepix-aatsr/src/test/java/org/esa/snap/idepix/aatsr/OrientationOpImageTest.java new file mode 100644 index 00000000..7e900704 --- /dev/null +++ b/idepix-aatsr/src/test/java/org/esa/snap/idepix/aatsr/OrientationOpImageTest.java @@ -0,0 +1,67 @@ +/* + * Copyright (c) 2022. Brockmann Consult GmbH (info@brockmann-consult.de) + * + * This program is free software; you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation; either version 3 of the License, or (at your option) + * any later version. + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for + * more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, see http://www.gnu.org/licenses/ + * + */ + +package org.esa.snap.idepix.aatsr; + +import org.esa.snap.core.datamodel.ProductData; +import org.esa.snap.core.util.ImageUtils; +import org.junit.Test; + +import java.awt.image.RenderedImage; + +import static org.junit.Assert.assertEquals; + +/** + * @author Marco Peters + */ +public class OrientationOpImageTest { + + + @Test + public void testOrientationImage() { + float[] lats = { + 56.0214f, 55.7504f, 55.4608f, 56.0214f, + 55.1538f, 55.7504f, 55.4608f, 55.1538f, + 57.3770f, 57.3770f, 57.1410f, 56.8849f, + 57.1410f, 56.6093f, 56.3146f, 56.8849f, + }; + float[] lons = { + -175.4993f, -177.0226f, -178.5242f, -175.4993f, + +179.9959f, -177.0226f, -178.5242f, -179.9959f, + -172.8034f, -174.4020f, -175.9795f, -174.4020f, + -177.5348f, -179.0667f, -175.9795f, +179.4245f + }; + + final RenderedImage latImage = ImageUtils.createRenderedImage(4, 4, ProductData.createInstance(lats)); + final RenderedImage lonImage = ImageUtils.createRenderedImage(4, 4, ProductData.createInstance(lons)); + final OrientationOpImage orientation = new OrientationOpImage(latImage, lonImage); + assertEquals(4, orientation.getWidth()); + assertEquals(4, orientation.getHeight()); + + assertEquals(+162.342737, orientation.getData().getSampleDouble(0,0,0), 1.0e-6); + assertEquals(-179.914133, orientation.getData().getSampleDouble(1,1,0), 1.0e-6); + assertEquals(90.0, orientation.getData().getSampleDouble(2,2,0), 1.0e-6); + assertEquals(-0.165765, orientation.getData().getSampleDouble(3,3,0), 1.0e-6); + } + + @Test + public void computeOrientation() { + assertEquals(-45.004363, OrientationOpImage.computeOrientation(1, 2, 2, 3f), 1.0e-6); + assertEquals(-26.917511, OrientationOpImage.computeOrientation(10f, 10f, 11f, 12f), 1.0e-6); + } + +} \ No newline at end of file diff --git a/idepix-aatsr/src/test/java/org/esa/snap/idepix/aatsr/RunOpMain.java b/idepix-aatsr/src/test/java/org/esa/snap/idepix/aatsr/RunOpMain.java new file mode 100644 index 00000000..eedc5d42 --- /dev/null +++ b/idepix-aatsr/src/test/java/org/esa/snap/idepix/aatsr/RunOpMain.java @@ -0,0 +1,52 @@ +/* + * Copyright (c) 2022. Brockmann Consult GmbH (info@brockmann-consult.de) + * + * This program is free software; you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation; either version 3 of the License, or (at your option) + * any later version. + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for + * more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, see http://www.gnu.org/licenses/ + * + */ + +package org.esa.snap.idepix.aatsr; + +import org.esa.snap.core.dataio.ProductIO; +import org.esa.snap.core.datamodel.Product; +import org.esa.snap.core.gpf.GPF; +import org.esa.snap.core.util.SystemUtils; + +import java.io.IOException; +import java.time.Duration; +import java.time.Instant; +import java.util.HashMap; +import java.util.logging.Level; + +/** + * @author Marco Peters + */ +public class RunOpMain { + + public static void main(String[] args) throws IOException { + SystemUtils.init3rdPartyLibs(RunOpMain.class); +// final Product aatsr = ProductIO.readProduct("H:/SENTINEL3/AATSR4RP/v2.0.5/ENV_AT_1_RBT____20021129T235200_20021130T013735_20210315T024827_6334_011_359______DSI_R_NT_004.SEN3/xfdumanifest.xml"); +// final Product aatsr = ProductIO.readProduct("H:\\related\\QA4EO\\AATSR4th Cloud Shadow\\ENV_AT_1_RBT____20020810T083508_20020810T102042_20210303T040313_6334_008_264______DSI_R_NT_004.dim"); + final Product aatsr = ProductIO.readProduct("H:\\related\\QA4EO\\AATSR4th Cloud Shadow\\ENV_AT_1_RBT____20090615T133401_20090615T151936_20210625T140648_6334_079_496______DSI_R_NT_004.SEN3"); + + Instant start = Instant.now(); + + final HashMap parameters = new HashMap<>(); + parameters.put("copySourceBands", false); + parameters.put("cloudTopHeight", 6000); + final Product shadowProduct = GPF.createProduct("Idepix.Aatsr", parameters, aatsr); + ProductIO.writeProduct(shadowProduct, "H:\\related\\QA4EO\\AATSR4th Cloud Shadow\\" + aatsr.getName() + "_shadowOnly.dim", "BEAM-DIMAP"); + Instant stop = Instant.now(); + SystemUtils.LOG.log(Level.INFO, "DURATION: " + Duration.between(start, stop)); + } +} diff --git a/pom.xml b/pom.xml index 83665b3a..45ccb6d0 100644 --- a/pom.xml +++ b/pom.xml @@ -35,8 +35,10 @@ + idepix-core + idepix-landsat8 idepix-meris