BookmarkSubscribeRSS Feed
Jamie_H
Fluorite | Level 6

Hi, 

 

I was wondering if anyone has tried to use the %CODEPOINT2GEOCODE SAS macro (to convert post codes in Great Britain to a usable format) in a locked down Enterprise Guide environment?

 

We dont have the ability to Pipe files and so the code falls over.  

 

I am able to upload all the appropriate postcode .CSV files to the server using the upload task (location = %let myWorkLocation = %sysfunc(pathname(work));  )  

 

Just curious to see if SAS or anyone else has created an amended version of the macro that would allow read-ins - read-outs to this location using methods that are allowed in EG?

 

Thanks

3 REPLIES 3
Jamie_H
Fluorite | Level 6

Just bumping this up to see if anyone know the answer?

Patrick
Opal | Level 21

@Jamie_H

This code appears to be a utility macro which you are allowed to modify as you wish.

Looking into the actual code the one bit which won't work with option NOXCMD is the pipe using DIR for the directory listing.

filename filelist pipe %unquote(%str(%'dir "&PATHIN\*" /b%'));

data _null_;
  infile filelist truncover end=finalObs;
  length fname $50;
  retain i 1;  /* Have to initialize counter    */
  input fname; /* Post code file name from pipe */

  /*--- Import only CSV files and ignore any others. */
  if index(fname, '.csv ') then do;
    /*--- Put name in numbered FNAMEi macro
          var and increment file counter. */
    call symputx('FNAME' || trim(left(put(i,8.))), trim(fname));
    i+1;
  end;
  /*--- If last file, put final count into FILETOTAL macro variable. */
  if finalObs then
    call symput('FILETOTAL', trim(left(put(i-1,8.))));
run;

You can also get such a directory listing using SAS code only.

To make things work for you change the above code bit in the macro with below code which should return the exactly same result.

data _null_;
  length fname $50;
  retain cnt 0;  /* Have to initialize counter    */

  /* Assigns a fileref to the directory and opens the directory */
  _rc=filename("filrf","&PATHIN");
  _did=dopen("filrf");

  /* Make sure directory can be opened */
  if _did eq 0 then
    do;
      put "Directory &pathin cannot be opened or does not exist";
      put "aborting job ... aborting job ... aborting job...";
      abort abend;
      stop;
    end;

  /* Loops through entire directory */
  do _i = 1 to dnum(_did);
    /* Retrieve name of each file */
    fname=dread(_did,_i);

    /*--- Import only CSV files and ignore any others. */
    if index(fname, '.csv ') then
      do;
        /*--- Put name in numbered FNAMEi macro
                 var and increment file counter. */
        cnt+1;
        call symputx('FNAME' || put(cnt,8. -l), trim(fname));
        output;
      end;
  end;

  /*--- If last file, put final count into FILETOTAL macro variable. */
  call symputx('FILETOTAL', put(cnt,8. -l));

  /* Closes the directory and clear the fileref */
  _rc=dclose(_did);
  _rc=filename("filrf");
run;

 

Tom
Super User Tom
Super User

So I pulled that macro from SAS's website. 

https://support.sas.com/downloads/download.htm?did=104271#

 

Looking quickly at it seems to be working way too hard.

 

So I removed all of the macro logic and just created a simple program that does the same thing with simple data steps.

 

I did leave two macro variable assignments at the top so that you can specify the directory where the source files are located and the name of the SAS dataset you want to create.

%let path=C:\Downloads\codepo_gb\Data\CSV\;
%let out=geo;

/*----------------------------------------------------------------------------
Code-Point Open postcode coordinates are in British National Grid
Eastings/Northings in the Ordnance Survey of Great Britain 1936 (OSGB36)
geodetic datum. They need to be changed to longitude/latitude in the World
Geodetic System 1984 (WGS84) datum.

This requires multiple conversions and a datum transformation:

1) Convert National Grid Easting/Northing coordinates to longitude/latitude
   within the OSGB36 geodetic datum.
2) Convert OSGB36 longitude/latitude/ellipsoid height to XYZ cartesian
   coordinates in the OSGB36 datum.
3) Transform the OSGB36 datum XYZ cartesian coordinates to XYZ cartesian
   coordinates in the WGS84 datum.
4) Convert the WGS84 XYZ cartesian coordinates to longitude/latitude in the
   WGS84 geodetic datum.
----------------------------------------------------------------------------*/

/* Read in the files */
data raw ;
  infile "&path.*.csv" dlm='2C0D'x dsd missover lrecl=10000;

/*--- If format of postcode file changes, this will have to change. */
  length PC $7 CY RH LH CC DC WC $9;
  input PC PQ EA NO CY RH LH CC DC WC;

/*--- See "Code-Point Open User guide and technical specification"
                  Ordnance Survey publication for variable details. */
  label PC = 'Postcode'
        PQ = 'Positional quality indicator'
        EA = 'Ordnance Survey grid easting'
        NO = 'Ordnance Survey grid northing'
        CY = 'Country code'
        RH = 'NHS Regional HS code'
        LH = 'NHS HA code'
        CC = 'Admin county code'
        DC = 'Admin district code'
        WC = 'Admin ward code'
  ;
/*--- Remove blanks and upcase postcode. */
  pc = upcase(compress(pc));
run;

/*-------------------------------------------------------------------------*
 |  STEP 1: Convert British National Grid Easting/Northing to
 |          longitude/latitude in OSGB36. See Appendix C.2 in "A Guide
 |          to Coordinate Systems in Great Britain" for methodology.
 *-------------------------------------------------------------------------*/
data step1 (drop=E0 N0 F0 a b Phi0 Lambda0 e2 Nu Rho n n2
                                 n3 PhiPrime DeltaPhi SumPhi M1 M2 M3 M4 M
                                 sinPhiPrime sinPhiPrime2 Eta2 Nu3 Nu5 Nu7
                                 tanPhiPrime tanPhiPrime2 tanPhiPrime4
                                 tanPhiPrime6 secPhiPrime VII VIII IX X XI
                                 XII XIIA deltaE)
;
    set raw (rename=(EA=Easting NO=Northing));

    /*--- Declare OSGB36 ellipsoid and grid constants from Appendix A.       */
    retain E0       400000      /* Easting of National Grid true origin (m)  */
           N0      -100000      /* Northing of Natioanl Grid true origin (m) */
           F0      0.9996012717 /* Scale factor on central meridian          */
           a       6377563.396  /* Airy 1830 ellipsoid semi-major axis (m)   */
           b       6356256.909  /* Airy 1830 ellipsoid semi-minor axis (m)   */
           Phi0    49           /* OSGB36 latitude of true origin (deg)      */
           Lambda0 -2           /* OSGB36 longitude of true origin (deg)     */
           e2                   /* Airy 1830 ellipsoid eccentricity squared  */
           Nu Rho               /* Conversion values (See Appendix C.1)      */
           n n2 n3              /* Constants for Equation C3                 */
    ;

    label Easting    = 'British National Grid Easting (meters)'
          Northing   = 'British National Grid Northing (meters)'
          Lat_OSGB36 = 'OSGB36 Latitude (radians)'
          Lon_OSGB36 = 'OSGB36 Longitude (radians)'
    ;

    /*--- First pass, convert angles to radians and compute constants. */
    if _n_ = 1 then do;
      Phi0    = Phi0    * atan(1) / 45;
      Lambda0 = Lambda0 * atan(1) / 45;
      e2      = (a**2 - b**2) / a**2;    /* Equation B1 */
      n       = (a - b) / (a + b);       /* Equation C1 */
      n2      = n**2;
      n3      = n**3;
    end;

    /*--- Reinitialize for each easting/northing to convert. */
    M        = 0;
    PhiPrime = Phi0;

    /*--- Converting grid easting/northing into OSGB36 longitude/latitude
          is an iterative process. Compute new latitude (PhiPrime) values
          until change is within the required tolerance (0.01). */
    do until ((Northing - N0 - M) < 0.01);
      PhiPrime = ((Northing - N0 - M) / a / F0) + PhiPrime; /* Equation C6 */
      DeltaPhi = PhiPrime - Phi0;
      SumPhi   = PhiPrime + Phi0;

      /*--- Equation C3 is broken into four parts below to compute M. */
      M1 = (1 + n + (1.25 * n2) + (1.25 * n3)) * DeltaPhi;
      M2 = (3 * n + 3 * n2 + 2.625 * n3) * sin(DeltaPhi) * cos(SumPhi);
      M3 = (1.875 * n2 + 1.875 * n3) * sin(2 * DeltaPhi) * cos(2 * SumPhi);
      M4 = 35 / 24 * n3 * sin(3 * DeltaPhi) * cos(3 * SumPhi);
      M  = b * F0 * (M1 - M2 + M3 - M4);
    end;

    /*--- Tolerance is within required specs. */
    PhiPrime = ((Northing - N0 - M) / a / F0) + PhiPrime; /* Equation C7 */

    /*--- Compute Nu, Rho and Eta2 (Eta-squared) from equations in Appendix C.1. */
    sinPhiPrime  = sin(PhiPrime);
    sinPhiPrime2 = sinPhiPrime ** 2;
    Nu           = a * F0 * ((1 - (e2 * sinPhiPrime2)) ** (-0.5));
    Rho          = a * F0 * (1 - e2) * ((1 - (e2 * sinPhiPrime2)) ** (-1.5));
    Eta2         = (Nu / Rho) - 1;  /* Equation C2 */

    /*--- Compute values used repeatedly in equations below. */
    Nu3          = Nu ** 3;
    Nu5          = Nu ** 5;
    Nu7          = Nu ** 7;
    tanPhiPrime  = tan(PhiPrime);
    tanPhiPrime2 = tanPhiPrime ** 2;
    tanPhiPrime4 = tanPhiPrime ** 4;
    tanPhiPrime6 = tanPhiPrime ** 6;
    secPhiPrime  = 1 / cos(PhiPrime);

    /*--- See Appendix C.2 for equations. */
    VII  = tanPhiPrime / 2 / Rho / Nu;
    VIII = tanPhiPrime / 24 / Rho / Nu3 *
           (5 + (3 * tanPhiPrime2) + Eta2 - (9 * tanPhiPrime2 * Eta2))
    ;
    IX   = tanPhiPrime / 720 / Rho / Nu5 *
           (61 + (90 * tanPhiPrime2) + (45 * tanPhiPrime4))
    ;
    X    = secPhiPrime / Nu;
    XI   = secPhiPrime / 6 / Nu3 * ((Nu / Rho) + (2 * tanPhiPrime2));
    XII  = secPhiPrime / 120 / Nu5 *
           (5 + (28 * tanPhiPrime2) + (24 * tanPhiPrime4))
    ;
    XIIA = secPhiPrime / 5040 / Nu7 *
           (61 + (662 * tanPhiPrime2) + (1320 * tanPhiPrime4) +
           (720 * tanPhiPrime6))
    ;

    /*--- Compute latitude and longtitude in OSGB36 geodetic datum. */
    deltaE     = Easting - E0;
    /* Equation C8 */
    lat_OSGB36 = PhiPrime - (VII * deltaE**2) + (VIII * deltaE**4) -
                (IX * deltaE**6)
    ;
    /* Equation C9 */
    lon_OSGB36 = Lambda0 + (X * deltaE) - (XI * deltaE**3) +
                (XII * deltaE**5) - (XIIA * deltaE**7)
    ;
run;

/*-----------------------------------------------------------------------*
 | STEP 2: Convert British National Grid longitude/latitude/height to
 |         XYZ cartesian coordinates within the OSGB36 datum.
 |         See Appendix B.1 in "A Guide to Coordinate Systems in
 |         Great Britain" for methodology
 *-----------------------------------------------------------------------*/
data step2 (drop=a b H e2 Nu);
    set step1;

    /*--- Specify conversion constants.                               */
    retain a  6377563.396 /* Airy 1830 ellipsoid semi-major axis (m)  */
           b  6356256.909 /* Airy 1830 ellipsoid semi-minor axis (m)  */
           H  0           /* Ellipsoid height (m): use zero           */
           e2             /* Airy 1830 ellipsoid eccentricity squared */
    ;

    label X_OSGB36 = 'OSGB36 X coordinate (meters)'
          Y_OSGB36 = 'OSGB36 Y coordinate (meters)'
          Z_OSGB36 = 'OSGB36 Z coordinate (meters)'
    ;

    /*--- On first pass calculate eccentricity squared. */
    if _n_ = 1 then
      e2 = (a**2 - b**2) / a**2  /* Equation B1 */
    ;

    /*--- Convert OSGB36 latitude/longitude to OSGB36 XYZ coordinates. */
    Nu       = a / (sqrt(1 - (e2 * sin(lat_OSGB36)**2)));    /* Equation B2 */
    X_OSGB36 = (Nu + H) * cos(lat_OSGB36) * cos(lon_OSGB36); /* Equation B3 */
    Y_OSGB36 = (Nu + H) * cos(lat_OSGB36) * sin(lon_OSGB36); /* Equation B4 */
    Z_OSGB36 = ((1 - e2) * Nu + H) * sin(lat_OSGB36);        /* Equation B5 */
run;


/*-----------------------------------------------------------------------*
 | Step 3: Transform cartesian XYZ coordinates from Ordnance Survey of
 |         Great Britain 1936 datum (OSGB36) to XYZ coordinates in
 |         World Geodetic System 1984 datum (WGS84). See Section 6.2 in
 |         "A Guide to Coordinate Systems in Great Britain".
 *-----------------------------------------------------------------------*/
data step3 (drop=Tx Ty Tz Rx Ry Rz s sFactor sec2rad);
  set step2;

  /*--- Specify Helmert constants for OSGB36-to-WGS84 transformation.
        If transforming in opposite direction, signs would be reversed. */
  retain Tx   446.448  /* Translation along X-axis (m) */
         Ty  -125.157  /* Translation along Y-axis (m) */
         Tz   542.060  /* Translation along Z-axis (m) */
         Rx   0.1502   /* Rotation about X-axis (sec)  */
         Ry   0.2470   /* Rotation about Y-axis (sec)  */
         Rz   0.8421   /* Rotation about Z-axis (sec)  */
         s   -20.4894  /* Scale change (ppm)           */
         sFactor       /* Scale factor (unitless)      */
  ;

  label X_WGS84 = 'WGS84 X coordinate (meters)'
        Y_WGS84 = 'WGS84 Y coordinate (meters)'
        Z_WGS84 = 'WGS84 Z coordinate (meters)'
  ;

  /*--- On first pass convert rotations to radians and compute scale factor. */
  if _n_ = 1 then do;
    sec2rad  = atan(1) / 3600 / 45;
    Rx       = Rx * sec2rad;
    Ry       = Ry * sec2rad;
    Rz       = Rz * sec2rad;
    sFactor  = 1 + (s / 1000000);
  end;

  /*--- Transform XYZ coordinates from OSGB36 to WGS84 datum using
        Equation 3 in "A Guide to Coordinate Systems in Great Britain."
          -         -   -    -   -               -   -          -
          | X_WGS84 |   | Tx |   | 1+s  -Rz   Ry |   | X_OSGB36 |
          | Y_WGS84 | = | Ty | + |  Rz  1+s  -Rx | . | Y_OSGB36 |
          | Z_WGS84 |   | Tz |   | -Ry   Rx  1+s |   | Z_OSGB36 |
          -         -   -    -   -               -   -          -   */
  X_WGS84 = Tx + (X_OSGB36 * sFactor) - (Y_OSGB36 * Rz)      + (Z_OSGB36 * Ry);
  Y_WGS84 = Ty + (X_OSGB36 * Rz)      + (Y_OSGB36 * sFactor) - (Z_OSGB36 * Rx);
  Z_WGS84 = Tz - (X_OSGB36 * Ry)      + (Y_OSGB36 * Rx)      + (Z_OSGB36 * sFactor);
run;

/*-----------------------------------------------------------------------*
 | Step 4: Convert cartesian XYZ coordinates to longitude/latitude in
 |         World Geodetic System 1984 datum (WGS84). See Appendix B.2
 |         in "A Guide to Coordinate Systems in Great Britain".
 |         A simple index using the postcode variable (pc) is added.
 *-----------------------------------------------------------------------*/
data &out (label="Postcodes from Ordnance Survey Code-Point Open (&SYSDATE9)"
             index=(pc /unique)
             drop=a b e2 Lambda Phi Rho PhiPrime DeltaPhi
                  Nu sinPhi rad2deg );
    set step3;
    retain a  6378137.000  /* WGS84 ellipsoid semi-major axis (m)   */
           b  6356752.3141 /* WGS84 ellipsoid semi-minor axis (m)   */
           e2              /* WGS84 ellipsoid eccentricity squared  */
           rad2deg         /* Radian-to-degree conversion factor    */
    ;

    label Y     = 'WGS84 Latitude (DD)'
          X     = 'WGS84 Longitude (DD)'
          Y_DMS = 'WGS84 Latitude (DMS)'
          X_DMS = 'WGS84 Longitude (DMS)'
          H     = 'WGS84 Ellipsoid Height (meters)'
    ;
    length Y_DMS X_DMS $ 15;

    /*--- On first pass calculate eccentricity squared. */
    if _n_ = 1 then do;
      e2      = (a**2 - b**2) / a**2;          /* Equation B1 */
      rad2deg = 45 / atan(1);
    end;

    /*--- Compute longitude and initial value of latitude. */
    Lambda = atan(Y_WGS84 / X_WGS84);          /* Equation B6 */
    Rho    = sqrt(X_WGS84**2 + Y_WGS84**2);
    Phi    = atan(Z_WGS84 / (Rho * (1 - e2))); /* Equation B7 */

    /*--- Iteratively adjust latitude until change is within tolerance. */
    do until (DeltaPhi < 0.00000001);
      sinPhi   = sin(Phi);
      Nu       = a / sqrt(1 - (e2 * sinPhi**2));             /* Equation B2 */
      PhiPrime = atan((Z_WGS84 + (e2 * Nu * sinPhi)) / Rho); /* Equation B8 */
      DeltaPhi = abs(Phi - PhiPrime);
      Phi      = PhiPrime;
    end;

    /*--- Final latitude change (DeltaPhi) is within desired tolerance.
          Compute ellipsoid height at the computed lat/lon position. */
    H = (Rho / cos(Phi)) - Nu;

    /*--- Convert lat/lon from radians to decimal degrees. */
    Y = Phi    * rad2deg;
    X = Lambda * rad2deg;

    /*--- Also display final location in degrees-minutes-seconds. */

    Y_DMS = put(int(Y), z2.) || 'B0'x || ' '
        || put(int((abs(Y)-abs(int(Y)))*60), z2.)||"' "
        || put(((abs(Y)-abs(int(Y)))*60-int((abs(Y)-abs(int(Y)))*60))*60,z5.2)||'"'
    ;
    X_DMS = put(int(X), z2.) || 'B0'x || ' '
        || put(int((abs(X)-abs(int(X)))*60), z2.)||"' "
        || put(((abs(X)-abs(int(X)))*60-int((abs(X)-abs(int(X)))*60))*60,z5.2)||'"'
    ;
run;

Although I am not if there is a way to replace the 'B0'x hexcode with an encoding setting independent way to represent the degree symbol in those last two statements.  Also I am not sure that the Z2. format is a wide enough to show the DEGREE values, especially if some of them are negative.  But perhaps since it is for UK postcodes it is enough?

SAS Innovate 2025: Save the Date

 SAS Innovate 2025 is scheduled for May 6-9 in Orlando, FL. Sign up to be first to learn about the agenda and registration!

Save the date!

SAS Enterprise Guide vs. SAS Studio

What’s the difference between SAS Enterprise Guide and SAS Studio? How are they similar? Just ask SAS’ Danny Modlin.

Find more tutorials on the SAS Users YouTube channel.

SAS Training: Just a Click Away

 Ready to level-up your skills? Choose your own adventure.

Browse our catalog!

Discussion stats
  • 3 replies
  • 1805 views
  • 0 likes
  • 3 in conversation