turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

Find a Community

- Home
- /
- SAS Programming
- /
- SAS/GRAPH and ODS Graphics
- /
- Re-project Lambert projection from imported shapef...

Topic Options

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

03-27-2017 11:48 AM

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!

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to levender8622

03-27-2017 12:37 PM

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?

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to ballardw

03-27-2017 01:05 PM

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;

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to levender8622

03-27-2017 01:59 PM

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.

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to levender8622

03-30-2017 11:54 AM

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;

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to levender8622

04-03-2017 12:08 PM

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.*