Re-project Lambert projection from imported shapefile

[levender8622]

I downloaded a shapefile (boundary file) from Statistics Canada website and try to map the province of Alberta. However the map is sort of rotated and I want to make it horizontal. The original projection in the shapefile according to Stats Can is Lambert conformal conic. I tried to use proc gproject to re-project but failed. How can I rotate the map so that the mapped province can have a horizontal display? Thanks!

[ballardw]

By "horizontal display" I am going to guess that you mean to the have the east-west southern border running roughly left-right on the display. Is that correct?

Do you have latitude and longitude in your data set?

Can you show the code you attempted and any error messages from the log?

[levender8622]

Thanks for your reply. Sorry that I didn't explain quite clearly. I don't have latitude and longitude in my dataset, I have x and y which seem to be projected lat and long already. My data look like this:

X                     Y                          SEGMENT    CDNAME              CDTYPE     CDUID    PRNAME    PRUID

4864674.58    1822138.7686     1                    Division No.  1      CDR           4801        Alberta       48

where x ranges from 4427326.16 to 5229748.67, and y ranges from 1616915.88 to 2965270.49. My codes are fairly simple and straightforward:

proc mapimport datafile="lcd_000b16a_e.shp" out=mymap; run;

data ab;
   set mymap;
   where PRNAME="Alberta";
run;

proc gmap data=ab map=ab;
   id CDUID;
   choro CDNAME / statistic=first;
run;

There's no error msg and the map looks like attachment 1. However what I'm looking for is a map like attachment 2 (I got this by projecting in another software and import the data to SAS). I tried gproject, hoping that I could project into UTM N18 using the following code, but the resulting map still looks like attachment 1...

data test;
   set ab;
   rename x=LONG y=LAT;
run;

proc gproject latlon
   project=proj4 to="EPSG:32618"
   data=test out=want;
   id CDUID;
run;

[ballardw]

From the documentation:

Note:   Projection is appropriate for map data sets in which the X and Y variable values represent longitude and latitude. Some of the map data sets that are supplied with SAS/GRAPH have already been projected; such data set should not be projected again. 

Your shapefile after import is going to be similar to one of the SAS supplied map sets without longitude and latitude. 

You may want to investigate the map libraries supplied with SAS. There may well be a canada map with what you are using. The version I have seems to be a bit outdate but Canada3 comes close to having what you show though it is unprojected and hence east-west looks funny. The Canada3 set I have has "Census Districts" represented by CDCODE variable if that is what you are looking for.

Though since it looks like you already had a successful import of another map set I am a bit confused why you wouldn't use the one you imported from the other source.

[GraphGuy]

I believe SAS' Proc Gproject can only work on unprojected lat/long values.

If your map is already projected, I don't think you'll be able to unproject & reproject it in SAS.

I would recommend using the new unprojected GfK maps we ship with SAS/Graph. These have both a projected X/Y, and an unprojected lat/long. You could create the desired Alberta map as follows:

data alberta_map; set mapsgfk.canada (where=(ID1='CA-48'));
run;

proc gproject data=alberta_map out=alberta_map latlong eastlong degrees;
id id;
run;

data alberta_data; set mapsgfk.canada_attr (where=(ID1='CA-48'));
cd_num=.; cd_num=scan(idname,-1,' ');
run;

legend1 label=(position=top 'Division No.') value=(justify=left)
position=(bottom center) shape=bar(.15in,.15in);

proc gmap data=alberta_data map=alberta_map all;
id id;
choro cd_num / discrete legend=legend1;
run;

[levender8622]

Thank you so much for all who share their valuable ideas to this question, all have been a great help. I myself also got help from SAS technical support which I'd like to share here, I hope it can be somewhat helpful to others as well.

Below are quote from Marcia Surratt from SAS support, and thank her again for the help.

Since the data is projected to have the curvature that you see, you can project the data to another projection. For example:

  /* Use the FROM value for Statistics Canada Lambert projection and the TO value for WGS 84 */

proc gproject data=ab out=new from="EPSG:3347" to="EPSG:4326" ;

  id csduid;

proc gmap data=new map=new;

  id csduid;

  choro cdname / statistic=first;

run;

quit;

Since that displays the map wider than the original, you can project the data to arbitrary cartesian to create a better look:

proc gproject data=new out=proj degree eastlong;

  id csduid;

proc gmap data=proj map=proj;

  id csduid;

  choro cdname / statistic=first;

run;

quit;

I also followed up on her as for how to figure out the "from" and "to" values, and below are her responses:

The SASHELP.PROJ4DEF contains the proj4 string for many known projections.  If you do not know the name of the projection, you can include each of the parameters for the proj4 string to include in the FROM or TO option on the PROC GPROJECT statement.

You can find information on the parameters for the PROJ4 string on the following page:

http://proj4.org/parameters.html#parameter-list

By using the values in the lcsd000b16a_e.prj file, you can determine the values to use for the strings and determine if one exists in SASHELP.PROJ4DEF.  Below you can see the values for EPSG:3347 (description NAD83 / Statistics Canada Lambert) from SASHELP.PROJ4DEF and the values from the shapefile .prj file:

/* from SASHELP.PROJ4DEF data set */

/*+proj      Projection name*/

+proj=lcc

/*+lat_1     Latitude of first standard parallel*/

+lat_1=49

/*+lat_2     Latitude of second standard parallel*/

+lat_2=77

/*+lat_0     Latitude of origin*/

+lat_0=63.390675

/*+lon_0     Central meridian*/

+lon_0=-91.86666666666666

/*+x_0       False easting*/

+x_0=6200000

/*+y_0       False northing*/

+y_0=3000000

/*+datum     Datum name*/

+datum=NAD83

/*+units     meters, US survey feet, etc.*/

+units=m

+no_defs

/*From lcsd000b16a_e.prj*/

PROJCS["PCS_Lambert_Conformal_Conic",

GEOGCS["GCS_North_American_1983",

/*+datum     Datum name*/

DATUM["D_North_American_1983",

SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],

UNIT["Degree",0.0174532925199433]],

/*+proj      Projection name*/

PROJECTION["Lambert_Conformal_Conic"],

/*+x_0       False easting*/

PARAMETER["False_Easting",6200000.0],

/*+y_0       False northing*/

PARAMETER["False_Northing",3000000.0],

/*+lon_0     Central meridian*/

PARAMETER["Central_Meridian",-91.86666666666666],

/*+lat_1     Latitude of first standard parallel*/

PARAMETER["Standard_Parallel_1",49.0],

/*+lat_2     Latitude of second standard parallel*/

PARAMETER["Standard_Parallel_2",77.0],

/*+lat_0     Latitude of origin*/

PARAMETER["Latitude_Of_Origin",63.390675],

/*+units     meters, US survey feet, etc.*/

UNIT["Meter",1.0]]

The value of WGS84 (EPSG:4326) is latitude and longitude geographic coordinates.