利用现有shapefile数据提取轮廓线

/ 0评 / 0

一种利用现有shapefile数据提取轮廓线的方法,并可根据所设定的范围进行截取。最近碰到两个问题,一个是我需要显示南中国海地区的矢量数据,但我又不想把全世界的海岸线数据都读进来,既浪费内存又浪费时间;另一个问题是我需要中国地区的国界线数据,而不要省界数据,找了好久都没找到这种shapefile数据。所以就只有自己动手了,利用cntry02.shp(世界政区)和china.shp(中国省界)两个文件,和利用shapefile数据提取轮廓线的方法,提取出了南中国海和中国国界线的数据,并保存为了shapefile文件,以后就可以随便用了。

75e1170fb011ef9634d1222b

 

;==========================================================

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

 

发表评论

电子邮件地址不会被公开。 必填项已用*标注