<jupyter_start><jupyter_text>Fit DSM to low resolution initial DEMThis notebook details how to estimate and apply the transform to fit A DSM to the low resolution initial DEM. This method is currently implemented in cars.  Notebook parametersThose parameters have to be completed to use the notebook.<jupyter_code># Path to the cars folder
cars_home = "TODO"
# Path to the directory containing the content.json file of the prepare step output
content_dir = "TODO"<jupyter_output><empty_output><jupyter_text>Imports<jupyter_code>### Trick to override cars verision
import sys
sys.path = [cars_home] + sys.path
import os
os.environ['OTB_APPLICATION_PATH'] = os.path.join(cars_home,'build','lib','otb','applications')+':'+os.environ['OTB_APPLICATION_PATH']
###
# Silent OTB info logs
os.environ['OTB_LOGGER_LEVEL']='WARNING'
import xarray as xr
import scipy as sp
import numpy as np
import math
import matplotlib.pyplot as plt
from cars.core import projection
from cars.steps import triangulation
from cars.steps import rasterization
from cars.conf import output_prepare
from cars.core import constants as cst
from scipy.signal import butter, lfilter, filtfilt, lfilter_zi<jupyter_output><empty_output><jupyter_text>Reading input data<jupyter_code>lowres_dsm_from_matches = xr.open_dataset(os.path.join(content_dir,'lowres_dsm_from_matches.nc'))
lowres_initial_dem = xr.open_dataset(os.path.join(content_dir,'lowres_initial_dem.nc'))
conf = output_prepare.read_preprocessing_content_file(os.path.join(content_dir,'content.json'))
img1 = conf['input']['img1']
img2 = conf['input']['img2']
dem = conf['input']['srtm_dir']
matches = np.load(os.path.join(content_dir,'matches.npy'))
disp_to_alt_ratio = conf['preprocessing']['output']['disp_to_alt_ratio']
srtm_dir = conf['input']['srtm_dir']<jupyter_output><empty_output><jupyter_text>Look for direction of oscillations We compute the direction in which oscillations might happen (increasing acquisition time direction)<jupyter_code>vec1 = projection.get_time_ground_direction(img1,dem=srtm_dir)
vec2 = projection.get_time_ground_direction(img2, dem=srtm_dir)
time_direction_vector = (vec1+vec2)/2
print("Direction img1 {} degree wrt horizontal axis<".format(180*math.atan2(vec1[1], vec1[0])/math.pi))
print("Direction img2 {} degree wrt horizontal axis".format(180*math.atan2(vec2[1], vec2[0])/math.pi))
print("Oscillation direction: {} ({} degree wrt horizontal axis)".format(time_direction_vector,180*math.atan2(time_direction_vector[1], time_direction_vector[0])/math.pi))<jupyter_output><empty_output><jupyter_text>Measuring initial difference with low resolution DEM<jupyter_code>dsm_diff = lowres_initial_dem[cst.RASTER_HGT]-lowres_dsm_from_matches[cst.RASTER_HGT]
(dsm_diff).plot(robust=True)<jupyter_output><empty_output><jupyter_text>Extracting difference signal (offset and oscillations) We project differences on the  `(x,y), time_direction_vector` axis.<jupyter_code>origin = [float(dsm_diff[cst.X][0].values), float(dsm_diff[cst.Y][0].values)]
x_values_2d, y_values_2d = np.meshgrid(dsm_diff[cst.X], dsm_diff[cst.Y])
curv_coords = projection.project_coordinates_on_line(x_values_2d, y_values_2d, origin, time_direction_vector)
curv_array = xr.DataArray(dsm_diff.values.ravel(), coords={"curv" : curv_coords.ravel()}, dims = ("curv"))
curv_array = curv_array.dropna(dim='curv')
curv_array = curv_array.sortby('curv')
curv_array.plot(figsize=(15,5))<jupyter_output><empty_output><jupyter_text>The we perform denoising by aggregating median along time axis.<jupyter_code>min_curv = np.min(curv_array.curv)
max_curv = np.max(curv_array.curv)
nbins = int(math.ceil((max_curv-min_curv)/(lowres_dsm_from_matches.attrs[cst.RESOLUTION])))
filtered_curv_array = curv_array.groupby_bins('curv',nbins).median()
filtered_curv_array = filtered_curv_array.rename({'curv_bins': 'curv'})
filtered_curv_array = filtered_curv_array.assign_coords({'curv' : np.array([d.mid for d in filtered_curv_array.curv.data])})<jupyter_output><empty_output><jupyter_text>We also compute the number of point in each median slot, and discard measurements for slots with insufficient number of points (< 100)<jupyter_code>filtered_curv_npoints = curv_array.groupby_bins('curv',nbins).count()
filtered_curv_npoints = filtered_curv_npoints.rename({'curv_bins': 'curv'})
filtered_curv_npoints = filtered_curv_npoints.assign_coords({'curv' : np.array([d.mid for d in filtered_curv_npoints.curv.data])})
filtered_curv_array = filtered_curv_array.where(filtered_curv_npoints > 100).dropna(dim='curv')<jupyter_output><empty_output><jupyter_text>Next, we will apply butterworth low pass filtering to extract low frequency of the difference signal.<jupyter_code>b, a = butter(3, 0.05)
zi = lfilter_zi(b, a)
z, _ = lfilter(b, a, filtered_curv_array.values, zi=zi*filtered_curv_array.values[0])
z2, _ = lfilter(b, a, z, zi=zi*z[0])
filtered_curv_array_lowpass = xr.DataArray(filtfilt(b, a,filtered_curv_array.values),coords=filtered_curv_array.coords)<jupyter_output><empty_output><jupyter_text>Display both curves (median difference and low passed filtered low pass median difference)<jupyter_code>xr.concat([filtered_curv_array,filtered_curv_array_lowpass],dim='tmp').isel(tmp=[0,1]).plot(hue='tmp',figsize=(15,5))<jupyter_output><empty_output><jupyter_text>Smoothed difference modellingWe will use cubic splines to model this smoothed difference signal. Our aim is to find the smoothest cubic splines that have a RMSE < 0.3 meters.<jupyter_code># Initialize s parameter
best_s = 100*len(filtered_curv_array_lowpass.curv)

# Compute first spline and RMSE
splines = sp.interpolate.UnivariateSpline(filtered_curv_array_lowpass.curv,filtered_curv_array_lowpass.values, ext=3, k=3, s=best_s)
estimated_correction = xr.DataArray(splines(filtered_curv_array_lowpass.curv),coords=filtered_curv_array_lowpass.coords)
rmse = (filtered_curv_array_lowpass-estimated_correction).std(dim='curv')

# Loop to find best s
s = [best_s]
rmses = [rmse]
target_rmse = 0.3
while rmse > target_rmse and best_s > 0.001:
    best_s/=2.
    splines = sp.interpolate.UnivariateSpline(filtered_curv_array_lowpass.curv,filtered_curv_array_lowpass.values, ext=3, k=3, s=best_s)
    estimated_correction = xr.DataArray(splines(filtered_curv_array_lowpass.curv),coords=filtered_curv_array_lowpass.coords)
    rmse = (filtered_curv_array_lowpass-estimated_correction).std(dim='curv')
    s.append(best_s)
    rmses.append(rmse)

# Display
fig, ax = plt.subplots()
ax.set_xlabel('value for s')
ax.set_ylabel('rmse')
ax.plot(rmses)
print("Best smoothing factor: {} (rmse = {} meters)".format(best_s, rmse.values))<jupyter_output><empty_output><jupyter_text>Now we can interpolate the spline to obtain the 1D correction<jupyter_code>splines = sp.interpolate.UnivariateSpline(filtered_curv_array_lowpass.curv,filtered_curv_array_lowpass.values, ext=3, k=3, s = best_s)
estimated_correction = xr.DataArray(splines(filtered_curv_array.curv),coords=filtered_curv_array.coords)
xr.concat([filtered_curv_array, filtered_curv_array_lowpass, estimated_correction, filtered_curv_array-estimated_correction],dim='tmp').isel(tmp=[0,1,2,3]).plot(hue='tmp',figsize=(15,5))<jupyter_output><empty_output><jupyter_text>Likewise, we can interpolate the splines in 2D to get the 2D z correction field<jupyter_code>estimated_correction_2d = xr.DataArray(splines(curv_coords),coords=dsm_diff.coords)
xr.concat((dsm_diff,estimated_correction_2d, lowres_initial_dem[cst.RASTER_HGT] - lowres_dsm_from_matches[cst.RASTER_HGT] - estimated_correction_2d), dim='tmp').plot(col='tmp', robust=True, figsize=(15,4))<jupyter_output><empty_output><jupyter_text>The mean column is also visualized (before and after) as well as the correction (in the center).  Applying correction in triangulation In this section we will demonstrate how to use the correction during triangulation, using matches triangulation. first we compute the initial points cloud from matches<jupyter_code>initial_points_cloud = triangulation.triangulate_matches(conf,matches)<jupyter_output><empty_output><jupyter_text>We can now estimate z correction for each triangulated point using our splines<jupyter_code>initial_points_cloud_z_correction = splines(projection.project_coordinates_on_line(initial_points_cloud[cst.X], initial_points_cloud[cst.Y], origin, time_direction_vector))<jupyter_output><empty_output><jupyter_text>We can then convert this z correction into a disparity correction:<jupyter_code>disp_correction = initial_points_cloud_z_correction/disp_to_alt_ratio<jupyter_output><empty_output><jupyter_text>Apply this disparity correction to matches:<jupyter_code>corrected_matches = np.copy(matches)
corrected_matches[:,2]= corrected_matches[:,2]-disp_correction[:,0]<jupyter_output><empty_output><jupyter_text>And triangulate again:<jupyter_code>corrected_points_cloud = triangulation.triangulate_matches(conf,corrected_matches)<jupyter_output><empty_output><jupyter_text>Now we can rasterize corrected cloud.<jupyter_code>startx = float(np.min(lowres_dsm_from_matches[cst.X]).values)-0.5*lowres_dsm_from_matches.resolution
starty = float(np.max(lowres_dsm_from_matches[cst.Y]).values)+0.5*lowres_dsm_from_matches.resolution
sizex = lowres_dsm_from_matches[cst.X].shape[0]
sizey = lowres_dsm_from_matches[cst.Y].shape[0]
corrected_lowres_dsm_from_matches = rasterization.simple_rasterization_dataset([corrected_points_cloud],lowres_dsm_from_matches.resolution,xstart=startx, ystart=starty, xsize=sizex, ysize=sizey, epsg=4326)<jupyter_output><empty_output><jupyter_text>Last, we display diff with low resolution intial dem:<jupyter_code>xr.concat((lowres_initial_dem[cst.RASTER_HGT]-lowres_dsm_from_matches[cst.RASTER_HGT],
           lowres_initial_dem[cst.RASTER_HGT]-corrected_lowres_dsm_from_matches[cst.RASTER_HGT]), dim='tmp').plot(col='tmp', robust=True,figsize=(10,4))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11310975609756098, 'fraction_numerical': 0.011382113821138212, 'mean_word_length': 9.72004357298475, 'pattern_counts': {'":': 0, '<': 84, '<?xml version=': 0, '>': 84, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '305223', 'license': 'Apache-2.0', 'n_tokens_mistral': 3653, 'n_tokens_neox': 3498, 'n_words': 591, 'path': 'notebooks/lowres_dem_fit.ipynb', 'repo_name': 'fchapoton/cars', 'size': '9840'}