Example - Rasterizing Point Data

[1]:
import json
from functools import partial

import geopandas
from shapely.geometry import box, mapping

from geocube.api.core import make_geocube
from geocube.rasterize import rasterize_points_griddata, rasterize_points_radial

%matplotlib inline

Load in geopackage data and add CRS

[2]:
gdf = geopandas.read_file(
    "../../test/test_data/input/time_vector_data.geojson",
    crs="epsg:4326"
)
[3]:
gdf.head()
[3]:
test_attr test_str_attr test_time_attr geometry
0 1.3 dcf86619 5/21/2016 10:09:21 AM -05:00 POINT (-47.26681 44.21932)
1 1.3 dcf86619 5/21/2016 10:09:21 AM -05:00 POINT (-47.26680 44.21932)
2 1.9 dcf86619 5/21/2016 10:09:21 AM -05:00 POINT (-47.26681 44.21932)
3 1.3 dcf86619 5/21/2016 10:09:21 AM -05:00 POINT (-47.26680 44.21932)
4 1.3 dcf86619 5/21/2016 10:09:21 AM -05:00 POINT (-47.26679 44.21932)

Convert to raster with GeoCube

See docs for make_geocube

Load into grid with griddata nearest resampling

[4]:
geo_grid = make_geocube(
    vector_data=gdf,
    measurements=['test_attr'],
    resolution=(-0.1, 0.00001),
    rasterize_function=rasterize_points_griddata,
)
[5]:
geo_grid
[5]:
<xarray.Dataset>
Dimensions:      (x: 12, y: 11)
Coordinates:
  * y            (y) float64 45.25 45.15 45.05 44.95 ... 44.55 44.45 44.35 44.25
  * x            (x) float64 -47.27 -47.27 -47.27 ... -47.27 -47.27 -47.27
    spatial_ref  int64 0
Data variables:
    test_attr    (y, x) float64 1.3 1.3 1.3 1.2 1.3 1.3 ... 2.3 1.3 1.9 1.3 1.3
Attributes:
    grid_mapping:  spatial_ref
[6]:
# mask nodata and plot
geo_grid.test_attr.where(geo_grid.test_attr!=geo_grid.test_attr.rio.nodata).plot()
[6]:
<matplotlib.collections.QuadMesh at 0x7f91dda7e710>
../_images/examples_rasterize_point_data_9_1.png

Load into grid with griddata cubic resampling

[7]:
geo_grid = make_geocube(
    vector_data=gdf,
    measurements=['test_attr'],
    resolution=(-0.1, 0.00001),
    rasterize_function=partial(rasterize_points_griddata, method="cubic"),
)
[8]:
# mask nodata and plot
geo_grid.test_attr.where(geo_grid.test_attr!=geo_grid.test_attr.rio.nodata).plot()
[8]:
<matplotlib.collections.QuadMesh at 0x7f91d9b1f2b0>
../_images/examples_rasterize_point_data_12_1.png

Load into user-defined grid with radial linear resampling

[9]:
geo_grid = make_geocube(
    vector_data=gdf,
    measurements=['test_attr'],
    geom=json.dumps(mapping(box(-48, 44, -47, 45))),
    output_crs="epsg:3857",
    resolution=(-300, 300),
    rasterize_function=rasterize_points_radial,
)
[10]:
# mask nodata and plot
geo_grid.test_attr.where(geo_grid.test_attr!=geo_grid.test_attr.rio.nodata).plot()
[10]:
<matplotlib.collections.QuadMesh at 0x7f91d9ac8320>
../_images/examples_rasterize_point_data_15_1.png

Load into user-defined grid with radial cubic resampling

Note: This example is simply to demonstrate the option and and would need to be adjusted for better results. As seen in the plot below, the values swelled quite a bit due to the input generating an ill-conditioned matrix.

[11]:
geo_grid = make_geocube(
    vector_data=gdf,
    measurements=['test_attr'],
    geom=json.dumps(mapping(box(-48, 44, -47, 45))),
    output_crs="epsg:3857",
    resolution=(-300, 300),
    rasterize_function=partial(rasterize_points_radial, method="cubic", filter_nan=True),
)
/home/snowal/miniconda/envs/geocube/lib/python3.6/site-packages/scipy/interpolate/rbf.py:268: LinAlgWarning: Ill-conditioned matrix (rcond=3.16678e-18): result may not be accurate.
  self.nodes = linalg.solve(self.A, self.di)
[12]:
# mask nodata and plot
geo_grid.test_attr.where(geo_grid.test_attr!=geo_grid.test_attr.rio.nodata).plot()
[12]:
<matplotlib.collections.QuadMesh at 0x7f91d99fcdd8>
../_images/examples_rasterize_point_data_18_1.png