Plotting raster and vector data together in python¶
In this tutorial we will learn how to plot raster and vector data together in python.
We will be using cartopy,rioxarray, xarray and geopandas in our tutorial
This tutorial will focus more on rioxrray and xarray component, if you want to learn more about geopandas read this tutorial.
What is xarray and rioarray?¶
xarray and rioxarray are two powerful libraries in Python for working with multidimensional numerical arrays, specifically in the context of earth observation data.
xarray is a library for working with labeled multi-dimensional arrays, similar to the functionality of pandas for one-dimensional data. It provides a convenient way to store and manipulate large, multi-dimensional arrays of numerical data, along with metadata such as labels for dimensions and coordinates. xarray also integrates with popular libraries for data analysis and visualization, such as pandas and matplotlib.
rioxarray is a library built on top of xarray, specifically for working with geospatial raster data. It provides a convenient way to read, write, and manipulate raster data stored in a variety of file formats, including geospatial raster formats like GeoTIFF and NetCDF. It also includes functionality for handling spatial metadata such as projection information and transforming between different projections. With rioxarray, users can work with raster data in a familiar, pandas-like way, using multi-dimensional arrays and metadata to analyze and visualize geospatial data.
What is Cartopy?¶
Cartopy is a Python library for creating maps and working with geospatial data. It provides a high-level interface for working with geospatial data, including the ability to manipulate and visualize data in different map projections. Cartopy also provides tools for working with shapefiles and other geospatial data formats, as well as geospatial transformations and aggregation.
One of the key features of Cartopy is its ability to work with a variety of map projections, including both cylindrical and conic projections. This allows users to visualize their data in a variety of ways, depending on the specifics of their analysis. Cartopy also integrates well with other popular Python libraries for data analysis, such as numpy and pandas, making it a convenient tool for working with geospatial data in the context of a wider data analysis workflow
import geopandas as gpd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import xarray
import rioxarray
coastlines=gpd.read_file("Coastlines/Global_coastlines_low_res.shp")
topography=rioxarray.open_rasterio("ETopo.tif")
topography
<xarray.DataArray (band: 1, y: 10801, x: 21601)>
[233312401 values with dtype=float32]
Coordinates:
* band (band) int64 1
* x (x) float64 -180.0 -180.0 -180.0 -179.9 ... 180.0 180.0 180.0
* y (y) float64 90.0 89.98 89.97 89.95 ... -89.97 -89.98 -90.0
spatial_ref int64 0
Attributes:
_FillValue: -3.4028234663852886e+38
scale_factor: 1.0
add_offset: 0.0This is a description of a xarray DataArray object representing a single band of a raster dataset. The data is represented by a 2D array with 10801 rows and 21601 columns, for a total of 2,333,124,01 values. The data type of the values is float32.
The coordinates section lists the variables that define the spatial reference of the data. x and y are the longitude and latitude coordinates of each cell in the raster, respectively, with a range from -180 to 180 degrees for longitude and from -90 to 90 degrees for latitude. The band coordinate refers to the specific band within the raster data, in this case band 1.
The spatial_ref coordinate is a single integer value with no additional information provided, likely referring to the coordinate reference system used for the data.
The Attributes section lists additional metadata associated with the data, including the FillValue (representing missing or invalid data), “scale_factor” (a multiplier for the data values), and “add_offset” (an additive value to be applied to the data values).
fig=plt.figure(figsize=[12,8])
ax = fig.add_axes([0,0,1,1],projection=ccrs.PlateCarree())
raster_image=topography.plot(ax=ax)
coastlines.plot(ax=ax,color='none', edgecolor="black",linewidth=0.4)
plt.show()
This code block is an example of how to create a map visualization of a raster dataset using the Cartopy library in Python.
The fig variable is assigned the result of the figure method from the matplotlib library, which creates a new figure for plotting. The figsize argument is set to [12, 8], which specifies the size of the figure in inches.
The ax variable is assigned the result of the add_axes method, which creates an axes object on the figure. The projection argument is set to ccrs.PlateCarree(), which specifies that the map should be in a Plate Carree projection (a simple cylindrical projection that is often used for global maps).
The raster_image variable is assigned the result of the plot method from the xarray library, which creates a raster image plot of the topography dataset. The ax argument is set to the axes object created earlier, which specifies that the plot should be added to the existing axes.
The coastlines.plot method is used to add coastlines to the plot, which are displayed with a black edge color and a linewidth of 0.4. The color argument is set to ‘none’, which specifies that the coastlines should not have a fill color.
Finally, the plt.show() method is called to display the plot on the screen
Changing colormap and colorbar in xarray plot¶
fig=plt.figure(figsize=[12,8])
ax = fig.add_axes([0,0,1,1],projection=ccrs.PlateCarree())
raster_image=topography.plot(ax=ax,
cmap='terrain',
vmax=6000,vmin=-4000,
label="Elevation with Coastlines",
cbar_kwargs={'orientation':'horizontal','shrink':0.5,
'label':'ELEVATION (m)'})
plt.title("Elevation with Coastlines")
coastlines.plot(ax=ax,color='none', edgecolor="black",linewidth=0.4)
plt.show()
The cmap argument is set to ‘terrain’, which specifies that the terrain color map should be used for the plot. The vmax and vmin arguments are set to 6000 and -4000 respectively, which specify the maximum and minimum values for the color map. This allows for control over the range of colors used to represent the data.
The label argument is set to “Elevation with Coastlines”, which provides a label for the plot in the legend. The cbar_kwargs argument is set to a dictionary that specifies properties of the color bar. The orientation key is set to ‘horizontal’, which specifies that the color bar should be displayed horizontally. The shrink key is set to 0.5, which specifies that the color bar should be reduced by half. The label key is set to ‘ELEVATION (m)’, which provides a label for the color bar.