CGMAP_GSHHS: problem with land/water colors [message #85515] |
Wed, 14 August 2013 12:08  |
Matteo
Messages: 28 Registered: August 2011
|
Junior Member |
|
|
Hi,
I noticed using CGMAP_GSHHS.PRO that if no ocean waters are present on a map, the LAND keyword for the land color does not work properly. Here's a snippet code that reproduces the problem (note that it is very similar to the one included in the examples of the CGMAP_GSHHS header). Just set your PATH to include the CGxxx routines and the correct directory for the gshhs.b file.
PRO test_gshhs
!path=[SET PATH HERE]
datafile='./gshhs_h.b' [SET PATH TO GSHHS_H.B HERE]
; include some water
map_limits = [29.,-95.,34.,-88]
; land only
;map_limits = [31.,-95.,34.,-88]
cgDisplay, 500, 350
pos = [0.1,0.1, 0.9, 0.8]
; set map projection
cgmap_set, limit = map_limits, /mercator, position=pos
; set a SKYBLUE background
cgColorfill, [pos[0], pos[0], pos[2], pos[2], pos[0]], $
[pos[1], pos[3], pos[3], pos[1], pos[1]], $
/Normal, Color='skyblue'
; issue CGMAP_GSHHS
cgMap_GSHHS, datafile, Fill=1, Level=4, Color='black', /Outline, $
Land='tan', Water_color='skyblue', NoClip=0
; overdraw state borders
cgmap_set, limit = map_limits, /mercator, /continents, position=pos, $
color='yellow', /usa, /noerase
END
Run the program first with map_limits set to include some Gulf waters around Lousiana (map_limits = [29.,-95.,34.,-88]). Then comment in the 'land only' option to exclude the water at the lower latitudes.
In the latter case I obtain a SKYBLUE rather than a TAN background, which is obviously wrong. Furthermore, the WATER_COLOR keyword only works for inland water bodies, therefore CGCOLORFILL must be used to start from a SKYBLUE background.
You will also notice that having to reissue CGMAP_SET at the end to draw state borders (yellow) over everything leads to imperfect overlap with the black GSHHS shorelines, expected because of the different resolution.
I would appreciate if anybody could take a look at it. Apologies to David if I'm mistaking the correct usage of his routine.
m
|
|
|
|
Re: CGMAP_GSHHS: problem with land/water colors [message #85520 is a reply to message #85515] |
Wed, 14 August 2013 12:38   |
David Fanning
Messages: 11724 Registered: August 2001
|
Senior Member |
|
|
Matteo writes:
> I noticed using CGMAP_GSHHS.PRO that if no ocean waters are present on a map, the LAND keyword for the land color does not work properly. Here's a snippet code that reproduces the problem (note that it is very similar to the one included in the examples of the CGMAP_GSHHS header). Just set your PATH to include the CGxxx routines and the correct directory for the gshhs.b file.
Well, it appears your observation is correct, but your explanation is
wrong. :-)
I have a bit of code in here to speed things up:
xy = Convert_Coord(lon, lat, /Data, /To_Normal)
check = ((xy[0,*] GE !X.Window[0]) AND (xy[0,*] LE !X.Window[1])) AND $
((xy[1,*] GE !Y.Window[0]) AND (xy[1,*] LE !Y.Window[1]))
usePolygon = Max(check)
Your particular example appears to be a perverse one, because no land
polygons ever set the usePolygon flag to 1, so land polygons are not
filled.
If you set usePolygon=1 just after this code and run your program, you
will see it does the right thing, but takes a LONG time to run. (I used
the intermediate resolution file. You seem to be using the high
resolution file, which I wouldn't recommend. That will take an EXTREMELY
long time to run!)
So, the problem is to understand why these land polygons aren't being
drawn. I presume it is because they have one or more vertices outside
the map boundary. I'll have to see what I can do about this. I don't
know how many polygons there are, but the term "millions" comes to mind.
I'll see if I can find the 3-4 I am looking for. It may be awhile before
I get back in touch. :-(
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.idlcoyote.com/
Sepore ma de ni thue. ("Perhaps thou speakest truth.")
|
|
|
Re: CGMAP_GSHHS: problem with land/water colors [message #85522 is a reply to message #85515] |
Wed, 14 August 2013 13:08   |
David Fanning
Messages: 11724 Registered: August 2001
|
Senior Member |
|
|
Matteo writes:
> Run the program first with map_limits set to include some Gulf
> waters around Lousiana (map_limits = [29.,-95.,34.,-88]). Then
> comment in the 'land only' option to exclude the water at the
> lower latitudes. In the latter case I obtain a SKYBLUE rather
> than a TAN background, which is obviously wrong. Furthermore,
> the WATER_COLOR keyword only works for inland water bodies,
> therefore CGCOLORFILL must be used to start from a SKYBLUE background.
OK, here is the problem. To speed things up I assume that I don't want
to draw polygons that are outside the plot area. So, I exclude any
polygon in which all of its vertices are outside the drawing area. The
polygon we want to draw in this zoomed-in case is, say, the United
States polygon. But, all of its vertices are outside the drawing area,
so it never gets drawn.
In this case, you can "fix" the problem by simply drawing with a tan
color instead of a skyblue color before drawing the inland water
polygons. In general, it would probably be a better idea to NOT use the
GSHHS data base for drawing such zoomed in plots. I think I would have
used cgMap_Continents instead, and used the GSHHS data for just the
inland waterways, if I needed that much detail. In other words, I think
this is a wrong tool for the job kind of problem. (Unless, of course,
you have a coffee machine nearby and you don't really care now long it
all takes.)
You have to "fill" oceans with the fill color before you start to draw
because the GSHHS data doesn't include ocean polygons. In other words,
the oceans are what's left (the negative space) after you fill all the
other polygons in the file.
> You will also notice that having to reissue CGMAP_SET at the end to
> draw state borders (yellow) over everything leads to imperfect
> overlap with the black GSHHS shorelines, expected because of
> the different resolution.
You don't have to do this, and I don't know why you think so. I would
just use cgMap_Continents and just draw the state outlines:
cgMap_Continents, /usa, color='yellow'
Here is the code I used to draw the "correct" plot in this particular
instance.
PRO test_gshhs
datafile='C:\IDL\data\gshhs\gshhs_2.2\gshhs\gshhs_i.b'
; include some water
;map_limits = [29.,-95.,34.,-88]
; land only
map_limits = [31.,-95.,34.,-88]
cgDisplay, 500, 350, /free
pos = [0.1,0.1, 0.9, 0.8]
; set map projection
cgmap_set, limit = map_limits, /mercator, position=pos
; set a SKYBLUE background
cgColorfill, [pos[0], pos[0], pos[2], pos[2], pos[0]], $
[pos[1], pos[3], pos[3], pos[1], pos[1]], $
/Normal, Color='tan'
; issue CGMAP_GSHHS
cgMap_GSHHS, datafile, Fill=1, Level=4, Color='black', $
Land='tan', Water_color='skyblue'
; overdraw state borders
cgMap_Continents, /usa, color='yellow'
END
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.idlcoyote.com/
Sepore ma de ni thue. ("Perhaps thou speakest truth.")
|
|
|
Re: CGMAP_GSHHS: problem with land/water colors [message #85523 is a reply to message #85522] |
Wed, 14 August 2013 13:33   |
David Fanning
Messages: 11724 Registered: August 2001
|
Senior Member |
|
|
David Fanning writes:
> I think I would have
> used cgMap_Continents instead, and used the GSHHS data for just the
> inland waterways, if I needed that much detail. In other words, I think
> this is a wrong tool for the job kind of problem. (Unless, of course,
> you have a coffee machine nearby and you don't really care now long it
> all takes.)
Here is how I would have written this code. It works with either of your
two map limits.
PRO test_gshhs
datafile='C:\IDL\data\gshhs\gshhs_2.2\gshhs\gshhs_i.b'
; include some water
map_limits = [29.,-95.,34.,-88]
; land only
;map_limits = [31.,-95.,34.,-88]
cgDisplay, 500, 350, /free
pos = [0.1,0.1, 0.9, 0.8]
; set map projection
cgmap_set, limit = map_limits, /mercator, position=pos
cgColorFill, Position=pos, Color='tg5'
cgMap_Continents, color='tan', /continents, /fill
; issue CGMAP_GSHHS
cgMap_GSHHS, datafile, Fill=1, Level=2, Color='black', $
Water_color='blu4'
; overdraw state borders
cgMap_Continents, /usa, color='yellow'
cgPolygon, Position=pos, Color='black'
END
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.idlcoyote.com/
Sepore ma de ni thue. ("Perhaps thou speakest truth.")
|
|
|
|
Re: CGMAP_GSHHS: problem with land/water colors [message #85526 is a reply to message #85525] |
Wed, 14 August 2013 15:41  |
David Fanning
Messages: 11724 Registered: August 2001
|
Senior Member |
|
|
Matteo writes:
> thank you very much. Every time I tried to respond I was being anticipated ;-). I realized I could start off a TAN background but the fact is that I need to batch process a lot of files containing flight transect to be overlaid to these maps (the map boundaries are automatically calculated). So I needed the code to just "know" where land and water are. Another solution could be using a land/water flag and switch CGCOLORFILL between TAN and SKYBLUE after checking if
there's any ocean pixel within the map boundaries. I guess it'd be similar in concept to your Convert_Coord solution. It is a bit annoying to see racial discrimination among inland water bodies (try map_limits = [29.,-125.,44.,-105], but it doesn't matter at this point.
>
Here is what I've finally come up with that I like quite a lot. It seems
to work correctly for all of your map limits.
PRO test_gshhs
datafile='C:\IDL\data\gshhs\gshhs_2.2\gshhs\gshhs_i.b'
; include some water
;map_limits = [29.,-95.,34.,-88]
; land only
;map_limits = [31.,-95.,34.,-88]
map_limits = [29.,-125.,44.,-105]
cgDisplay, 500, 350, /free
pos = [0.1,0.1, 0.9, 0.8]
; set map projection
cgmap_set, limit = map_limits, /mercator, position=pos
cgColorFill, Position=pos, Color='tg5'
cgMap_Continents, color='tan', /continents, /fill, /hires
; issue CGMAP_GSHHS
cgMap_GSHHS, datafile, Fill=1, Level=2, Color='black', $
Water_color='blu4'
; overdraw state borders
cgMap_Continents, /usa, color='yellow'
;cgMap_Continents, /continents, color='black', /hires
cgMap_GSHHS, datafile, color='black'
cgPolygon, Position=pos, Color='black'
END
Cheers,
Dvaid
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.idlcoyote.com/
Sepore ma de ni thue. ("Perhaps thou speakest truth.")
|
|
|