.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "example_gallery/vector_synergy/intersect_geometries.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_example_gallery_vector_synergy_intersect_geometries.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_example_gallery_vector_synergy_intersect_geometries.py:


Geometry intersections
======================

Obtain the cells that intersect with vector geometries.

DEM data source: https://www.eea.europa.eu/data-and-maps/data/copernicus-land-monitoring-service-eu-dem

Introduction
------------

Vector geometries can be used to obtain the ID of the intersecting grid cells.
These vector geometries can be either in the form of Shapely (Multi) Geometries or GeoPandas GeoSeries.
A list of geometries can be supplied. This list can contain a mix of geometry types.

.. Tip ::
    While it is technically possible to stack Multi-Geometries (for example a MultiPolygon of MultiPolygons), this is not good for performance of ``intersect_geometries``.
    Try to avoid this if possible.

In this example, GeoPandas is used instead of Shapely, for it also allows us to conveniently read in the vector data.

.. Note ::
    GeoPandas is not installed with GridKit by default, while Shapely is.

Reading and preparing the data
------------------------------

After reading the streams and lakes vector files, they are concatenated into a single GeoDataFrame.
This allows us to obtain the intersecting cells in one function call.
Of course it is also an option to obtain the intersecting cell IDs for the strams and lakes separately and combine them after using `numpy.vstack`.
With that approach it might be worthwhile to call a `numpy.unique` on axis 0 afterwards to remove duplicate IDs.

.. GENERATED FROM PYTHON SOURCE LINES 34-48

.. code-block:: Python


    import geopandas
    import matplotlib.pyplot as plt
    import pandas

    from gridkit import read_raster

    rivers = geopandas.read_file("../../tests/data/streams.gpkg")
    lakes = geopandas.read_file("../../tests/data/lakes.gpkg")
    dem = read_raster("../../tests/data/alps_dem.tiff")

    water_bodies = pandas.concat([rivers, lakes]).reset_index().to_crs(dem.crs)

    river_cell_ids = dem.intersect_geometries(water_bodies.geometry)







.. GENERATED FROM PYTHON SOURCE LINES 49-55

Visualization
-------------
For a clear plot, we only need to plot the DEM near the vecor geometries as a backround.
For that reason, the dem is cropped to the extent of the water bodies with a little buffer.
Since ``intersect_geometries`` does it's own cropping, cropping the grid before ``intersect_geometries`` does not increase performance.
A plot showing DEM, the geometries (orange) and the selected cells (red) is then generated as follows:

.. GENERATED FROM PYTHON SOURCE LINES 56-84

.. code-block:: Python


    dem = dem.crop(water_bodies.total_bounds, bounds_crs=water_bodies.crs, buffer_cells=20)


    def generate_plot():
        fig, ax = plt.subplots()
        ax.imshow(dem, extent=dem.mpl_extent)

        for poly in dem.to_shapely(river_cell_ids):
            x, y = poly.exterior.xy
            ax.fill(x, y, alpha=0.5, color="red")

        for geom in water_bodies.geometry:
            try:
                x, y = geom.exterior.xy
                ax.fill(x, y, alpha=0.5, color="orange")
            except AttributeError:
                x, y = geom.xy
                ax.plot(x, y, c="orange")

        ax.set_xlabel("lon")
        ax.set_ylabel("lat")
        return ax


    generate_plot()
    plt.show()




.. image-sg:: /example_gallery/vector_synergy/images/sphx_glr_intersect_geometries_001.png
   :alt: intersect geometries
   :srcset: /example_gallery/vector_synergy/images/sphx_glr_intersect_geometries_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 85-86

A zoom-in of a section of the plot allows us to see more clearly what cells are selected.

.. GENERATED FROM PYTHON SOURCE LINES 86-91

.. code-block:: Python


    ax = generate_plot()
    ax.set_xlim(30040, 30160)
    ax.set_ylim(167230, 167325)
    plt.show()



.. image-sg:: /example_gallery/vector_synergy/images/sphx_glr_intersect_geometries_002.png
   :alt: intersect geometries
   :srcset: /example_gallery/vector_synergy/images/sphx_glr_intersect_geometries_002.png
   :class: sphx-glr-single-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 4.464 seconds)


.. _sphx_glr_download_example_gallery_vector_synergy_intersect_geometries.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: intersect_geometries.ipynb <intersect_geometries.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: intersect_geometries.py <intersect_geometries.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: intersect_geometries.zip <intersect_geometries.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_