1

I am trying to export raster values using multipolygon shapefile in python. I have found the answer here, but the calculation there is not valid for multipolygon. Could please someone guide me, how i should correct the code in order to have not polygon but multipolygon datatype in calculation.

My code is below:

import rasterio
from rasterio.mask import mask
import geopandas as gpd
import numpy as np
from rasterio import Affine
from shapely.geometry import mapping

shapefile = gpd.read_file(r'/Users..../polygon_sector.shp')
geoms = shapefile.geometry.values
geometry = geoms[0] # shapely geometry


transform to GeJSON format


geoms = [mapping(geoms[0])]


extract the raster values within the polygon


with rasterio.open("/Users/.../map_reclass.tif") as src:
     out_image, out_transform = mask(src, geoms, crop=True)


no data values of the original raster


no_data=src.nodata
print(no_data)


extract the values of the masked array


data = out_image[0,:,:]


extract the row, columns of the valid values


row, col = np.where(data != no_data)
rou = np.extract(data != no_data, data)


affine import Affine


T1 = out_transform * Affine.translation(0.5, 0.5)  # reference the pixel centre
rc2xy = lambda r, c: (c, r) * T1


d = gpd.GeoDataFrame({'col':col,'row':row,'ROU':rou})


coordinate transformation


d['x'] = d.apply(lambda row: rc2xy(row.row,row.col)[0], axis=1)
d['y'] = d.apply(lambda row: rc2xy(row.row,row.col)[1], axis=1)


geometry


from shapely.geometry import Point
d['geometry'] =d.apply(lambda row: Point(row['x'], row['y']), axis=1)


save to a shapefile


d.to_file(r'/Users/y.../result_full.shp', driver='ESRI Shapefile') ```
Lusia
  • 101
  • 7

1 Answers1

3

I found the solution, basically going backwards and deleted the part, where i start defining precisely the polygon. As result the code looks the following:

import rasterio
from rasterio.mask import mask
import geopandas as gpd
import numpy as np
from rasterio import Affine
from shapely.geometry import mapping

shapefile = gpd.read_file(r'/Users..../polygon_sector.shp')
geoms = shapefile.geometry.values


with rasterio.open("/Users/.../map_reclass.tif") as src:
     out_image, out_transform = mask(src, geoms, crop=True)


no data values of the original raster


no_data=src.nodata
print(no_data)


extract the values of the masked array


data = out_image[0,:,:]


extract the row, columns of the valid values


row, col = np.where(data != no_data)
rou = np.extract(data != no_data, data)


affine import Affine


T1 = out_transform * Affine.translation(0.5, 0.5)  # reference the pixel centre
rc2xy = lambda r, c: (c, r) * T1


d = gpd.GeoDataFrame({'col':col,'row':row,'ROU':rou})


coordinate transformation


d['x'] = d.apply(lambda row: rc2xy(row.row,row.col)[0], axis=1)
d['y'] = d.apply(lambda row: rc2xy(row.row,row.col)[1], axis=1)


geometry


from shapely.geometry import Point
d['geometry'] =d.apply(lambda row: Point(row['x'], row['y']), axis=1)


save to a shapefile


d.to_file(r'/Users/y.../result_full.shp', driver='ESRI Shapefile') ```
Lusia
  • 101
  • 7