(转载)利用现有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 =