Vector Source#
The wradlib.io.VectorSource class is designed to conveniently handle Vector Data (eg. shapefiles). It originates from the wradlib.zonalstats module but moved to wradlib.io.gdal for better visibility.
managing georeferenced data (grid points or grid polygons, zonal polygons),
output to vector and raster files available within ogr/gdal
geopandas dataframe connector
import wradlib as wrl
import matplotlib.pyplot as plt
import matplotlib as mpl
import warnings
import numpy as np
warnings.filterwarnings("ignore")
The wradlib.io.VectorSource class handles point or polygon vector data by wrapping ogr.DataSource with special functions.
The following example shows how to create different VectorSource objects:
from osgeo import osr
# create gk2 projection osr object
proj_gk2 = osr.SpatialReference()
proj_gk2.ImportFromEPSG(31466)
# Setting up DataSource
box0 = np.array(
[
[2600000.0, 5630000.0],
[2600000.0, 5640000.0],
[2610000.0, 5640000.0],
[2610000.0, 5630000.0],
[2600000.0, 5630000.0],
]
)
box1 = np.array(
[
[2610000.0, 5630000.0],
[2610000.0, 5640000.0],
[2620000.0, 5640000.0],
[2620000.0, 5630000.0],
[2610000.0, 5630000.0],
]
)
box2 = np.array(
[
[2600000.0, 5640000.0],
[2600000.0, 5650000.0],
[2610000.0, 5650000.0],
[2610000.0, 5640000.0],
[2600000.0, 5640000.0],
]
)
box3 = np.array(
[
[2610000.0, 5640000.0],
[2610000.0, 5650000.0],
[2620000.0, 5650000.0],
[2620000.0, 5640000.0],
[2610000.0, 5640000.0],
]
)
point0 = np.array(wrl.georef.get_centroid(box0))
point1 = np.array(wrl.georef.get_centroid(box1))
point2 = np.array(wrl.georef.get_centroid(box2))
point3 = np.array(wrl.georef.get_centroid(box3))
# creates Polygons in Datasource
poly = wrl.io.VectorSource(
np.array([box0, box1, box2, box3]), trg_crs=proj_gk2, name="poly"
)
# creates Points in Datasource
point = wrl.io.VectorSource(
np.vstack((point0, point1, point2, point3)), trg_crs=proj_gk2, name="point"
)
Warning 1: Layer poly has no spatial index, DROP SPATIAL INDEX failed.
Warning 1: Layer point has no spatial index, DROP SPATIAL INDEX failed.
print(poly)
<wradlib.VectorSource>
Type: Polygon
Geometries: 4
Let’s have a look at the data, which will be exported as numpy arrays. The property data exports all available data as numpy arrays:
numpy access#
print(poly.data)
print(point.data)
[[[2600000.0 5630000.0]
[2600000.0 5640000.0]
[2610000.0 5640000.0]
[2610000.0 5630000.0]
[2600000.0 5630000.0]]
[[2610000.0 5630000.0]
[2610000.0 5640000.0]
[2620000.0 5640000.0]
[2620000.0 5630000.0]
[2610000.0 5630000.0]]
[[2600000.0 5640000.0]
[2600000.0 5650000.0]
[2610000.0 5650000.0]
[2610000.0 5640000.0]
[2600000.0 5640000.0]]
[[2610000.0 5640000.0]
[2610000.0 5650000.0]
[2620000.0 5650000.0]
[2620000.0 5640000.0]
[2610000.0 5640000.0]]]
[[2605000.0 5635000.0]
[2615000.0 5635000.0]
[2605000.0 5645000.0]
[2615000.0 5645000.0]]
geopandas access#
poly.geo.explore()
Make this Notebook Trusted to load map: File -> Trust Notebook
point.geo.loc[slice(0, 2)]
| index | geometry | |
|---|---|---|
| 0 | 0 | POINT (2605000 5635000) |
| 1 | 1 | POINT (2615000 5635000) |
| 2 | 2 | POINT (2605000 5645000) |
point.geo.loc[[0, 1, 3]]
| index | geometry | |
|---|---|---|
| 0 | 0 | POINT (2605000 5635000) |
| 1 | 1 | POINT (2615000 5635000) |
| 3 | 3 | POINT (2615000 5645000) |
point.geo.query("index in (0, 2)")
| index | geometry | |
|---|---|---|
| 0 | 0 | POINT (2605000 5635000) |
| 2 | 2 | POINT (2605000 5645000) |
fig = plt.figure()
ax = fig.add_subplot(111)
poly.geo.plot(column="index", ax=ax)
point.geo.plot(ax=ax)
<Axes: >
Now, with the DataSource being created, we can add/set attribute data of the features:
# add attribute
poly.set_attribute("mean", np.array([10.1, 20.2, 30.3, 40.4]))
point.set_attribute("mean", np.array([10.1, 20.2, 30.3, 40.4]))
Attributes associated with features can also be retrieved:
# get attributes
print(poly.get_attributes(["mean"]))
# get attributes filtered
print(poly.get_attributes(["mean"], filt=("index", 2)))
[[10.1, 20.2, 30.3, 40.4]]
[[30.3]]
Currently data can also be retrieved by:
index -
wradlib.io.gdal.VectorSource.get_data_by_idx(),attribute -
wradlib.io.gdal.VectorSource.get_data_by_att()andgeometry -
wradlib.io.gdal.VectorSource.get_data_by_geom().
Using the property mode the output type can be set permanently.
get_data_by_idx#
point.get_data_by_idx([0, 2])
array([[2605000.0, 5635000.0],
[2605000.0, 5645000.0]], dtype=object)
point.get_data_by_idx([0, 2], mode="geo")
| index | mean | geometry | |
|---|---|---|---|
| 0 | 0 | 10.1 | POINT (2605000 5635000) |
| 2 | 2 | 30.3 | POINT (2605000 5645000) |
get_data_by_att#
point.get_data_by_att("index", [0, 2])
array([[2605000.0, 5635000.0],
[2605000.0, 5645000.0]], dtype=object)
point.get_data_by_att("index", [0, 2], mode="geo")
| index | mean | geometry | |
|---|---|---|---|
| 0 | 0 | 10.1 | POINT (2605000 5635000) |
| 2 | 2 | 30.3 | POINT (2605000 5645000) |
get_data_by_geom#
# get OGR.Geometry
geom0 = poly.get_data_by_idx([0], mode="ogr")[0]
# get geopandas Geometry
geom1 = poly.get_data_by_idx([0], mode="geo")
point.get_data_by_geom(geom=geom0)
array([[2605000.0, 5635000.0]], dtype=object)
point.get_data_by_geom(geom=geom0, mode="ogr")
array([<osgeo.ogr.Geometry; proxy of None >], dtype=object)
point.get_data_by_geom(geom=geom1, mode="geo")
| index | mean | geometry | |
|---|---|---|---|
| 0 | 0 | 10.1 | POINT (2605000 5635000) |
Finally, we can export the contained data to OGR/GDAL supported vector and raster files:
# dump as 'ESRI Shapefile', default
poly.dump_vector("test_poly.shp")
point.dump_vector("test_point.shp")
# dump as 'GeoJSON'
poly.dump_vector("test_poly.geojson", driver="GeoJSON")
point.dump_vector("test_point.geojson", driver="GeoJSON")
# dump as 'GTiff', default
poly.dump_raster("test_poly_raster.tif", attr="mean", pixel_size=100.0)
# dump as 'netCDF'
poly.dump_raster("test_poly_raster.nc", driver="netCDF", attr="mean", pixel_size=100.0)
0...10...20...30...40...50...60...70...80...90...100 - done.
0
reload geojson#
point2 = wrl.io.VectorSource("test_point.geojson")
poly2 = wrl.io.VectorSource("test_poly.geojson")
fig = plt.figure()
ax = fig.add_subplot(111)
poly2.geo.plot(column="index", ax=ax)
point2.geo.plot(ax=ax)
<Axes: >
reload raster geotiff#
import xarray as xr
ds = xr.open_dataset("test_poly_raster.tif")
ds.band_data[0].plot()
<matplotlib.collections.QuadMesh at 0x7adac930a900>
reload raster netcdf#
ds = xr.open_dataset("test_poly_raster.nc")
ds.Band1.plot()
<matplotlib.collections.QuadMesh at 0x7adb1fae82d0>
