lat long and GPROJECT and GINSIDE

[Matthijs]

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;

[ballardw]

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?

[Matthijs]

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;

[DaveShea · Accepted solution]

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;

 **************************************************************;
 *Give some of our values retrieved from MAPSGFK.WORLD more    ;
 *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;

[Matthijs]

Thank you so much!