I will start with #2 and #3. There is a tool out there called geocube: https://corteva.github.io/geocube/stable/index.html. It will handle all of the rasterization for you in a simpler process and handles multiple columns of data.
https://corteva.github.io/geocube/stable/examples/grid_to_vector_map.html
import geopandas as gpd
import numpy as np
from shapely.geometry import Polygon
from geocube.api.core import make_geocube
xmin,ymin,xmax,ymax = 1,1,5,5
lenght = 1
wide = 1
cols = list(range(int(np.floor(xmin)), int(np.ceil(xmax)), wide))
rows = list(range(int(np.floor(ymin)), int(np.ceil(ymax)), lenght))
rows.reverse()
polygons = []
for x in cols:
for y in rows:
polygons.append( Polygon([(x,y), (x+wide, y), (x+wide, y-lenght), (x, y-lenght)]) )
g = gpd.GeoDataFrame(
{"data1": list(range(len(polygons))), "data2": list(range(10, 10+len(polygons)))},
geometry=polygons,
crs={"init": "epsg:4326"}
)
cube = make_geocube(vector_data=g, resolution=(1, -1))
This will provide the data in a gridded xarray interface with the proper geospatial metadata set.
<xarray.Dataset>
Dimensions: (x: 4, y: 4)
Coordinates:
* y (y) float64 0.5 1.5 2.5 3.5
* x (x) float64 4.5 3.5 2.5 1.5
spatial_ref int64 0
Data variables:
data1 (y, x) float64 15.0 11.0 7.0 3.0 14.0 ... 1.0 12.0 8.0 4.0 0.0
data2 (y, x) float64 25.0 21.0 17.0 13.0 24.0 ... 22.0 18.0 14.0 10.0
The data looks like:
You can then export each column of the data to a tif with:
cube.data1.rio.to_raster("data1.tif")
1) Where in the code I am doing the error(s)
The main issue is that you are not adding the transform
for the raster that defines the resolution of the grid cells and its location.
4) why does the g.geometry column contain 5 coordinate pairs, first and last being the same, and not 4?)
This is because shapely is closing the loop of the polygon automatically for you.
Looks like you're only storing the geometry field in the polygons_buffered GDF. You can apply the buffer function directly to the original GDF and either modify the geometry field in place or store in another gdf (which requires you to store more than just the computed field in the new gdf, which is where you had issues with your original code). Since you're exporting immediately as a SHP, I'd go with the former option:
polygons['geometry'] = polygons['geometry'].buffer(polygons['Distance'])
polygons.head()
polygons.plot()
Best Answer
Using a mask and assignment you can achieve this:
Full MWE: