Download notebook (.ipynb)

Parameter use_crs#

The use_crs parameter is supported by all geometry layers which have the map parameter, i.e. by layers which can be used for visualization of Geopandas GeoDataFrame.

Those layers are: point, path, polygon, text, label and map.

By default, Lets-Plot converts all coordinates in every GeoDataFrame to degrees of longitude and latitude and later projects these map coordinates to screen coordinates using Mercator projection.

Sometimes it might be necessary to use some other coordinate reference system (CRS) instead of the default WGS84 / Mercator.

If this is the case, you can use the use_crs parameter to specify EPSG code of CRS to project GeoDataFrame coordinates to.

import numpy as np

from lets_plot import *

LetsPlot.setup_html()

Let’s say we have a raster image (for example, a satellite or aerial photograph, or a DEM image) that covers a square lend area with:

  • left-bottom corner at coordinates (25_000, 4_500_000)

  • right-top corner at coordinates (525_000, 5_000_000)

In practice, the image coordinates and the image coordinate reference system (CRS) are known from a metadata accompanying the image.

Here, for the sake of simplicity, we will use just a simple generated image 2 x 2 pix.

We will also presume that our image coordinates use UTM zone 33N CRS (aka EPSG:32633).

Thus, the image coordinates are given in meters instead of more familiar longitude and latitude.

# image 2 x 2 pix
image_arr = np.array([
    [[150, 0, 0], [0, 150, 0]],
    [[0, 0, 150], [150, 150, 0]]
])

image_CRS = "EPSG:32633"

# Image coordinates
min_x = 25_000
max_x = 525_000
min_y = 4_500_000
max_y = 5_000_000

ggplot() + geom_imshow(image_arr, extent=[min_x, max_x, min_y, max_y])

Combining image and map#

Step 1 here would be to acquire a Geopandas GeoDataFrame containing map of the area of interest.

In this demo we will use GeoDataFrame containing a map of Italy which we will obtain using Lets-Plot geocoding module.

from lets_plot.geo_data import *

The geodata is provided by © OpenStreetMap contributors and is made available here under the Open Database License (ODbL).

italy = geocode_countries('italy').inc_res().get_boundaries()

ggplot() + geom_map(map=italy)

italy.crs

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

The “Italy map” and our image use different CRS-s#

The image uses UTM CRS with coordinates in meters and the “Italy map” uses the EPSG:4326 CRS with coordinates in degrees of longitude and latitude.

Due to this discrepancy we can not right away combine these two layes in one plot.

The use_crs parameter#

The solution here is actually pretty straightforward: we will just let the map layer know which exactly CRS the other plot layer uses (i.e. which CRS our image uses).

ggplot() + geom_map(map=italy, use_crs=image_CRS) + labs(x="", y="")

# Image and map
( ggplot()  
  + geom_imshow(image_arr, extent=[min_x, max_x, min_y, max_y])
  + geom_map(map=italy, use_crs=image_CRS) 
  + labs(x="", y="")
)

# Select better theme and colors
( ggplot()
  + geom_imshow(image_arr, extent=[min_x, max_x, min_y, max_y])
  + geom_map(map=italy, use_crs=image_CRS, color="#a6b6ba", size=.8)
  + scale_x_continuous(format=".3~s")
  + scale_y_continuous(format=".3~s")
  + labs(x="", y="")
  + theme_bw() + flavor_solarized_dark()
  + ggsize(500, 500)
)