Using clipping in GPROJECT resulting in rotated maps

[DaveW]

Hi.

I'm trying to use the clipping functionality in PROC GPROJECT but when I do the resulting map (created using PROC GMAP) is rotated counter clockwise.  The amount of rotation appears to be proportional to the amount of clipping but occurs even when I set the bounds beyond the latitude and longitude that occurs in the un-projected file.

I've run the example found in the documentation (SAS/GRAPH(R) 9.2: Reference, Second Edition) and it works with no issues.

Of course the simplest solution is to just sub-set the boundary file before I project it but that loses some of the nice properties of clipping.

Has anyone seen this before?

David

[GraphGuy]

Are you using one of the maps we ship with SAS?

If so, can you show us the code you're using?

Is your map in westlong radians? If not, are you using the 'eastlong' and 'degrees' options?

Make sure you don't have the longmin and longmax reversed.

What projection are you using? Have you tried it with project=none ?

[DaveW]

Hi Robert.

No, I'm not using one of the maps that ships with SAS but rather ESRI boundary files (which I know are in lat/long) supplied by a third party, an example of which can be found here:  http://www.abs.gov.au/ausstats/subscriber.nsf/log?openagent&1270055001_sa4_2011_aust_shape.zip&1270.....

I am using the EAST and DEGREES options as the longitude values are in degrees and increase to the east.  I also don't have the LONGMIN and LONGMAX reversed (SAS gives an error if this is the case).

The map produced by PROC GMAP looks fine if I don't clip (I.E. exclude LONGMIN and LONGMAX from the PROC GPROJECT statement) for both the Albers' and gnomonic projections (for some reason Lambert's projection rotates the whole map 180 degrees but that's an issue for another day).  And the final map looks fine (apart from the east-west stretch of an unprojected map) if I use the PROJECT=NONE option with clipping.  It's when I try to use both that things get interesting:

• when using Albers' projection the map rotates counterclockwise (see image);
• when using the gnomonic projection I get a boundary file with just 4 rows with missing X or Y data (and no warning or error in the log)!?!

This code...

gproject data=sa4_bound out=crop_1 project=albers east degrees; id sa4_code11; run;

results in

This code...

gproject data=sa4_bound out=crop_1 project=albers east degrees longmin=106 longmax=159; id sa4_code11; run;

results in

Finally if I use no projection, this code

gproject data=sa4_bound out=crop_1 project=none east degrees longmin=106 longmax=159; id sa4_code11; run;

results in

[Darrell_sas]

It appears rotated in 9.3 and does not in 9,4.  I'll have to investigate further, but it may be a defect that was fixed.

[GraphGuy]

Ahh - that might explain it! There have been several fixes and enhancements to proc gproject in 9.4.

[Darrell_sas]

Remove the "east' to fix the problem.

First, turn the x to "west".

    data sa4_bound;

         set sa4_bound;

         x= -x;

   run;

Then, remove the "east" from the GPROJECT.

proc gproject data=sa4_bound out=crop_1 project=albers /*east*/ degrees

     longmin=106 longmax=159; id sa4_code11; run;

It will now work.

[DaveW]

Hi.

Thanks for your response.  You are right that I'm using 9.3 but unfortunately the suggested fix isn't working for me.  Running the following code I'm still getting the rotation...

sa4_bound_2;

     set sa4_bound;

     x = -x;

run;

gproject data=sa4_bound_2 out=crop_1 project=albers degrees longmin=-159 longmax=-106;

     id sa4_code11;

run;

gmap map=crop_1 data=crop_1 all;

     id sa4_code11;

     choro sa4_code11 / levels=1 woutline=1 nolegend;

run;

[DaveW]

And if I clip the altered boundary file with PROJECT=NONE I get a mirror image.

[Darrell_sas]

You changed

longmin=106 longmax=159

to

longmin=-159 longmax=-106

[Darrell_sas]

Oh, and it is backwards because you reversed it with the x=-x and you are NOT projecting it with PROJECT=NONE.

[DaveW]

OK, so the option change from EAST to WEST (or rather allowing the default of WEST to apply) doesn't account for the change in X values if there is no projection?

[DaveW]

The maximum and minimum longitude have to be changed because the values on the boundary file have changed.  If I don't make this change then no clipping occurs and I get a message saying that the minimum value was larger than the largest value on the file.

[Darrell_sas · Accepted solution]

Then split the Proc GPROJECT into two pieces.   And it is a good idea to reduce the map size with Proc GREDUCE.

proc mapimport out=sa4_bound

infile='C:\Public\Userdata\aust_shape\SA4_2011_AUST.shp';

run;

proc gproject data=sa4_bound out=crop_2 project=none degrees east 

longmin=106 longmax=159;
id sa4_code11;

run;

proc gproject data=crop_2 out=crop_1 project=albers degrees east;

id sa4_code11;

run;

proc greduce data=crop_1 out=crop_1; 

id sa4_code11;

run;

proc gmap map=crop_1(where=(density<3)) data=crop_1 all;

id sa4_code11;

choro sa4_code11 / levels=1 woutline=1 nolegend;

run;

[DaveW]

That works... thanks for your help!

Now I'll just have to hang out for our move to 9.4 for the ultimate fix.