| Re: cgmap_gshhs.pro minarea issue [message #86597 is a reply to message #86594] |
Thu, 21 November 2013 18:05   |
pvelissariou
Messages: 4
Registered: November 2013
|
Junior Member |
|
|
On Thursday, November 21, 2013 7:09:28 PM UTC-5, Panagiotis Velissariou wrote:
> Apparently, in the recent versions (>= 2.2) of gshhs database
>
> the units of the header.area changed from 1/10 km^2 to 1/10 m^2.
>
> For cgmap_gshhs to work properly the line:
>
> polygonArea = header.area * 0.1 (ok for gshhs < 2.2)
>
> should be changed to:
>
> polygonArea = header.area * 1.0e-7 (for gshhs >= 2.2)
David,
Thank you for the reply. You are right.
The problem is that from version 2.2 and on they have introduced a
magnification factor for the area, see the header structure below:
struct GSHHS { /* Global Self-consistent Hierarchical High-resolution Shorelines */
int id; /* Unique polygon id number, starting at 0 */
int n; /* Number of points in this polygon */
int flag; /* = level + version << 8 + greenwich << 16 + source << 24 + river << 25 + p << 26 */
/* flag contains 6 items, as follows:
* low byte: level = flag & 255: Values: 1 land, 2 lake, 3 island_in_lake, 4 pond_in_island_in_lake
* 2nd byte: version = (flag >> 8) & 255: Values: Should be 9 for GSHHS release 9
* 3rd byte: greenwich = (flag >> 16) & 3: Values: 0 if Greenwich nor Dateline are crossed,
* 1 if Greenwich is crossed, 2 if Dateline is crossed, 3 if both is crossed.
* 4th byte: source = (flag >> 24) & 1: Values: 0 = CIA WDBII, 1 = WVS
* 4th byte: river = (flag >> 25) & 1: Values: 0 = not set, 1 = river-lake and GSHHS level = 2 (or WDBII class 0)
* 4th byte: area magnitude scale p (as in 10^p) = flag >> 26. We divide area by 10^p.
*/
int west, east, south, north; /* min/max extent in micro-degrees */
int area; /* Area of polygon in km^2 * 10^p for this resolution file */
int area_full; /* Area of corresponding full-resolution polygon in km^2 * 10^p */
int container; /* Id of container polygon that encloses this polygon (-1 if none) */
int ancestor; /* Id of ancestor polygon in the full resolution set that was the source of this polygon (-1 if none) */
};
I have modified the code in cgmap_gshhs.pro as follows:
; Discriminate polygons based on header information.
IF version LT 9 THEN BEGIN
area_mag = 1.0e-1 ; km^2 / 10 -> km^2
ENDIF ELSE BEGIN
area_mag = 10.0 ^ ( - ISHFT(f, -26) ); km^2 / 10^p -> km^2
ENDELSE
polygonArea = header.area * area_mag
and it seems that gives the correct results.
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
