| Trouble plotting GSHHS shapefile [message #88150] |
Mon, 24 March 2014 14:46  |
timothyja123
Messages: 57
Registered: February 2013
|
Member |
|
|
Hi Guys,
I'm having trouble plotting the GSHHS shapefile as found here: http://www.ngdc.noaa.gov/mgg/shorelines/data/gshhg/latest/
I believe its go something to do with the way I'm using mapset but I cant figure it out basically my code seems to have trouble drawing the antarctic region. The result also means I cant use polyfill to colour the land as it just draws over everything due to the same issues with plot.
I've replaced my plotting code with a call to the coyote drawshapes routine but it has the same issue. Anyway here is my code any help is greatly appreaciated.
pro plot_gshhs_shapefile
limit = [-90, 0, 90, 360]
POLON = 0; limit[1] + 0.5*(limit[3] - limit[1])
POLAT = 0 ; default
map_set, POLAT, POLON, 0, /cylindrical, ISOTROPIC=1, LIMIT=limit, noborder=1, noerase=1, xmargin = 0, ymargin = 0
Stack = Scope_Traceback(/Structure)
Filename = Stack[N_elements(Stack) - 2].Filename
If (Arg_Present(BaseName)) then Begin
BaseName = File_BaseName(Filename)
EndIf
sourcePath = File_DirName(Filename, _Extra = Extra)
colour = '476A34'
file = filepath('GSHHS_i_L1.shp', root_dir=sourcePath, subdir=['data','gshhg_shapefiles','GSHHS_shp','i'])
;drawshapes, file
oGshhsShape = obj_new('idlffshape', file)
oGshhsShape -> getproperty, n_entities = num_ent
entities = oGshhsShape -> GetEntity(/All, /Attributes)
; Add the data from the shape file.
FOR i = 0L, N_Elements(entities)-1 do begin
entity = entities[i]
lats = reform((*(entity.vertices))[1,*])
lons = reform((*(entity.vertices))[0,*])
parts = *entity.parts
FOR iPart = 0, n_elements(parts)-1 DO BEGIN
startVertIdx = parts[iPart]
endVertIdx = (iPart EQ entity.n_parts-1) ? (N_ELEMENTS(lats)-1) : (parts[iPart+1]-1)
part_lons = lons[startVertIdx:endVertIdx]
part_lats = lats[startVertIdx:endVertIdx]
;polyfill, part_lons, part_lats, /data, color=colour
plots, part_lons, part_lats, /data
ENDFOR
PTR_FREE, entity.vertices
PTR_FREE, entity.measure
PTR_FREE, entity.parts
PTR_FREE, entity.part_types
PTR_FREE, entity.attributes
ENDFOR
end
|
|
|
|
| Re: Trouble plotting GSHHS shapefile [message #88151 is a reply to message #88150] |
Mon, 24 March 2014 19:28  |
timothyja123
Messages: 57
Registered: February 2013
|
Member |
|
|
ok so it turned out I needed to use map_proj_init() to use the correct projection but now how do I move the map so its centered in a different location like I can do with map_set?
I tried setting the center_longitude param but then again things are not plotted correctly.
pro plot_gshhs_shapefile
limit = [-90, 0, 90, 360]
POLON = limit[1] + 0.5*(limit[3] - limit[1])
POLAT = 0 ; default
junk = MAP_PROJ_INIT( 'Cylindrical', ELLIPSOID=24, CENTER_LONGITUDE=POLON)
uv_box = junk.uv_box
Plot, uv_box[[0, 2]], uv_box[[1, 3]], Position=[0.1,0.1,0.9,0.9], $
/Nodata, XStyle=5, YStyle=5, /NoErase
Stack = Scope_Traceback(/Structure)
Filename = Stack[N_elements(Stack) - 2].Filename
If (Arg_Present(BaseName)) then Begin
BaseName = File_BaseName(Filename)
EndIf
sourcePath = File_DirName(Filename, _Extra = Extra)
colour = '476A34'
file = filepath('GSHHS_i_L1.shp', root_dir=sourcePath, subdir=['data','gshhg_shapefiles','GSHHS_shp','i'])
oGshhsShape = obj_new('idlffshape', file)
oGshhsShape -> getproperty, n_entities = num_ent
entities = oGshhsShape -> GetEntity(/All, /Attributes)
; Add the data from the shape file.
FOR i = 0L, N_Elements(entities)-1 do begin
entity = entities[i]
lats = reform((*(entity.vertices))[1,*])
lons = reform((*(entity.vertices))[0,*])
parts = *entity.parts
FOR iPart = 0, n_elements(parts)-1 DO BEGIN
startVertIdx = parts[iPart]
endVertIdx = (iPart EQ entity.n_parts-1) ? (N_ELEMENTS(lats)-1) : (parts[iPart+1]-1)
part_lons = lons[startVertIdx:endVertIdx]
part_lats = lats[startVertIdx:endVertIdx]
result = MAP_PROJ_FORWARD(part_lons, part_lats, MAP_STRUCTURE=junk)
polyfill, result, color=colour
plots, result
ENDFOR
PTR_FREE, entity.vertices
PTR_FREE, entity.measure
PTR_FREE, entity.parts
PTR_FREE, entity.part_types
PTR_FREE, entity.attributes
ENDFOR
end
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
