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.0
This 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()