一种利用现有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 =
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_STRUCTURE)
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