Rasterize a vector

This example demonstrates the rasterizing of a vector using rasterize().

We open a raster and vector.

import geoutils as gu

filename_rast = gu.examples.get_path("everest_landsat_b4")
filename_vect = gu.examples.get_path("everest_rgi_outlines")
rast = gu.open_raster(filename_rast)
vect = gu.open_vector(filename_vect)

Let’s plot the raster and vector.

rast.rst.plot(cmap="Purples")
vect.vct.plot(ref=rast, fc="none", ec="k", lw=2)

First option: using the raster as a reference to match, we rasterize the vector in any projection and georeferenced grid. We simply have to pass the raster as single argument to rasterize(). See Match-reference operations for more details.

vect_rasterized = vect.vct.rasterize(rast)
vect_rasterized.rst.plot(ax="new", cmap="viridis")

By default, rasterize() will burn the index of the vector’s features in their geometry. We can specify the in_value to burn a single value, or any iterable with the same length as there are features in the vector. An out_value can be passed to burn outside the geometries.

vect_rasterized = vect.vct.rasterize(rast, in_value=1)
vect_rasterized.rst.plot(ax="new")

Note

If the rasterized in_value is fixed to 1 and out_value to 0 (default), then rasterize() is creating a boolean mask. This is equivalent to using create_mask(), and will return a boolean xarray.DataArray.

vect_rasterized

<xarray.DataArray (y: 655, x: 800)> Size: 4MB
array([[1., 1., 1., ..., 1., 1., 1.],
       [1., 1., 1., ..., 1., 1., 1.],
       [1., 1., 1., ..., 1., 1., 1.],
       ...,
       [0., 0., 0., ..., 0., 1., 1.],
       [0., 0., 0., ..., 0., 1., 1.],
       [0., 0., 0., ..., 1., 1., 1.]], shape=(655, 800))
Coordinates:
  * y            (y) float64 5kB 3.108e+06 3.108e+06 ... 3.089e+06 3.089e+06
  * x            (x) float64 6kB 4.78e+05 4.78e+05 ... 5.02e+05 5.02e+05
    spatial_ref  int64 8B 0

xarray.DataArray

• y: 655
• x: 800

• 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 ... 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0

  array([[1., 1., 1., ..., 1., 1., 1.],
         [1., 1., 1., ..., 1., 1., 1.],
         [1., 1., 1., ..., 1., 1., 1.],
         ...,
         [0., 0., 0., ..., 0., 1., 1.],
         [0., 0., 0., ..., 0., 1., 1.],
         [0., 0., 0., ..., 1., 1., 1.]], shape=(655, 800))

• Coordinates: (3)
  • y
    (y)
    float64
    3.108e+06 3.108e+06 ... 3.089e+06

    array([3108125., 3108095., 3108065., ..., 3088565., 3088535., 3088505.],
          shape=(655,))

  • x
    (x)
    float64
    4.78e+05 4.78e+05 ... 5.02e+05

    array([478015., 478045., 478075., ..., 501925., 501955., 501985.], shape=(800,))

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["WGS 84 / UTM zone 45N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",87],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32645"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    projected_crs_name :

    WGS 84 / UTM zone 45N

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    87.0

    false_easting :

    500000.0

    false_northing :

    0.0

    scale_factor_at_central_meridian :

    0.9996

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 45N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",87],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32645"]]

    GeoTransform :

    478000.0 30.0 0.0 3108140.0 0.0 -30.0

    array(0)

Second option: we can pass any georeferencing parameter to rasterize(). Any unpassed attribute will be deduced from the vector itself, except from the raster shape that will default to 1000 x 1000.

# vect_rasterized = vect.vct.rasterize(res=500)
# vect_rasterized.rst.plot()

Important

The raster shape or resolution are the only unknown arguments to rasterize a vector, one or the other can be passed.