map creation

[Garik]

Hi. I have confronted to a problem on map creation. I need to add the shops on the map.

the code that creates the map is below. 

So, please give me a hint how to add the coordinates on the map. 

Thanks beforehand. 

%annomac;
%maplabel (mapsgfk.armenia, mapsgfk.armenia_attr, work.labelout, idname, id,
           font=Albany AMT, color=black, size=1.5, hsys=3);
proc gmap map=mapsgfk.armenia data=mapsgfk.armenia;
	*where id eq "AM-01";
   id  id;
   choro id / nolegend annotate=labelout coutline=gray stat=first;
run;
quit;

proc gproject data=myxlfl.sc_coord_idd

 out = sc_coord_id_XY;
 id Segment;
run; 

I have tried to project the coordinates in such way but it was wrong.

 The coordinates of the shops is here. 

ID	Zone	X	Y
AM-01	Yerevan	40.19741295	44.5678061
AM-10	Vayots Dzor	39.76049643	45.33073236
AM-10	Vayots Dzor	39.68937627	45.47180898
AM-10	Vayots Dzor	39.83534526	45.66534589
AM-11	Tavush	41.22988562	44.83745333
AM-01	Yerevan	40.1416615	44.52101183

[ChrisNZ]

You could just get the closest X and Y for your LAT on LONG.

data SHOPS;
 input ROW ID $ LAT LONG;
cards;
1 AM-01 40.19741295 44.5678061
2 AM-10 39.76049643 45.33073236
3 AM-10 39.68937627 45.47180898
4 AM-10 39.83534526 45.66534589
5 AM-11 41.22988562 44.83745333
6 AM-01 40.1416615 44.52101183
run;

proc sql; 
  create table ANNO as
  select SHOPS.ROW, SHOPS.ID, la.Y ,lo.X  
         , 'label' as FUNCTION, put(ROW,z2.) as TEXT, '2' as XSYS, '2' as YSYS
  from SHOPS  
      ,MAPSGFK.ARMENIA la   
      ,MAPSGFK.ARMENIA lo   
  where SHOPS.ID=la.ID 
    and   SHOPS.ID=lo.ID 
  group by 1,2 
  having abs(SHOPS.LAT-la.LAT)  =min(abs(SHOPS.LAT -la.LAT)) 
     and abs(SHOPS.LONG-lo.LONG)=min(abs(SHOPS.LONG-lo.LONG)) ;
quit;

proc gmap  map=MAPSGFK.ARMENIA data=MAPSGFK.ARMENIA; 
   id  id;
   choro id / nolegend annotate=ANNO coutline=gray stat=first;
run;
quit;

[ChrisNZ]

Moved to "SAS/GRAPH" board

[Garik]

Hi. many thanks for help.  one more question.

The code you have written, creates "ANNO" data set, in which X and Y ccordinates exist. So, please explain, what is the format of the "X" and  "Y" . 

I tried to use proc gproject, but in vain.

1

1

AM-01

0.0022782902

-0.006334428

label

01

2

2

2

2

AM-10

-0.005401111

0.0039082507

label

02

2

2

3

3

AM-10

-0.006660053

0.0058276305

label

03

2

2

4

4

AM-10

-0.004002435

0.0083556018

label

04

2

2

5

5

AM-11

0.0202396895

-0.002652059

label

05

2

2

6

6

AM-01

0.0012979944

-0.006904374

label

06

2

2

[ChrisNZ]

Not too sure, they are the coordinates used by the SAS maps, but as long as you have a way to map X and Y to LAT and LON as I did, I doesn't really matter for what you want to do.

[Garik]

So, to be more clear. I have a bunch of coordinates of shops, as I wrote before. I want to understand how  you converted the coordinates into the coordinates that SAS uses (X, Y ). are there another methods besides the way you used for converting. the way of converting the coordinates that you used is not familiar to me.

I know several formats of coordiantes such as decimal, cartesian, degree, polar which are explaind in sas documents. but the said is not known. 

could I use 

proc gprogect for converting the coordinates.

Thanks.

[Reeza]

@Garik I think you're overthinking this. You need the latitude and longitude, X=Longitude and Y=Latitude

If you're source data is projected onto a different coordinate system then you need some modifications or GPROJECT. But try it with just your latitude and longitude first.

[ChrisNZ]

My method is really simple: find the same latitude in the map file as what you have, and get the matching Y.

A slightly more sophisticated method, since for most countries the projection is essentially flat: you can probably create a first degree polynomial equation to transform LAT/LON into Y/X just by looking the four columns in the maps file.

If you want to project, see the post by @ballardw for a third method to match LAT/LON to X/Y.

[ballardw]

You projection is incomplete. You want to reference the map data set to have the "bounds" of the map area and some other options.

I don't know how/what you pasted for the Proc Gproject code but I can't copy and paste to edit to something close to what you need.

From the online documentation example for projecting an annotate data set:

proc gproject data=cities
              degrees
              latlon
              out=citiesp
              parmentry=us48p
              dupok;
id;
run;

Data=cities would be the equivalent of your coordinate set

Degrees says that the data for latitude and longitude are measured in degrees not radians

Latlon says to use the Lat and Long variables from the MAP data set for projection

Out= output data set

Parmentry= the map data set with X, Y and Lat Long variables

Dupok and option not to remove duplicate projections.

Id does not have a variable because you want to project points (cities) not polygons that require identification of area units.

[GraphGuy]

One problem - I think the x & y in your shops data is reversed (which is an easy mistake to make when the lat & long values are very close to being the same, like they are in Armenia).

When working with maps in SAS, you have latitude and longitude coordinates, and then once you 'project' the coordinates you have y and x values. If you're working with both a map, and point data (such as your shop locations), then you have to project the map and the points using exactly the same projection parameters. One way to do that is by using the Proc Gproject's parmout= and parmin= options.

Here is some SAS code that does what I think you are wanting ...

data my_armenia; set mapsgfk.armenia;
run;

proc gproject data=my_armenia out=my_armenia latlong eastlong degrees
parmout=work.projparm;
id id;
run;

data my_shops;
input ID $ 1-5 Zone $ 6-17 lat long;
datalines;
AM-01 Yerevan 40.19741295 44.5678061
AM-10 Vayots Dzor 39.76049643 45.33073236
AM-10 Vayots Dzor 39.68937627 45.47180898
AM-10 Vayots Dzor 39.83534526 45.66534589
AM-11 Tavush 41.22988562 44.83745333
AM-01 Yerevan 40.1416615 44.52101183
;
run;

proc gproject data=my_shops out=my_shops latlong eastlong degrees
parmin=work.projparm parmentry=my_armenia;
id;
run;

data my_shops; set my_shops;
length function $8 style $35 color $12;
xsys='2'; ysys='2'; hsys='3'; when='a';
function='pie'; rotate=360; size=1.5;
style='psolid'; color='cyan'; output;
style='pempty'; color='black'; output;
run;

legend1 label=none position=(bottom left) across=3 mode=share offset=(5,5);

title1 ls=1.5 h=18pt "Shops in Armenia";
proc gmap data=mapsgfk.armenia_attr map=my_armenia all anno=my_shops;
id id;
choro idname / legend=legend1 coutline=grayaa;
run;