Line rasterization to align the other raster

[guopeng-jiang]

I try to rasterize a road network spatial line (GeoDataFrame) to align my DEM raster (xarray.DataArray) so I can do raster calculations.

from geocube.api.core import make_geocube
from shapely.geometry import mapping,box

bbox = masked_dem.rio.bounds()

roads_data = make_geocube(
    vector_data=roads,
    measurements=['lane_count'],
    output_crs="epsg:2193",
    resolution=masked_dem.rio.resolution(),
    geom=mapping(box(*bbox)),
    fill=0
)

But I get this overflow error below. I don't know why it happens and how to fix it.

---------------------------------------------------------------------------
OverflowError                             Traceback (most recent call last)
Cell In[289], line 9
      5 bbox = masked_dem.rio.bounds()
      6 # bbox = gpd.GeoSeries(box(*masked_dem.rio.bounds()), crs=2193)
      7 # geom = bbox.__geo_interface__["features"][0]["geometry"]
----> 9 roads_data = make_geocube(
     10     vector_data=roads,
     11     measurements=['lane_count'],
     12     output_crs="epsg:2193",
     13     resolution=masked_dem.rio.resolution(),
     14     geom=mapping(box(*bbox)),
     15     # geom=json.dumps(geom),
     16     fill=0
     17 )

File /opt/miniconda/lib/python3.11/site-packages/geocube/api/core.py:99, in make_geocube(vector_data, measurements, datetime_measurements, output_crs, resolution, align, geom, like, fill, group_by, interpolate_na_method, categorical_enums, rasterize_function)
     35 """
     36 Rasterize vector data into an ``xarray`` object.  Each attribute will be a data
     37 variable in the :class:`xarray.Dataset`.
   (...)
     95 
     96 """
     97 geobox_maker = GeoBoxMaker(output_crs, resolution, align, geom, like)
---> 99 return VectorToCube(
    100     vector_data=vector_data,
    101     geobox_maker=geobox_maker,
    102     fill=fill,
    103     categorical_enums=categorical_enums,
    104 ).make_geocube(
    105     measurements=measurements,
    106     datetime_measurements=datetime_measurements,
    107     group_by=group_by,
    108     interpolate_na_method=interpolate_na_method,
    109     rasterize_function=rasterize_function,
    110 )

File /opt/miniconda/lib/python3.11/site-packages/geocube/vector_to_cube.py:79, in VectorToCube.__init__(self, vector_data, geobox_maker, fill, categorical_enums)
     59 """
     60 Initialize the GeoCube class.
     61 
   (...)
     76 
     77 """
     78 self._vector_data = load_vector_data(vector_data)
---> 79 self._geobox = geobox_maker.from_vector(self._vector_data)
     80 self._grid_coords = affine_to_coords(
     81     self._geobox.affine, self._geobox.width, self._geobox.height
     82 )
     83 self._fill = fill if fill is not None else numpy.nan

File /opt/miniconda/lib/python3.11/site-packages/geocube/geo_utils/geobox.py:190, in GeoBoxMaker.from_vector(self, vector_data)
    188     geom_crs = CRS(crs_input)
    189     geopoly = Geometry(self.geom, crs=geom_crs)
--> 190 return GeoBox.from_geopolygon(geopoly, self.resolution, crs, self.align)

File /opt/miniconda/lib/python3.11/site-packages/odc/geo/geobox.py:576, in GeoBox.from_geopolygon(geopolygon, resolution, crs, align, shape, tight, anchor, tol)
    573 else:
    574     geopolygon = geopolygon.to_crs(crs)
--> 576 return GeoBox.from_bbox(
    577     geopolygon.boundingbox,
    578     crs,
    579     shape=shape,
    580     resolution=resolution,
    581     anchor=anchor,
    582     tol=tol,
    583     tight=tight,
    584 )

File /opt/miniconda/lib/python3.11/site-packages/odc/geo/geobox.py:500, in GeoBox.from_bbox(bbox, crs, tight, shape, resolution, anchor, tol)
    498     offy, ny = snap_grid(bbox.bottom, bbox.top, ry, None, tol=tol)
    499 else:
--> 500     offx, nx = snap_grid(bbox.left, bbox.right, rx, _snap.x, tol=tol)
    501     offy, ny = snap_grid(bbox.bottom, bbox.top, ry, _snap.y, tol=tol)
    503 affine = Affine.translation(offx, offy) * Affine.scale(rx, ry)

File /opt/miniconda/lib/python3.11/site-packages/odc/geo/math.py:191, in snap_grid(x0, x1, res, off_pix, tol)
    188     return x1, max(nx, 1)
    190 off = off_pix * abs(res)
--> 191 _tx, nx = _snap_edge(x0 - off, x1 - off, res, tol)
    192 return _tx + off, nx

File /opt/miniconda/lib/python3.11/site-packages/odc/geo/math.py:160, in _snap_edge(x0, x1, res, tol)
    158 if res > 0:
    159     return _snap_edge_pos(x0, x1, res, tol)
--> 160 _tx, nx = _snap_edge_pos(x0, x1, -res, tol)
    161 tx = _tx + nx * (-res)
    162 return tx, nx

File /opt/miniconda/lib/python3.11/site-packages/odc/geo/math.py:150, in _snap_edge_pos(x0, x1, res, tol)
    148 assert res > 0
    149 assert x1 >= x0
--> 150 _x0 = floor(maybe_int(x0 / res, tol))
    151 _x1 = ceil(maybe_int(x1 / res, tol))
    152 nx = max(1, _x1 - _x0)

OverflowError: cannot convert float infinity to integer

Datasets:

roads

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 14882 entries, 0 to 14881
Data columns (total 8 columns):
 #   Column      Non-Null Count  Dtype   
---  ------      --------------  -----   
 0   name        0 non-null      float64 
 1   hway_num    50 non-null     object  
 2   name_id     0 non-null      float64 
 3   lane_count  14882 non-null  int64   
 4   way_count   0 non-null      float64 
 5   status      0 non-null      float64 
 6   surface     14882 non-null  object  
 7   geometry    14882 non-null  geometry
dtypes: float64(4), geometry(1), int64(1), object(2)
memory usage: 930.3+ KB

masked_dem

<xarray.DataArray (band: 1, y: 2563, x: 2089)>
array([[[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]]], dtype=float32)
Coordinates:
  * band         (band) int64 1
  * x            (x) float64 1.09e+06 1.091e+06 ... 1.716e+06 1.717e+06
  * y            (y) float64 5.517e+06 5.517e+06 ... 4.749e+06 4.748e+06
    spatial_ref  int64 0
Attributes:
    AREA_OR_POINT:       Area
    STATISTICS_MAXIMUM:  3576
    STATISTICS_MEAN:     -9999
    STATISTICS_MINIMUM:  0
    STATISTICS_STDDEV:   -9999
    scale_factor:        1.0
    add_offset:          0.0
    long_name:           nztm
    _FillValue:          nan

[snowman2]

Does this work:

roads_data = make_geocube(
    vector_data=roads,
    measurements=['lane_count'],
    like=masked_dem,
    fill=0
)

[guopeng-jiang]

Yes it works, thank you :-)