1 #-*- coding: utf-8 -*-
2 importrasterio3 importjson4 from rasterio.mask importmask5
6 defClip(Tiff, Geodata):7 rasterfile =Tiff8 geoms = json.loads(Geodata) #解析string格式的geojson数据
9 rasterdata =rasterio.open(rasterfile)10 member = 0.0 #记录总人数
11 Gridnumber = 0 #记录相交区域总像素点数
12 Transform = rasterdata._transform #得到影像六参数
13
14 for i in range(len(geoms['features'])):15 geo = [geoms['features'][i]['geometry']]16 #掩模得到相交区域
17 out_image, out_transform = mask(rasterdata, geo, all_touched=True, crop=True, nodata=rasterdata.nodata)18 out_list =out_image.tolist()19 out_list =out_list[0]20
21 for k inrange(len(out_list)):22 for j inrange(len(out_list[k])):23 if out_list[k][j] >=0:24 member +=out_list[k][j]25 Gridnumber += 1
26
27 #人数单位为万人,小数点后保留两位
28 print(round(member / 2500, 2))29 #面积单位为平方公里,小数点后保留两位
30 print(round(Gridnumber * Transform[0] * Transform[3] / 250000, 2))31
32 if __name__ == '__main__':33 geojson = "{\"type\":\"FeatureCollection\",\"features\":[{\"type\":\"Feature\",\"geometry\":{\"type\":\"Polygon\",\"coordinates\":[[[114.4,21.2],[118.25777710514137,21.2],[118.25890798273747,21.299120898984143],[118.25711853337589,21.398166452889882],[118.25240412231551,21.49706137436084],[118.24476226308434,21.59573049143711],[118.23419263830453,21.694098805135695],[118.22069711944573,21.792091546888855],[118.20427978543194,21.88963423579395],[118.1849469400189,21.98665273562841],[118.16270712785808,22.08307331158244],[118.13757114915893,22.178822686664432],[118.10955207286088,22.27382809773303],[118.07866524822043,22.368017351110097],[118.04492831472203,22.461318877731173],[118.00836121021473,22.553661787787973],[117.96898617717864,22.644975924820585],[117.926827767026,22.735191919215197],[117.88191284233778,22.82424124106683],[117.83427057694053,22.912056252364096],[117.78393245372968,22.998570258456652],[117.73093226014316,23.083717558764647],[117.67530608