comp.lang.idl-pvwave archive: archive » cgmap

Home » Public Forums » archive » cgmap_gshhs.pro minarea issue

Show: Today's Messages :: Show Polls :: Message Navigator
E-mail to friend

cgmap_gshhs.pro minarea issue [message #86594]

Thu, 21 November 2013 16:09

pvelissariou
Messages: 4
Registered: November 2013

Junior Member

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)

Report message to a moderator

Re: cgmap_gshhs.pro minarea issue [message #86595 is a reply to message #86594]

Thu, 21 November 2013 16:33

David Fanning
Messages: 11724
Registered: August 2001

Senior Member

pvelissariou@fsu.edu writes:

> 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)

Here is what the creators of the data base say:

"The area of small (< 0.1 km^2) polygons got truncated to 0.
This would cause gshhs to consider them as lines (borders or
rivers) instead of polygons. Furthermore, the areas were
recomputed using the WGS-84 ellipsoid as the previous
area values were based on a spherical calculation. Thanks
to José Luis García Pallero for pointing this out. We now
store the area with a magnitude scale tuned to each polygon."

I don't really know what "tuned to each polygon" means, but I'm not
convinced it means what you seem to think it means.

Does anyone else know anything about this?

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.")

Report message to a moderator

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.

Report message to a moderator

Re: cgmap_gshhs.pro minarea issue [message #86598 is a reply to message #86597]

Thu, 21 November 2013 20:10

David Fanning
Messages: 11724
Registered: August 2001

Senior Member

Panagiotis Velissariou writes:

> 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:

Yes, thanks. I found this independently and have already made the
changes you suggest. I just have to do a bit of testing in the morning,
and I'll check it into the Library. Appreciate your help with this. :-)

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.")

Report message to a moderator

Re: cgmap_gshhs.pro minarea issue [message #86609 is a reply to message #86594]

Fri, 22 November 2013 07:57

Takis.Velissariou
Messages: 1
Registered: November 2013

Junior Member

On Thursday, November 21, 2013 7:09:28 PM UTC-5, pvelis...@fsu.edu 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,

I saw the modified code in gmap_gshhs.pro and I think that:
magnitude = ISHFT(f, -26) AND 1B
should change to:
magnitude = ISHFT(f, -26) AND 255B

also
polygonArea = header.area * 0.1 * 10^magnitude
should change to:
polygonArea = header.area * 0.1 (for version lt 9)
and
polygonArea = header.area * 10^(-magnitude) (for version gt 9)

Report message to a moderator

Re: cgmap_gshhs.pro minarea issue [message #86610 is a reply to message #86609]

Fri, 22 November 2013 08:22

David Fanning
Messages: 11724
Registered: August 2001

Senior Member

Takis.Velissariou@deep-c.org writes:

>
> On Thursday, November 21, 2013 7:09:28 PM UTC-5, pvelis...@fsu.edu 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,
>
> I saw the modified code in gmap_gshhs.pro and I think that:
> magnitude = ISHFT(f, -26) AND 1B
> should change to:
> magnitude = ISHFT(f, -26) AND 255B
>

Well, here is the header I am using to figure this out:

http://www.idlcoyote.com/misc/gshhs.h

I presume the magnitude is the 6th byte in the flag (not the 4th as the
documentation indicates). I see no evidence that we are to use a
negative magnitude power. My way of pulling out the flag value is
consistent with how I find the other byte values in the flag.

> also
> polygonArea = header.area * 0.1 * 10^magnitude
> should change to:
> polygonArea = header.area * 0.1 (for version lt 9)
> and
> polygonArea = header.area * 10^(-magnitude) (for version gt 9)

Since I set magnitude=0 if magnitude is not found in the file, and 10^0
is 1, I think my way of doing this is equivalent to yours.

Do you have a polygon you KNOW the area of that we can test?

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.")

Report message to a moderator

Re: cgmap_gshhs.pro minarea issue [message #86611 is a reply to message #86609]

Fri, 22 November 2013 08:39

David Fanning
Messages: 11724
Registered: August 2001

Senior Member

Takis.Velissariou@deep-c.org writes:

> I saw the modified code in gmap_gshhs.pro and I think that:
> magnitude = ISHFT(f, -26) AND 1B
> should change to:
> magnitude = ISHFT(f, -26) AND 255B

I think I am doing this correctly:

IDL> print, binary(f)
0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 1 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1
IDL> print, binary(ishft(f,-26))
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

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.")

Report message to a moderator

Re: cgmap_gshhs.pro minarea issue [message #86612 is a reply to message #86611]

Fri, 22 November 2013 08:46

David Fanning
Messages: 11724
Registered: August 2001

Senior Member

David Fanning writes:

>
> Takis.Velissariou@deep-c.org writes:
>
>> I saw the modified code in gmap_gshhs.pro and I think that:
>> magnitude = ISHFT(f, -26) AND 1B
>> should change to:
>> magnitude = ISHFT(f, -26) AND 255B
>
> I think I am doing this correctly:
>
> IDL> print, binary(f)
> 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 1 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1
> IDL> print, binary(ishft(f,-26))
> 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

I've changed my mind. I think you are right about this. :-)

It should be AND 255B.

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.")

Report message to a moderator

Re: cgmap_gshhs.pro minarea issue [message #86613 is a reply to message #86609]

Fri, 22 November 2013 08:50

David Fanning
Messages: 11724
Registered: August 2001

Senior Member

Takis.Velissariou@deep-c.org writes:

> also
> polygonArea = header.area * 0.1 * 10^magnitude
> should change to:
> polygonArea = header.area * 0.1 (for version lt 9)
> and
> polygonArea = header.area * 10^(-magnitude) (for version gt 9)

I still don't know what to make of this. The documentation is screwy.

At one place in the header file you find this:

"area magnitude scale p (as in 10^p) = flag >> 26.
We divide area by 10^p."

Then, in the next couple of lines, you find this *twice*:

"Area of polygon in km^2 * 10^p for this resolution file"

We need to know the area of known polygons to be able to tell.

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.")

Report message to a moderator

Re: cgmap_gshhs.pro minarea issue [message #86614 is a reply to message #86594]

Fri, 22 November 2013 09:34

pvelissariou
Messages: 4
Registered: November 2013

Junior Member

On Thursday, November 21, 2013 7:09:28 PM UTC-5, pvelis...@fsu.edu 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,

From gshhs.c file we have:
line = (h.area) ? 0 : 1; /* Either Polygon (0) or Line (1) (if no area) */
if (version >= 9) { /* Variable magnitude for area scale */
m = h.flag >> 26;
scale = pow (10.0, (double)m); /* Area scale */
}
area = h.area / scale; /* Now im km^2 */
f_area = h.area_full / scale; /* Now im km^2 */

so, in the cgmap_gshhs.pro
we should have something like:
magnitude = 10.0 ^ ( - ISHFT(f, -26) ) for version >= 9
and
magnitude = 0.1 for version < 9

and
polygonArea = header.area * magnitude

in order to produce the correct results for the polygon area

Report message to a moderator

Re: cgmap_gshhs.pro minarea issue [message #86615 is a reply to message #86613]

Fri, 22 November 2013 09:37

David Fanning
Messages: 11724
Registered: August 2001

Senior Member

David Fanning writes:

>
> Takis.Velissariou@deep-c.org writes:
>
>> also
>> polygonArea = header.area * 0.1 * 10^magnitude
>> should change to:
>> polygonArea = header.area * 0.1 (for version lt 9)
>> and
>> polygonArea = header.area * 10^(-magnitude) (for version gt 9)
>
> I still don't know what to make of this. The documentation is screwy.
>
> At one place in the header file you find this:
>
> "area magnitude scale p (as in 10^p) = flag >> 26.
> We divide area by 10^p."
>
> Then, in the next couple of lines, you find this *twice*:
>
> "Area of polygon in km^2 * 10^p for this resolution file"
>
> We need to know the area of known polygons to be able to tell.

OK, I agree with you about dividing by the magnitude, too. :-)

I have an example using the Australian polygon, and I get the correct
area for the land mass only by dividing: header.area / 10^magnitude.

I'm currently having trouble with v. 2.1 of the GSHHS data. Clearly the
header is different from the one for v 2.0 and 2.2, but I can't figure
out what it is. :-(

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.")

Report message to a moderator

Re: cgmap_gshhs.pro minarea issue [message #86621 is a reply to message #86594]

Sat, 23 November 2013 20:03

pvelissariou
Messages: 4
Registered: November 2013

Junior Member

On Thursday, November 21, 2013 7:09:28 PM UTC-5, pvelis...@fsu.edu 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,

I think the following changes are required:

the line: river = ISHFT(f, -25) AND 255B
should be: river = ISHFT(f, -25) AND 1B

also, header.area is of type LONG and 10^magnitude is of type INT
therefore the calculated polygonarea is wrong, I have modified the
following lines as:

from:
IF version LE 8 THEN BEGIN
polygonArea = header.area * 0.1 ; km^2
ENDIF ELSE BEGIN
polygonArea = header.area / 10^magnitude ; km^2
ENDELSE
to:
IF version LT 9 THEN BEGIN
polygonArea = Double(header.area) * 0.1 ; km^2
ENDIF ELSE BEGIN
polygonArea = Double(header.area) / 10.0^magnitude ; km^2
ENDELSE

Report message to a moderator

Re: cgmap_gshhs.pro minarea issue [message #86622 is a reply to message #86621]

Sat, 23 November 2013 21:15

David Fanning
Messages: 11724
Registered: August 2001

Senior Member

pvelissariou@fsu.edu writes:

> I think the following changes are required:
>
> the line: river = ISHFT(f, -25) AND 255B
> should be: river = ISHFT(f, -25) AND 1B
>
> also, header.area is of type LONG and 10^magnitude is of type INT
> therefore the calculated polygonarea is wrong, I have modified the
> following lines as:
>
> from:
> IF version LE 8 THEN BEGIN
> polygonArea = header.area * 0.1 ; km^2
> ENDIF ELSE BEGIN
> polygonArea = header.area / 10^magnitude ; km^2
> ENDELSE
> to:
> IF version LT 9 THEN BEGIN
> polygonArea = Double(header.area) * 0.1 ; km^2
> ENDIF ELSE BEGIN
> polygonArea = Double(header.area) / 10.0^magnitude ; km^2
> ENDELSE

OK, I'm going to assume you can make more sense of that documentation
than I can. ;-)

I'll check this in tomorrow with the other changes. Thanks for your
help.

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.")

Report message to a moderator

Previous Topic:	IDL 8.2, read pixel value along a surface
Next Topic:	Exelis code library down - someone has G

-=] Back to Top [=-

[ Syndicate this forum (XML) ] [

] [

]

Current Time: Wed Oct 08 15:37:26 PDT 2025

Total time taken to generate the page: 0.00693 seconds