工作要解決的問題：
專案根據衛星影像和演算法得到同一區域兩個時間的變化情況，可以得到乙個包含變化圖斑的.shp檔案，但是該區域中不同區域的性質不同，需要進行區別處理。
依據是甲方提供的乙個區域分割的圖斑.shp檔案。
最後，需要輸出的是甲方每個圖斑中真實的變化情況(即前shp檔案和後shp檔案求交集)，還要將兩個shp中的圖斑屬性進行綜合。
最先參考文章
ogr layer intersection
之後(最終使用方案)
intersect shapefiles using shapely
之後還參考
fiona使用者手冊
首先碰到的是兩個shp無交集，但是軟體做有，最後發現是兩個shp的座標參考係不同
在進行圖斑相交處理過程中會碰到
topologyexception: input geom 1
is invalid: ring self-intersection at or near
嘗試此種方法無果。
多邊形自相交處理-selfintersection
最後解決方案參考文章：
修復不合法的polygon
編碼問題
最後綜合的屬性中，中文列印為亂碼，因此將encoding設定為utf-8即可解決
fiona使用過程中
with fiona.
open
('oriented-ccw.shp'
,'w'
, crs=source.crs,
# 該引數報錯
driver=source.driver,
schema=sink_schema,
)as sink:
# 值後嘗試此種
sink = fiona.
open
('/tmp/foo.shp'
,'w'
,**source.meta)
最後傳參使用 crs_wkt 進行傳遞成功
有序字典的合併
merge 兩個 ordereddict 到乙個 ordereddict 中的可用和不可用方法
import time
import fiona
import rtree
defintersect_shp_shp
(manager_shp_path,input_shp_path,output_shp_path,min_area=
none):
""" manager_shp_path : 第一shp檔案路徑 該專案中是 林地小班
input_shp_path : 第二shp檔案路徑 該專案中是 演算法產生的變化圖斑
output_shp_path : 輸出shp檔案路徑
min_area : 用於圖斑面積篩選(面積小於該值的圖斑，將會pass，也提高演算法執行時間)
return none
"""start_=time.time(
)with fiona.
open
(manager_shp_path,
'r',encoding=
'utf-8'
)as layer1:
with fiona.
open
(input_shp_path,
'r',encoding=
'utf-8'
)as layer2:
if layer1.crs[
'init'
]!=layer2.crs[
'init']:
raise exception(
'shapefile的座標系不同，無法進行計算！請先轉換座標系！'
)# 合併兩shp的schema
schema3 = layer2.schema.copy(
) schema3[
'properties'
].update(layer1.schema[
'properties'])
# 新增面積屬性
schema3[
'properties'][
'area']=
'float'
# schema3['geometry']='polygon'
with fiona.
open
(output_shp_path, mode=
'w', schema=schema3,driver=
'esri shapefile'
,crs_wkt=layer2.meta[
'crs_wkt'
],encoding=
'utf-8'
)as layer3:
# 建立rtree索引
print
('開始建立索引......'
) index = rtree.index.index(
) manager_num=
0for feat1 in layer1:
fid =
int(feat1[
'id'])
geom1 = shape(feat1[
'geometry'])
index.insert(fid, geom1.bounds)
manager_num+=
1print
('林地小班數量共有：%d 個，建立索引共耗時%.2f s'
%(manager_num,time.time(
)-start_)
)# 執行相交運算
intersect_execute(layer1,layer2,layer3,index,min_area)
print
('完成本次運算共耗時%.2f s'
%(time.time(
)-start_,))
defintersect_execute
(layer1,layer2,layer3,index,min_area)
:print
('開始進行交集運算......'
) result_num,polygon_num,small_num=0,
0,0for feat2 in layer2:
polygon_num+=
1 geom2 = shape(feat2[
'geometry'])
if min_area and geom2.areasmall_num+=
1continue
# 檢測合法性，並進行合法化(主要針對 有洞的polygon)
ifnot geom2.is_valid:
geom2=geom2.
buffer(0
)for fid in
list
(index.intersection(geom2.bounds)):
feat1 = layer1[fid]
geom1 = shape(feat1[
'geometry'])
ifnot geom1.is_valid:
geom1=geom1.
buffer(0
)if geom1.intersects(geom2)
:# 合併屬性
props = feat2[
'properties'
].copy(
) props.update(feat1[
'properties'])
intersect=geom1.intersection(geom2)
if min_area and intersect.areacontinue
# 給area屬性賦值
props[
'area'
]=intersect.area
try: layer3.write(
) result_num+=
1except exception as e:
print
(e)print
('輸入圖斑數量 %d 個，小面積的圖斑有 %d 個，得到符合要求的圖斑共： %d 個'
%(polygon_num,small_num,result_num,))
if __name__ ==
"__main__"
: intersect_shp_shp(
)
shp與geojson格式轉換
有兩種方法，第一種是用arcgismap自帶的toolbox裡的工具，路徑為 system toolboxes conversion tools json json to features與features to json。這裡shp轉json一般不會報錯。主要說明一下json轉shp。1.首先你的...
相交 光線與球
光線與球之間的相交很容易計算。如果光線的兩端分別是 x1,y1,z1 和 x2,y2,z2 則第一步是將光線引數化 x x1 x2 x1 t x1 it y y1 y2 y1 t y1 jt z z1 z2 z1 t z1 kt 其中0 t 1 中心在 l,m,n 半徑為r的球由下式給出 x l 2...
IDL實現向量 shp 裁剪柵格TASK（一
隨著envi idl版本的更新，idl對向量和柵格資料的處理也變得越來越簡單化。其提供了很多方便的介面，使得使用者呼叫和學習練習便捷成為了可能。最近接觸idl，發現好多網上的 都是延後的，新的介面 理解和編寫起來都比較方便，尤其是在做大量資料研究和應用時，使用批處理的方式顯得尤其重要。新的介面還在摸...