Obsidian | Level 7

## lat long and GPROJECT and GINSIDE

Hi,

I have a set with just latitude and longitude. I want two things:

1. Check if lat long is in the Netherlands

If so

2. In witch provice

So i've tried a couple of things.

https://stackoverflow.com/questions/42077817/map-latitude-and-longitude-to-state-in-sas (and tried to rebuild it to NL)

https://stackoverflow.com/questions/40025725/sas-base-getting-country-using-longitude-and-latitude (and tried to link it to province)

It looks like I don't understand the procedure. From the seond link I do get the country. But the generated x and y do not match the x and y in the mapsgfk.netherlands set.

Who can help me out.

``````/* GPS location in the Netherlands (Noord Holland)*/

data gps;
input lat lon;
datalines;
52.5155 4.99475
;
run;``````
1 ACCEPTED SOLUTION

Accepted Solutions
Lapis Lazuli | Level 10

## Re: lat long and GPROJECT and GINSIDE

Hi Matthijs,

I too am trying to learn how to map Latitude/Longitude points on maps with SAS. You are not alone in finding it confusing

I have been working on a similar bit of code to yours as a way of learning how this works and I am pretty sure that it does what you need it to do; that is, take a pair of Latitude/Longitude values and find out both the Country and the Province/Region that the point belongs to.

The SAS code that I have included below is a bit long-winded, but I hope that the step-by-step approach and the comments help you. I used Google Maps to take a stroll around The Netherlands, starting with your example location, and found a few more places of interest and noted their Latitude/Longitude. Then I got SAS to do the hard work.

Best of luck,

Downunder Dave
New Zealand

``````******************************************************************************;
*Step 1 of 3 - Read our Latitude/Longitude pairs into a SAS dataset.          ;
******************************************************************************;

******************************************************************************;
*Read in our Latitude/Longitude values stored in Degrees.                     ;
******************************************************************************;
Data GPS_Locations;
Input Latitude
Longitude
@20 Location_Name \$35.
;

**************************************************************;
*Copy Latitude into column Y                                  ;
**************************************************************;
Y=Latitude;

**************************************************************;
*Copy Longitude into column X                                 ;
**************************************************************;
X=Longitude;

**************************************************************;
*Create a unique record identifier.                           ;
**************************************************************;
Location_ID=_N_;

Datalines;
52.5155   4.99475  Baanstee-West
51.770285 4.070793 Melissant
51.985550 5.898843 Arnhem Railway Station
52.103446 4.350538 National Car Museum Den Haag
52.963962 4.771389 Dutch Navy Museum
51.495199 3.617950 Middelburg Railway Station
53.079383 6.582999 Federal Government Office Vries
;
Run;

******************************************************************************;
*Step 2 of 3 - Which Country does each of our Latitude/Longitudes live in ?   ;
******************************************************************************;

******************************************************************************;
*Project the X/Y values so that we can get the Country that each point is in. ;
*Keep the original Latitude/Longitude values safe in their own columns.       ;
*The PARMIN and PARMENTRY options tell SAS to look in the MAPSGFK.PROJPARM    ;
*dataset where PROJ_MAP="WORLD" and that row tells Proc GProject how to do    ;
*this particular projection. In this case we will:                            ;
*1)Use MILLER2 projection with options:                                       ;
*2)Use the options DEGREES EASTLONG DUPOK                                     ;
*                                                                             ;
*This configuration could have been manually coded in the Proc GProject       ;
*statement.                                                                   ;
******************************************************************************;
Proc GProject   Data     =WORK.GPS_Locations
Out      =WORK.Projected_For_Country
ParmIn   =MAPSGFK.PROJPARM
ParmEntry=WORLD
;
ID Location_ID;
Run;

******************************************************************************;
*Which Country does each projected X/Y combination sit inside ?               ;
******************************************************************************;
Proc GInside    Data=WORK.Projected_For_Country
Map =MAPSGFK.WORLD
Out =WORK.GPS_Location_Countries
;
**************************************************************;
*I am really not at all sure what this bit does.......sorry ! ;
**************************************************************;
ID ID IDName ISO;
Run;

******************************************************************************;
*Step 3 of 3 - Which Region does each of our Latitude/Longitudes live in ?    ;
*In the next steps we are going to run another pair of GProject and GInside   ;
*PROCs in order to get the more detailed Region that our points are found in. ;
******************************************************************************;

******************************************************************************;
*Protect some columns from the upcoming GProject and GInside calls by renaming;
*the columns we want to protect.                                              ;
******************************************************************************;
Data GPS_Location_Countries;
Set GPS_Location_Countries;

**************************************************************;
*Re-set our X/Y values to be their original Longitude/Latitude;
**************************************************************;
Y = Latitude;
X = Longitude;

**************************************************************;
*useful names.                                                ;
**************************************************************;
Country_Code       = ID;
Country_Name       = IDName;
Country_ISO_Code   = ISO;

Keep   X
Y
Country_Code
Country_Name
Country_ISO_Code
Latitude
Longitude
Location_Name
Location_ID
;
Run;

******************************************************************************;
*Now project our X/Y values again,this time to get the Region each point is in;
*                                                                             ;
*This time the PARMIN and PARMENTRY options tell SAS to look in the           ;
*MAPSGFK.PROJPARM dataset where PROJ_MAP="NETHERLANDS_ALL" and that row tells ;
*Proc GProject how to do this particular projection. In this case we will:    ;
*1)Use MILLER2 projection with options:                                       ;
*2)Use the options DEGREES EASTLONG DUPOK                                     ;
*                                                                             ;
*This configuration could have been manually coded in the Proc GProject       ;
*statement.                                                                   ;
******************************************************************************;
Proc GProject   Data     =WORK.GPS_Location_Countries
Out      =Projected_For_Region
ParmIn   =MAPSGFK.PROJPARM
ParmEntry=NETHERLANDS_ALL
;
ID Location_ID;
Run;

******************************************************************************;
*Which Region of Netherlands does each projected X/Y combination sit inside ? ;
******************************************************************************;
Proc GInside    Data=Projected_For_Region
Map =MAPSGFK.NETHERLANDS_ALL
Out =GPS_Location_Regions
;
ID ID IDName Segment;

Run;

******************************************************************************;
*Now tidy up the columns and column-order in our final dataset to make it read;
*properly, left to right.                                                     ;
******************************************************************************;
Data GPS_Location_Regions;
Informat  Latitude
Longitude
Location_Name
Location_ID
X
Y
Country_Code
Country_Name
Country_ISO_Code
Map_Region_Code
Map_Region_Name
SEGMENT
;

Set GPS_Location_Regions;

Map_Region_Code = ID;
Map_Region_Name = IDName;

Keep    Latitude
Longitude
Location_Name
Location_ID
X
Y
Country_Code
Country_Name
Country_ISO_Code
Map_Region_Code
Map_Region_Name
SEGMENT
;

Run;``````

4 REPLIES 4
Super User

## Re: lat long and GPROJECT and GINSIDE

It would probably help to show the actual code you were using.

When you say "the generated x and y do not match the x and y in the mapsgfk.netherlands set." do you mean that a projected point that you suspect may be in the area does not have a matching (x,y) in the mapsgfk.Netherlands_all data set? Since the coordinates in the map data set would be boundaries of areas then it is very likely that a projected point will not exactly match.

Or did mean that you attempted to map a projected point and it did not appear on the map where expected?

Obsidian | Level 7

## Re: lat long and GPROJECT and GINSIDE

Hi,

Good point. I can't share my actual data, but can provide the used query(s).

``````DATA geo_events2;
SET WORK.GEO_EVENTS;
x = lon;
y = lat;
rownum = _n_;
RUN;

proc gproject data=work.geo_events2 out=projected
parmin=mapsgfk.projparm parmentry=world;
id rownum;
run;

proc ginside data=projected map=mapsgfk.world out=countries;
id id idname segment;
run;``````

So with this code I rename my lat long to x and y and create and eventid (rownumber).

With the next step the 'Projected' set contains an x and y.

With proc ginside I can find the country(s) belonging to de gps-data.

Although I can see the Netherlands now, it's not providing me with the province in Netherlands (duh... It's not in mapsgfk.world)

So I looked for the Netherlands in the same library. There's a sets with lat long and x and y. BUT non of the x and y's in the previeus set is near/close to the province x and y. So another ginside has no point because all the x and y are outside the x and y in the mapsgfk.netherlands set.

Perhaps there's a simpler way of getting my question answered then the above code. What I want is to add country and province to my latlong data in geo_events.

p.s

Also tried the code below... Same results, different library

``````goptions reset=global border;
data gpscounties;
set work.geo_events;
x=long*arcos(-1)/180;
x=x*(-1);
y=lat*arcos(-1)/180;
run;
proc ginside data=gpscounties map=mapssas.world out=gpscounties;
id id cont;
run;``````
Lapis Lazuli | Level 10

## Re: lat long and GPROJECT and GINSIDE

Hi Matthijs,

I too am trying to learn how to map Latitude/Longitude points on maps with SAS. You are not alone in finding it confusing

I have been working on a similar bit of code to yours as a way of learning how this works and I am pretty sure that it does what you need it to do; that is, take a pair of Latitude/Longitude values and find out both the Country and the Province/Region that the point belongs to.

The SAS code that I have included below is a bit long-winded, but I hope that the step-by-step approach and the comments help you. I used Google Maps to take a stroll around The Netherlands, starting with your example location, and found a few more places of interest and noted their Latitude/Longitude. Then I got SAS to do the hard work.

Best of luck,

Downunder Dave
New Zealand

``````******************************************************************************;
*Step 1 of 3 - Read our Latitude/Longitude pairs into a SAS dataset.          ;
******************************************************************************;

******************************************************************************;
*Read in our Latitude/Longitude values stored in Degrees.                     ;
******************************************************************************;
Data GPS_Locations;
Input Latitude
Longitude
@20 Location_Name \$35.
;

**************************************************************;
*Copy Latitude into column Y                                  ;
**************************************************************;
Y=Latitude;

**************************************************************;
*Copy Longitude into column X                                 ;
**************************************************************;
X=Longitude;

**************************************************************;
*Create a unique record identifier.                           ;
**************************************************************;
Location_ID=_N_;

Datalines;
52.5155   4.99475  Baanstee-West
51.770285 4.070793 Melissant
51.985550 5.898843 Arnhem Railway Station
52.103446 4.350538 National Car Museum Den Haag
52.963962 4.771389 Dutch Navy Museum
51.495199 3.617950 Middelburg Railway Station
53.079383 6.582999 Federal Government Office Vries
;
Run;

******************************************************************************;
*Step 2 of 3 - Which Country does each of our Latitude/Longitudes live in ?   ;
******************************************************************************;

******************************************************************************;
*Project the X/Y values so that we can get the Country that each point is in. ;
*Keep the original Latitude/Longitude values safe in their own columns.       ;
*The PARMIN and PARMENTRY options tell SAS to look in the MAPSGFK.PROJPARM    ;
*dataset where PROJ_MAP="WORLD" and that row tells Proc GProject how to do    ;
*this particular projection. In this case we will:                            ;
*1)Use MILLER2 projection with options:                                       ;
*2)Use the options DEGREES EASTLONG DUPOK                                     ;
*                                                                             ;
*This configuration could have been manually coded in the Proc GProject       ;
*statement.                                                                   ;
******************************************************************************;
Proc GProject   Data     =WORK.GPS_Locations
Out      =WORK.Projected_For_Country
ParmIn   =MAPSGFK.PROJPARM
ParmEntry=WORLD
;
ID Location_ID;
Run;

******************************************************************************;
*Which Country does each projected X/Y combination sit inside ?               ;
******************************************************************************;
Proc GInside    Data=WORK.Projected_For_Country
Map =MAPSGFK.WORLD
Out =WORK.GPS_Location_Countries
;
**************************************************************;
*I am really not at all sure what this bit does.......sorry ! ;
**************************************************************;
ID ID IDName ISO;
Run;

******************************************************************************;
*Step 3 of 3 - Which Region does each of our Latitude/Longitudes live in ?    ;
*In the next steps we are going to run another pair of GProject and GInside   ;
*PROCs in order to get the more detailed Region that our points are found in. ;
******************************************************************************;

******************************************************************************;
*Protect some columns from the upcoming GProject and GInside calls by renaming;
*the columns we want to protect.                                              ;
******************************************************************************;
Data GPS_Location_Countries;
Set GPS_Location_Countries;

**************************************************************;
*Re-set our X/Y values to be their original Longitude/Latitude;
**************************************************************;
Y = Latitude;
X = Longitude;

**************************************************************;
*useful names.                                                ;
**************************************************************;
Country_Code       = ID;
Country_Name       = IDName;
Country_ISO_Code   = ISO;

Keep   X
Y
Country_Code
Country_Name
Country_ISO_Code
Latitude
Longitude
Location_Name
Location_ID
;
Run;

******************************************************************************;
*Now project our X/Y values again,this time to get the Region each point is in;
*                                                                             ;
*This time the PARMIN and PARMENTRY options tell SAS to look in the           ;
*MAPSGFK.PROJPARM dataset where PROJ_MAP="NETHERLANDS_ALL" and that row tells ;
*Proc GProject how to do this particular projection. In this case we will:    ;
*1)Use MILLER2 projection with options:                                       ;
*2)Use the options DEGREES EASTLONG DUPOK                                     ;
*                                                                             ;
*This configuration could have been manually coded in the Proc GProject       ;
*statement.                                                                   ;
******************************************************************************;
Proc GProject   Data     =WORK.GPS_Location_Countries
Out      =Projected_For_Region
ParmIn   =MAPSGFK.PROJPARM
ParmEntry=NETHERLANDS_ALL
;
ID Location_ID;
Run;

******************************************************************************;
*Which Region of Netherlands does each projected X/Y combination sit inside ? ;
******************************************************************************;
Proc GInside    Data=Projected_For_Region
Map =MAPSGFK.NETHERLANDS_ALL
Out =GPS_Location_Regions
;
ID ID IDName Segment;

Run;

******************************************************************************;
*Now tidy up the columns and column-order in our final dataset to make it read;
*properly, left to right.                                                     ;
******************************************************************************;
Data GPS_Location_Regions;
Informat  Latitude
Longitude
Location_Name
Location_ID
X
Y
Country_Code
Country_Name
Country_ISO_Code
Map_Region_Code
Map_Region_Name
SEGMENT
;

Set GPS_Location_Regions;

Map_Region_Code = ID;
Map_Region_Name = IDName;

Keep    Latitude
Longitude
Location_Name
Location_ID
X
Y
Country_Code
Country_Name
Country_ISO_Code
Map_Region_Code
Map_Region_Name
SEGMENT
;

Run;``````

Obsidian | Level 7

## Re: lat long and GPROJECT and GINSIDE

Thank you so much!
Discussion stats
• 4 replies
• 2853 views
• 1 like
• 3 in conversation