.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "example_gallery/raster_operations/parial_overlap.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_example_gallery_raster_operations_parial_overlap.py: .. _example partial overlap: Partial overlap =============== Merge rasters with different spatial extent Introduction ------------ Since each BoundedGrid in GridKit has information on the bounds, multiple grids can be combined in a way where their spatial extend is taken into account. Basic mathematical operations were shown in the example :ref:`ndvi.py `, but in that example the bands were fully overlapping. In this example we are combining grids that are partially overlapping. To combine the grids we will be taking the mean value for the cells that are shared between the grids. In order to keep the scope of this test limited, three slices of tha same geotiff are used. In most real cases, it is likely that these would be separate datasets that don't perfectly match. Note that in these cases, all datasets need to be resampled onto the same grid. Also, since all slices in this example are from the same dataset, the values in the overlapping sections will be identical, making for seamless fusing of the grids. If the values in the overlapping area differ a lot between the datasets, the seams won't be as smooth. That said, let's read in the data. The source of the DEM data used in this example is: https://www.eea.europa.eu/data-and-maps/data/copernicus-land-monitoring-service-eu-dem .. GENERATED FROM PYTHON SOURCE LINES 32-42 .. code-block:: Python from gridkit import read_raster path_dem = "../../tests/data/alps_dem.tiff" dem1 = read_raster(path_dem, bounds=(29040, 167500, 29100, 167620)) dem2 = read_raster(path_dem, bounds=(29070, 167470, 29150, 167580)) dem3 = read_raster(path_dem, bounds=(29080, 167520, 29160, 167600)) .. GENERATED FROM PYTHON SOURCE LINES 44-46 Let's group them for convenience .. GENERATED FROM PYTHON SOURCE LINES 47-50 .. code-block:: Python dem_slices = [dem1, dem2, dem3] .. GENERATED FROM PYTHON SOURCE LINES 51-53 And then plot these data slices individually to see what we are dealing with .. GENERATED FROM PYTHON SOURCE LINES 54-64 .. code-block:: Python import matplotlib.pyplot as plt fig, axes = plt.subplots(1, 3) for ax, dem in zip(axes, dem_slices): ax.imshow(dem, cmap="terrain", extent=dem.mpl_extent) fig.tight_layout() plt.show() .. image-sg:: /example_gallery/raster_operations/images/sphx_glr_parial_overlap_001.png :alt: parial overlap :srcset: /example_gallery/raster_operations/images/sphx_glr_parial_overlap_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 65-68 It is hard to see from these three separate images which parts of the slices overlap. If we visualize the ``count`` per cell, we can get a better idea on the distribution of the data. .. GENERATED FROM PYTHON SOURCE LINES 69-78 .. code-block:: Python import gridkit dem_count = gridkit.count(dem_slices) plt.imshow(dem_count, extent=dem_count.mpl_extent) plt.colorbar() plt.show() .. image-sg:: /example_gallery/raster_operations/images/sphx_glr_parial_overlap_002.png :alt: parial overlap :srcset: /example_gallery/raster_operations/images/sphx_glr_parial_overlap_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 79-83 Now we have an idea on what the spatial distribution of the data looks like, let's take the mean of the datasets to combine them into one. For each cell, this will take the mean value of all grids that cover it. .. GENERATED FROM PYTHON SOURCE LINES 84-91 .. code-block:: Python dem_mean = gridkit.mean(dem_slices) plt.imshow(dem_mean, cmap="terrain", extent=dem_mean.mpl_extent) plt.colorbar() plt.show() .. image-sg:: /example_gallery/raster_operations/images/sphx_glr_parial_overlap_003.png :alt: parial overlap :srcset: /example_gallery/raster_operations/images/sphx_glr_parial_overlap_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 92-97 Beautiful, The data is merged! In this case of course all values of the slices are the same, so the seams are perfectly smooth. We can see the seems if we offset the values of the different slices a bit. .. GENERATED FROM PYTHON SOURCE LINES 98-106 .. code-block:: Python dem_mean = gridkit.mean([dem1 - 100, dem2, dem3 + 100]) plt.imshow(dem_mean, cmap="terrain", extent=dem_mean.mpl_extent) plt.colorbar() plt.show() .. image-sg:: /example_gallery/raster_operations/images/sphx_glr_parial_overlap_004.png :alt: parial overlap :srcset: /example_gallery/raster_operations/images/sphx_glr_parial_overlap_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 107-119 Now the seams become visible. If this is a problem for your dataset, consider cleaning the data before merging or smoothing the data after merging. For smoothing the data afterwards, we can use scipy's gaussian_filter. Note that scipy treats our grid object as a numpy array, so we need to update the existing grid with the array returned from the scipy function. A larger convolution window (sigma) will result in more aggressive filtering. This will result in smoother seams between slices, but also remove high-frequency signal from your data. When merging tiles therefore, a tradeoff will need to be reached that will differ per usecase. .. GENERATED FROM PYTHON SOURCE LINES 120-131 .. code-block:: Python from scipy.ndimage import gaussian_filter filtered_data = gaussian_filter(dem_mean, sigma=1) dem_filtered = dem_mean.update(filtered_data) plt.imshow(dem_filtered, cmap="terrain", extent=dem_mean.mpl_extent) plt.colorbar() plt.show() .. image-sg:: /example_gallery/raster_operations/images/sphx_glr_parial_overlap_005.png :alt: parial overlap :srcset: /example_gallery/raster_operations/images/sphx_glr_parial_overlap_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 132-150 .. Tip :: You may have noticed that area's containing NaNs get exaggerated after filtering. Scipy's implementation of ``gaussian_filter`` does not take NaNs into account. This means each pixel where the kernal contains a NaN will become NaN too. If this is a problem for your dataset, have a look at `astropy's implementation `_. .. Nodata values ------------- For each grid, the nodata value of that grid is also taken into account. This means that values outside of the original extent of each respective data grid are not considered in the calculation. If the data grids have different nodata_values, this still works, for each grid's nodata_value is considered separately. In that scenario the resulting product will have the nodata_value of the first grid. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.668 seconds) .. _sphx_glr_download_example_gallery_raster_operations_parial_overlap.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: parial_overlap.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: parial_overlap.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: parial_overlap.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_