| Adding shapefiles to globe model [message #49964] |
Thu, 24 August 2006 10:06 |
Pete Warner
Messages: 14
Registered: July 2006
|
Junior Member |
|
|
I've been moving from direct graphics to object graphics in IDL and
stumbled onto the problem of displaying maps in object graphics. I
mashed out some code to read a shapefile and display the contents on a
globe. With the hope that this helps someone as much as this board has
helped me I am posting my solution. Of course I have the ulterior
motive of having someone find mistakes or suggest improvements.
;###################################################
;function to convert latitude, longitude value
;in degrees to cartesian coordinates on a globe
;radius = the radius of the sphere the points are on
function toglobe, locations, RADIUS=radius
compile_opt idl2
if n_elements(radius) eq 0 then radius = 1.0d else $
radius = double(radius)
size_arr = size( locations )
locations[1,*] = 90.0d - locations[1,*] ;convert latitude to
colatitude
locations[0,*] = 90.0d - locations[0,*]
locations = locations * !DTOR ;convert latitude, longitude to radians
cartesian = radius * $
[cos(locations[0,*])*sin(locations[1,*]), $
cos(locations[1,*]), $
sin(locations[0,*]) * sin(locations[1,*])]
return, cartesian
end
;###############################################
function read_shapefile, file, RADIUS = radius, SKIP=skip
compile_opt idl2
radius = (n_elements(radius) gt 0) ? radius : 1.0
skip = (n_elements(skip) gt 0) ? skip : -1
model = Obj_new('IDLgrModel')
shapefile = Obj_New('IDLffShape', file)
shapefile->GetProperty, N_ENTITIES=num_ent
; Read all the entities.
FOR i=0, (num_ent-1) DO BEGIN
;Read the entity x
ent = shapefile->GetEntity(i)
for j = 0, (ent.n_parts-1) do begin
if ( where( skip eq i) eq -1 ) then begin ;skip selected shapes
index_begin = (*ent.parts)[j]
index_end = (j ne ent.n_parts - 1) ? $
(*ent.parts)[j+1] - 1 : ent.n_vertices-1
current_verts = (*ent.vertices)[*,index_begin:index_end]
points = toglobe( current_verts, RADIUS=radius )
oshape = obj_new('IDLgrPolyline', points)
model -> add, oshape
endif
endfor
shapefile->DestroyEntity, ent ;Clean-up of pointers
ENDFOR
obj_destroy, shapefile
return, model
end
;###############################################
example use:
countryFile = Filepath(Subdirectory=['resource', 'maps', 'shape'], $
'cntry02.shp')
ocountries = read_shapefile(countryFile)
xobjview, ocountries
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
