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') ```

1 Answer 1

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') ```

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.