Re: inverse mapproject -C confusion

Hi Joaquim,

I'm probably still confused about the projections, but I don't understand the different behavior (-C and -R) of the inverse projection between -Je and -Jm. The inverse projection appeared to work as I expected using the Azimuth Equidistance (-Je) projection with "-C", but, with the same set Cartesian (km) points, the Mercator (-Jm) inverse projection only works without "-C" and with a modified range.

I'm trying to create a geographic grid covering a 100-km radius circle centered at -122,38 with 1-km spacing. Is there a better way to generate a geographic grid (data table) with arbitrary grid spacing specified in km?

I started by creating x-y points for a 150x150 (km) grid using perl (example uses 10-km spacing).

perl -e 'for $x (-15..15) {for $y (-15..15) {print $x*10," ",$y*10,"\n";}}' > grid_pts_km.xy

I expected the -C option to put 0,0 km at the projection center (-122,38). This seems to have worked with -Je, range set to encompass expected range of inverse projected coordinates:

mapproject grid_pts_km.xy -Fk -I -C -Je-122/38/1:1 -R-123.75/-120.25/36.5/39.5 > grid_pts_ll_Je.xy;

To get this to work with the -Jm inverse projection, I had to omit the -C option and change the range so that the West and South values were the projection center coordinates (i.e. I wanted 0,0 km to be -122,38 degrees):

mapproject grid_pts_km.xy -Fk -I    -Jm-122/38/1:1 -R-122.00/-120.25/38.0/39.5 > grid_pts_ll_Jm.xy;

The attached map compares the two grids.

Thank you,

Jason

From: Joaquim Luis [mailto:[email protected]]
Sent: Tuesday, January 14, 2014 7:00 PM
To: World-Wide GMT Usage and Help Mailing List
Cc: Altekruse, Jason [FWLA]
Subject: Re: [GMT-HELP] inverse mapproject -C confusion

Jason,

How were your x,y (cartesian) coordinates obtained? What you need to do is to inverse project using the same projection used to by x,y in first place.

mapproject -I -Jm... inverse project coordinates that are in Mercator projection, not any km(etric) units

Joaquim


Hi,

I'm trying to convert a table of x and y coordinates in km to geographic coordinates using mapproject -I. My grid of points is centered at W122 N38 (0,0). If I use an Azimuthal Equidistance projection with the -C option, (0,0) in km projects to (-122,38) as I expected. When I try a Mercator projection with the -C option, x projects to -122 but y is projected to 0 degrees latitude. I have to exclude the -C option and adjust the Range so that the lower left corner is where I want the center of the grid (-122,38).

Inverse projection (expect: 0,0 km --> -122,38 degrees)
$ echo 0 0 | mapproject -I -Fk -C -Je-122/38/1:1 -R-123.5/-120.5/36.5/39.5
-122   38

$ echo 0 0 | mapproject -I -Fk -C -Jm-122/38/1:1 -R-123.5/-120.5/36.5/39.5
-122   0

$ echo 0 0 | mapproject -I -Fk    -Jm-122/38/1:1 -R-123.5/-120.5/36.5/39.5
-123.5 36.5

$ echo 0 0 | mapproject -I -Fk    -Jm-122/38/1:1 -R-122/-120.5/38/39.5
-122   38

The forward projection behaves similarly (expect: -122,38 degrees --> 0,0 km):
$ echo -122 38 | mapproject    -Fk -C -Je-122/38/1:1 -R-123.5/-120.5/36.5/39.5
0      0
...
$ echo -122 38 | mapproject    -Fk    -Jm-122/38/1:1 -R-122/-120.5/38/39.5

0                     0

Am I using the -C (or -Fk, -J, -R,) option incorrectly, or misunderstanding the projections (or making a poor choice)? Is there a better way to create a geographic grid with 1 km spacing?

I'm using GMT 5.1.1_r12604 on a virtual Ubuntu machine - I'm visiting another office and haven't managed to navigate svn through their proxy to update.

I'm ultimately trying to create a 100x100 km grid of points with 1 km spacing in geographic coordinates for use in a seismic hazard calculation, and will eventually calculate rate, beta, and cell area (in units of equatorial square degrees) at each point. I started by generating the x,y table with awk/perl.

Thanks for your help,

Jason


Jason Altekruse
Project Geologist

Fugro Consultants, Inc.
Lakewood, CO 80401

To unsubscribe, send the message "signoff gmt-help" to [email protected]<mailto:[email protected]>


To unsubscribe, send the message "signoff gmt-help" to [email protected]