Raster to regular points

This example demonstrates the conversion of a raster regular-grid values to a point cloud using geoutils.Raster.to_points().

We open a raster.

import geoutils as gu

filename_rast = gu.examples.get_path("exploradores_aster_dem")
rast = gu.Raster(filename_rast)
rast = rast.crop([rast.bounds.left, rast.bounds.bottom, rast.bounds.left + 500, rast.bounds.bottom + 500])

Let’s plot the raster.

rast.plot(cmap="terrain")

We convert the raster to points. By default, this returns a vector with column geometry burned.

pc = rast.to_pointcloud()
pc

PointCloud(
  ds=              b1                geometry
       0    2787.605957  POINT (627175 4834055)
       1    2791.827393  POINT (627205 4834055)
       2    2796.588623  POINT (627235 4834055)
       3    2802.831299  POINT (627265 4834055)
       4    2809.177734  POINT (627295 4834055)
       ..           ...                     ...
       284  2765.868408  POINT (627535 4833575)
       285  2768.098389  POINT (627565 4833575)
       286  2770.859375  POINT (627595 4833575)
       287  2773.972168  POINT (627625 4833575)
       288  2777.362305  POINT (627655 4833575)

       [289 rows x 2 columns]
  crs=PROJCS["WGS 84 / UTM zone 18S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-75],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32718"]]
  bounds=BoundingBox(left=np.float64(627175.0), bottom=np.float64(4833575.0), right=np.float64(627655.0), top=np.float64(4834055.0)))

We plot the point vector.

pc.plot(ax="new", cmap="terrain", legend=True)