(转载)利用现有shapefile数据提取轮廓线

(转载)利用现有shapefile数据提取轮廓线
2008-01-08 12:37一种利用现有shapefile数据提取轮廓线的方法，并可根据所设定的范围进行截取。

最近碰到两个问题，一个是我需要显示南中国海地区的矢量数据，但我又不想把全世界的海岸线数据都读进来，既浪费内存又浪费时间；另一个问题是我需要中国地区的国界线数据，而不要省界数据，找了好久都没找到这种shapefile数据。

所以就只有自己动手了，利用cntry02.shp(世界政区)和
china.shp(中国省界)两个文件，和利用shapefile数据提取轮廓线的方法，提取出了南中国海和中国国界线的数据，并保存为了shapefile文件，以后就可以随便用
了。

;======================================== ==================
pro buildshp, maskimage, outlineshpfile, needrange
;trick
maskimage =
congrid(maskimage,800,800,/interp,/center)
window, /pixmap, xsize=800, ysize=800
contour, maskimage, path_xy=xy,
path_info=info, $
nlevels=1, closed=1, $
xmargin=[0,0], ymargin=[0,0], xstyle=4, ystyle=4
isoShp = obj_new('IDLffShape', outlineshpfile, /update, entity_type=5)
;
; 添加属性
;
isoShp->IDLffShape::AddAttribute, 'LEVEL', 3, 20
value = info.value
uniqV = value[uniq(value, sort(value))]
for i=0, n_elements(UniqV)-1 do begin
Index = where(info.VALUE eq UniqV[i])
x = 0.0
y = 0.0
n_parts = n_elements(Index)
parts = 0
n_vertices = 0
for j=0, n_parts-1 do begin
Pos = index[j]
LL = double(xy[*,
info[Pos].offset:(info[Pos].offset+info[Pos].N-1)])
LL[0,*] = LL[0,*] *
(needrange[2]-needrange[0]) + needrange[0]
LL[1,*] = LL[1,*] *
(needrange[3]-needrange[1]) + needrange[1]
x = [x, transpose(LL[0,*]), (transpose(LL[0,*]))[0]]
y = [y, transpose(LL[1,*]), (transpose(LL[1,*]))[0]]
parts = [parts,
parts[j]+n_elements(LL[0,*])+1]
n_vertices = n_vertices +
n_elements(LL[0,*]) + 1
endfor
vertices = dblarr(2, n_vertices)
vertices[0,*] = x[1:*]
vertices[1,*] = y[1:*]
parts = parts[0:(n_elements(parts)-2)] ;
; Create structure for new entity.
;
entNew = {IDL_SHAPE_ENTITY}
;
; Define the values for the new entity.
;
entNew.SHAPE_TYPE = 5L
entNew.N_VERTICES = long(n_vertices)
entNew.VERTICES = ptr_new(vertices)
entNew.N_PARTS =
long(n_elements(parts))
entNew.PARTS = ptr_new(parts)
isoShp->IDLffShape::PutEntity, entNew
attrNew =
IsoShp->IDLffShape::GetAttributes(/ATTRIBUTE_STR UCTURE)
attrNew.ATTRIBUTE_0 = long(i)
isoShp->IDLffShape::SetAttributes, i, attrNew
isoShp->DestroyEntity, entNew
endfor
isoShp->Close
obj_destroy, isoShp
end
;=========================================== ===============
function buildmask, shpfile, needrange
;xsize, ysize可用来控制输出的精度,尺寸上限受系统的限制
xsize = (needrange[2]-needrange[0])*100 < 2000
ysize = (needrange[3]-needrange[1])*100 < 2000
; 创建背板，图形在背板中绘制
window, 1, /pixmap, xsize=xsize, ysize=ysize
wset, 1
oshape = obj_new('IDLffShape', shpfile)
oshape->getproperty, n_entities=nentity
for i=0, nentity-1 do begin
ent = oshape->IDLffShape::getentity(i)
if ptr_valid(ent.parts) then begin
cuts = [*ent.parts, ent.n_vertices]
for j=0, ent.n_parts-1 do begin
x =
(*ent.vertices)[0,cuts[j]:cuts[j+1]-1]
y =
(*ent.vertices)[1,cuts[j]:cuts[j+1]-1]
x =
(x-needrange[0])/(needrange[2]-needrange[0])*xsize
y =
(y-needrange[1])/(needrange[3]-needrange[1])*ysize
polyfill, x, y, /device
endfor
endif
oshape->destroyentity, ent
endfor
obj_destroy, oshape
; 拷贝背板中的图像
maskimage = tvrd()
return, maskimage
end
;=========================================== ===============
pro shpoutline
;
;=====STEP 1=====
;在pixmap中绘制所需区域的填充图
;
;源shapefile文件
shpfile =
'C:\RSI\IDL63\resource\maps\shape\cntry02.shp'
;所需区域的范围
needrange = double([105.,0.,122.,23.])
maskimage = buildmask(shpfile, needrange)
;
;=====STEP 2=====
;使用IDL提供的contour方法绘制填充图的等值线，再将等值线保存为shapefile格式。

;
;输出的shapefile文件
outlineshpfile = 'd:\outline.shp'
buildshp, maskimage, outlineshpfile, needrange
end。