.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/geo/ex_geohash.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. or to run this example in your browser via Binder .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_geo_ex_geohash.py: ******* Geohash ******* Geohashing is a geocoding method used to encode geographic coordinates (latitude and longitude) into a short string of digits and letters. This string, also known as a geohash, delineates a rectangular area on a map, called a cell. The longer the string, the smaller the cell and the more precise the location. This library provides two ways to work with geohashes: * The :py:class:`pyinterp.GeoHash` class allows for an object-oriented manipulation of a geohash. * The :py:mod:`pyinterp.geohash` module provides functions for vectorized computation of geohashes. In general, you should use the `pyinterp.GeoHash` class when you want to manipulate a few geohashes, and the `pyinterp.geohash` module when you want to process a large number of geohashes stored in `numpy` arrays. The following sections will demonstrate how to use both approaches. Geohash Grid ============ First, let's define a few functions to visualize geohash grids. These functions are for illustrative purposes and are not part of the `pyinterp` library. .. GENERATED FROM PYTHON SOURCE LINES 30-42 .. code-block:: Python import timeit import cartopy.crs import matplotlib.colors import matplotlib.patches import matplotlib.pyplot import numpy import pandas import pyinterp .. GENERATED FROM PYTHON SOURCE LINES 43-44 Writing a visualization routine for GeoHash grids. .. GENERATED FROM PYTHON SOURCE LINES 44-123 .. code-block:: Python def _sort_colors(colors): """Sort colors by hue, saturation, value and name in descending order.""" by_hsv = sorted( (tuple(matplotlib.colors.rgb_to_hsv(matplotlib.colors.to_rgb(color))), name) for name, color in colors.items()) return [name for hsv, name in reversed(by_hsv)] def _plot_box(ax, code, color, caption=True): """Plot a GeoHash bounding box.""" box = pyinterp.GeoHash.from_string(code.decode()).bounding_box() x0 = box.min_corner.lon x1 = box.max_corner.lon y0 = box.min_corner.lat y1 = box.max_corner.lat dx = x1 - x0 dy = y1 - y0 box = matplotlib.patches.Rectangle((x0, y0), dx, dy, alpha=0.5, color=color, ec='black', lw=1, transform=cartopy.crs.PlateCarree()) ax.add_artist(box) if not caption: return rx, ry = box.get_xy() cx = rx + box.get_width() * 0.5 cy = ry + box.get_height() * 0.5 ax.annotate(code.decode(), (cx, cy), color='w', weight='bold', fontsize=16, ha='center', va='center') def plot_geohash_grid(precision, polygon=None, caption=True, color_list=None, inc=7): """Plot geohash bounding boxes.""" color_list = color_list or matplotlib.colors.CSS4_COLORS fig = matplotlib.pyplot.figure(figsize=(24, 12)) ax = fig.add_subplot(1, 1, 1, projection=cartopy.crs.PlateCarree()) if polygon is not None: box = polygon.envelope() if isinstance( polygon, pyinterp.geodetic.Polygon) else polygon ax.set_extent( [ box.min_corner.lon, box.max_corner.lon, box.min_corner.lat, box.max_corner.lat, ], crs=cartopy.crs.PlateCarree(), ) else: box = None colors = _sort_colors(color_list) ic = 0 codes = pyinterp.geohash.bounding_boxes(polygon, precision=precision) color_codes = {codes[0][0]: colors[ic]} for item in codes: prefix = item[precision - 1] if prefix not in color_codes: ic += inc color_codes[prefix] = colors[ic % len(colors)] _plot_box(ax, item, color_codes[prefix], caption) ax.stock_img() ax.coastlines() ax.grid() .. GENERATED FROM PYTHON SOURCE LINES 124-125 Bounds of geohash with a precision of 1 character. .. GENERATED FROM PYTHON SOURCE LINES 125-127 .. code-block:: Python plot_geohash_grid(1) .. image-sg:: /auto_examples/geo/images/sphx_glr_ex_geohash_001.png :alt: ex geohash :srcset: /auto_examples/geo/images/sphx_glr_ex_geohash_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 128-129 Bounds of the geohash ``d`` with a precision of two characters. .. GENERATED FROM PYTHON SOURCE LINES 129-131 .. code-block:: Python plot_geohash_grid(2, polygon=pyinterp.GeoHash.from_string('d').bounding_box()) .. image-sg:: /auto_examples/geo/images/sphx_glr_ex_geohash_002.png :alt: ex geohash :srcset: /auto_examples/geo/images/sphx_glr_ex_geohash_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 132-133 Bounds of the geohash ``dd`` with a precision of three characters. .. GENERATED FROM PYTHON SOURCE LINES 133-135 .. code-block:: Python plot_geohash_grid(3, polygon=pyinterp.GeoHash.from_string('dd').bounding_box()) .. image-sg:: /auto_examples/geo/images/sphx_glr_ex_geohash_003.png :alt: ex geohash :srcset: /auto_examples/geo/images/sphx_glr_ex_geohash_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 136-137 Bounds of the geohash ``dds`` with a precision of four characters. .. GENERATED FROM PYTHON SOURCE LINES 137-140 .. code-block:: Python plot_geohash_grid(4, polygon=pyinterp.GeoHash.from_string('dds').bounding_box()) .. image-sg:: /auto_examples/geo/images/sphx_glr_ex_geohash_004.png :alt: ex geohash :srcset: /auto_examples/geo/images/sphx_glr_ex_geohash_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 141-147 Object-Oriented Interface ========================= The :py:class:`pyinterp.GeoHash` class allows for an object-oriented manipulation of a geohash. You can create a geohash from a longitude and latitude, and then access its properties. .. GENERATED FROM PYTHON SOURCE LINES 147-153 .. code-block:: Python code = pyinterp.GeoHash(-67.5, 22.5, precision=4) print(f'The geohash is {code!s}') print(f'The precision of the geohash is {code.precision()}') print(f'The number of bits is {code.number_of_bits()}') print(f'The center of the cell is at {code.center()}') .. rst-class:: sphx-glr-script-out .. code-block:: none The geohash is ds00 The precision of the geohash is 4 The number of bits is 20 The center of the cell is at (-67.3242, 22.5879) .. GENERATED FROM PYTHON SOURCE LINES 154-155 You can also use this class to get the neighboring geohash codes. .. GENERATED FROM PYTHON SOURCE LINES 155-157 .. code-block:: Python [str(item) for item in code.neighbors()] .. rst-class:: sphx-glr-script-out .. code-block:: none ['ds01', 'ds03', 'ds02', 'debr', 'debp', 'd7zz', 'dkpb', 'dkpc'] .. GENERATED FROM PYTHON SOURCE LINES 158-165 Vectorized Interface ==================== On the other hand, when you want to encode a large volume of data, you should use functions that work on `numpy` arrays. Generation of dummy data .. GENERATED FROM PYTHON SOURCE LINES 165-171 .. code-block:: Python SIZE = 1_000_000 generator = numpy.random.Generator(numpy.random.PCG64(0)) lon = generator.uniform(-180, 180, SIZE) lat = generator.uniform(-80, 80, SIZE) measures = generator.random(SIZE) .. GENERATED FROM PYTHON SOURCE LINES 172-174 Encoding the data is done using the :py:func:`pyinterp.geohash.encode` function. .. GENERATED FROM PYTHON SOURCE LINES 174-177 .. code-block:: Python codes = pyinterp.geohash.encode(lon, lat, precision=4) print(codes) .. rst-class:: sphx-glr-script-out .. code-block:: none [b'mng5' b'fhv0' b'87d7' ... b'9hfj' b'3yvj' b'5ghs'] .. GENERATED FROM PYTHON SOURCE LINES 178-184 As you can see, the resulting codes are encoding as numpy byte arrays. This is to save memory and speed up the processing. This algorithm is very fast, which makes it possible to process a lot of data quickly. The following benchmark measures the time it takes to encode the coordinates into geohashes. .. GENERATED FROM PYTHON SOURCE LINES 184-193 .. code-block:: Python elapsed = timeit.timeit('pyinterp.geohash.encode(lon, lat, precision=12)', number=50, globals={ 'pyinterp': pyinterp, 'lon': lon, 'lat': lat }) / 50 print(f'Elapsed time: {elapsed:.6f} seconds') .. rst-class:: sphx-glr-script-out .. code-block:: none Elapsed time: 0.036425 seconds .. GENERATED FROM PYTHON SOURCE LINES 194-196 The inverse operation is also possible using the :py:func:`pyinterp.geohash.decode` function. .. GENERATED FROM PYTHON SOURCE LINES 196-198 .. code-block:: Python lon, lat = pyinterp.geohash.decode(codes) .. GENERATED FROM PYTHON SOURCE LINES 199-201 You can also use the :py:func:`pyinterp.geohash.transform` to transform geohashes from one precision to another. .. GENERATED FROM PYTHON SOURCE LINES 201-204 .. code-block:: Python codes = pyinterp.geohash.transform(codes, precision=1) print(codes) .. rst-class:: sphx-glr-script-out .. code-block:: none [b'0' b'1' b'2' b'3' b'4' b'5' b'6' b'7' b'8' b'9' b'b' b'c' b'd' b'e' b'f' b'g' b'h' b'j' b'k' b'm' b'n' b'p' b'q' b'r' b's' b't' b'u' b'v' b'w' b'x' b'y' b'z'] .. GENERATED FROM PYTHON SOURCE LINES 205-208 .. code-block:: Python codes = pyinterp.geohash.transform(codes, precision=3) print(codes) .. rst-class:: sphx-glr-script-out .. code-block:: none [b'000' b'001' b'004' ... b'zzv' b'zzy' b'zzz'] .. GENERATED FROM PYTHON SOURCE LINES 209-212 The :py:func:`pyinterp.geohash.bounding_boxes` function allows calculating the geohash codes contained in a box or a polygon. This function allows, for example, to obtain all the geohash codes present on the Mediterranean Sea. .. GENERATED FROM PYTHON SOURCE LINES 212-284 .. code-block:: Python MEDITERRANEAN_SEA = [(-1.43504, 35.38124), (-1.68901, 35.18381), (-1.93947, 35.18664), (-2.18994, 35.18945), (-2.44041, 35.19223), (-2.69089, 35.19498), (-3.19185, 35.20043), (-3.44234, 35.20312), (-4.19382, 35.21103), (-4.44432, 35.21363), (-4.69483, 35.21619), (-4.94221, 35.42047), (-5.18954, 35.62438), (-5.18619, 35.82515), (-5.18272, 36.02538), (-5.17911, 36.22509), (-4.92467, 36.42114), (-4.41929, 36.61311), (-4.16856, 36.60978), (-3.91785, 36.60641), (-3.66714, 36.60302), (-3.41643, 36.59959), (-1.64252, 37.35734), (-0.61682, 38.11193), (-0.05974, 39.61617), (-0.05188, 39.80296), (0.47608, 40.34624), (1.0068, 40.88237), (1.26823, 41.05687), (3.17405, 43.12814), (3.95946, 43.44167), (4.21176, 43.4303), (4.46402, 43.41885), (5.4578, 43.20103), (7.00043, 43.46998), (7.2684, 43.62615), (8.32692, 44.07445), (8.59664, 44.22593), (8.86679, 44.37629), (9.11878, 44.36187), (12.30323, 45.45208), (12.82897, 45.57254), (13.35502, 45.69103), (13.60653, 45.6726), (14.56628, 45.29165), (14.77339, 44.96545), (16.35039, 43.4468), (16.83351, 43.25904), (17.31681, 43.07169), (18.26742, 42.54036), (18.50146, 42.36803), (18.73581, 42.19562), (19.45557, 41.83742), (22.59567, 40.39801), (23.86967, 40.66103), (24.38198, 40.79528), (25.14395, 40.91539), (25.39346, 40.90312), (25.87922, 40.72341), (27.3756, 40.65146), (27.3629, 40.4972), (27.35049, 40.34214), (28.89327, 36.67855), (30.6478, 36.80647), (30.89745, 36.80017), (31.14709, 36.79388), (31.63983, 36.61137), (32.1328, 36.42873), (34.38535, 36.54538), (34.64117, 36.70704), (34.89068, 36.70094), (35.89499, 36.84229), (36.14443, 36.83607), (36.13812, 36.67067), (35.84461, 35.31519), (35.83989, 35.14037), (35.83535, 34.96453), (35.83096, 34.78764), (35.82673, 34.60973), (35.82266, 34.43077), (35.56521, 34.07323), (35.30493, 33.52656), (35.04582, 32.96983), (34.78535, 32.21131), (34.53086, 31.82665), (34.27678, 31.43758), (34.02492, 31.24227), (33.77499, 31.24401), (33.52506, 31.24576), (33.27513, 31.24751), (33.02519, 31.24927), (32.77526, 31.25102), (32.52532, 31.25279), (32.27538, 31.25455), (30.52777, 31.46553), (29.27227, 30.87518), (29.0223, 30.87679), (28.27422, 31.08288), (27.52624, 31.28865), (26.5284, 31.49606), (25.53067, 31.70343), (25.28069, 31.70546), (24.03548, 32.11364), (23.28805, 32.31846), (20.02259, 30.93524), (19.77083, 30.7316), (19.51917, 30.52699), (19.26761, 30.32143), (19.01759, 30.32277), (18.51909, 30.5327), (17.77237, 30.94967), (17.0241, 31.15993), (16.02586, 31.37173), (15.77581, 31.37348), (15.52992, 31.78362), (13.79181, 32.80886), (12.7914, 32.81862), (12.54129, 32.82104), (12.04387, 33.02684), (11.54655, 33.23236), (10.30532, 33.84446), (10.05858, 34.04584), (10.06218, 34.24361), (7.87667, 36.9875), (7.12516, 37.00275), (6.87464, 37.00779), (5.86671, 36.83664), (5.61618, 36.84139), (5.36563, 36.84612), (5.11508, 36.85081), (4.86453, 36.85549), (3.61163, 36.87846), (2.60396, 36.70273), (2.35336, 36.70699), (1.34599, 36.52884), (0.07965, 35.95831), (-0.67601, 35.77038), (-1.43504, 35.38124)] polygon = pyinterp.geodetic.Polygon( [pyinterp.geodetic.Point(lon, lat) for lon, lat in MEDITERRANEAN_SEA]) .. GENERATED FROM PYTHON SOURCE LINES 285-288 .. code-block:: Python precision = 4 plot_geohash_grid(precision, polygon=polygon, caption=False) .. image-sg:: /auto_examples/geo/images/sphx_glr_ex_geohash_005.png :alt: ex geohash :srcset: /auto_examples/geo/images/sphx_glr_ex_geohash_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 289-296 Density calculation =================== Finally, we will calculate the density of points in each geohash cell. For this, we will use the :py:func:`pyinterp.geohash.encode` function to get the geohash of each point, and then we will group the points by geohash and count the number of points in each cell. .. GENERATED FROM PYTHON SOURCE LINES 296-305 .. code-block:: Python df = pandas.DataFrame({ 'lon': lon, 'lat': lat, 'measures': measures, 'geohash': pyinterp.geohash.encode(lon, lat, precision=3) }) df.set_index('geohash', inplace=True) df = df.groupby('geohash').count()['measures'].rename('count').to_frame() .. GENERATED FROM PYTHON SOURCE LINES 306-308 Then, we calculate the density of points in each cell by dividing the number of points by the area of the cell in square kilometers. .. GENERATED FROM PYTHON SOURCE LINES 308-311 .. code-block:: Python df['density'] = df['count'] / ( pyinterp.geohash.area(df.index.values.astype('S')) / 1e6) .. GENERATED FROM PYTHON SOURCE LINES 312-315 Finally, we can visualize the density of points on a map. For this, we will use the :py:func:`pyinterp.geohash.to_xarray` function to convert the geohash and the density into an :py:class:`xarray.DataArray`. .. GENERATED FROM PYTHON SOURCE LINES 315-328 .. code-block:: Python array = pyinterp.geohash.to_xarray(df.index.values.astype('S'), df.density.values) array = array.where(array != 0, numpy.nan) fig = matplotlib.pyplot.figure(figsize=(10, 5)) ax = fig.add_subplot(111, projection=cartopy.crs.PlateCarree()) array.plot(ax=ax, transform=cartopy.crs.PlateCarree(), cmap='viridis', robust=True) ax.coastlines() ax.gridlines(draw_labels=True) ax.set_global() .. image-sg:: /auto_examples/geo/images/sphx_glr_ex_geohash_006.png :alt: ex geohash :srcset: /auto_examples/geo/images/sphx_glr_ex_geohash_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 6.202 seconds) .. _sphx_glr_download_auto_examples_geo_ex_geohash.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/CNES/pangeo-pyinterp/master?urlpath=lab/tree/notebooks/auto_examples/geo/ex_geohash.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: ex_geohash.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: ex_geohash.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: ex_geohash.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_