1

I have the following file.grd with the content shown below. This file is an output of the program CALPOST, which reads the outputs of the program CALPUFF. I can read the file.grd using R and raster:

library(raster)
r1 <- raster("file.grd")
class       : RasterLayer 
dimensions  : 29, 29, 841  (nrow, ncol, ncell)
resolution  : 3, 3  (x, y)
extent      : -43.5, 43.5, -43.5, 43.5  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : /home/sergio/Documents/CALPUFF_EPA/POST/CALPOST/rank(1)_so2_1hr_conc.grd 
names       : rank.1._so2_1hr_conc 
values      : 0, 17.56  (min, max)
spplot(r1, scales = list(draw = T))
origin(r1)
[1] 1.5 1.5

Figure

However, the raster does not have coord. ref. and it is centered at 0,0. It seems that the units ar km.

CALPUFF needs a control file, shown here. CALPOST also needs a control file, shown here.

The original data is Lambert Conformal Conic, with center coordinates at 33.039S and 71.626W. It is gridded data, with DX 3 km, but when I rewrite the origin, nothing changes:

origin(r1) <- c(-71.626, - 33.039)
origin(r1)
[1] 0.374 -0.039

I dont know what I'm missing. Below the data file.grd and also on this link

DSAA
          29          29
     -42.000      42.000
     -42.000      42.000
  0.0000E+00  0.1756E+02
 0.0000E+00  2.3093E-05  1.8688E-04  3.8684E-04  6.0536E-04  9.0230E-04  1.2761E-03  1.7696E-03  2.3521E-03  3.0307E-03  3.8392E-03  4.6814E-03  5.5433E-03  6.4334E-03  7.2733E-03  8.0669E-03  9.1398E-03  9.9640E-03  1.0748E-02  1.1413E-02  1.1886E-02  1.2337E-02  1.6720E-02  2.1224E-02  2.5071E-02  2.7608E-02  2.8188E-02  2.6685E-02  2.3531E-02
 0.0000E+00  3.7950E-05  1.7743E-04  3.7044E-04  6.2073E-04  9.4318E-04  1.3566E-03  1.9102E-03  2.5731E-03  3.3713E-03  4.2576E-03  5.2858E-03  6.3014E-03  7.3816E-03  8.4255E-03  9.6190E-03  1.0808E-02  1.1902E-02  1.2875E-02  1.3681E-02  1.4455E-02  1.8404E-02  2.5050E-02  3.1541E-02  3.6823E-02  3.9920E-02  4.0072E-02  3.7277E-02  3.2298E-02
 0.0000E+00  5.0970E-05  1.5543E-04  3.6608E-04  6.9616E-04  1.0156E-03  1.4173E-03  2.0417E-03  2.7907E-03  3.6841E-03  4.7173E-03  5.8408E-03  7.1526E-03  8.4093E-03  9.7571E-03  1.1212E-02  1.2794E-02  1.4171E-02  1.5365E-02  1.6447E-02  1.9782E-02  2.5817E-02  3.4849E-02  4.3411E-02  5.0012E-02  5.3339E-02  5.2602E-02  4.8120E-02  4.1067E-02
 0.0000E+00  3.8226E-05  1.6111E-04  4.3967E-04  7.2468E-04  1.0982E-03  1.5254E-03  2.1317E-03  2.9607E-03  3.9636E-03  5.1389E-03  6.4682E-03  8.0991E-03  9.6253E-03  1.1280E-02  1.3046E-02  1.5664E-02  1.6765E-02  1.8440E-02  2.5743E-02  3.2626E-02  3.7106E-02  4.4654E-02  5.4959E-02  6.2388E-02  6.5503E-02  6.3573E-02  5.7275E-02  4.9131E-02
 0.0000E+00  2.8361E-05  2.0757E-04  4.8299E-04  7.8377E-04  1.1616E-03  1.6423E-03  2.2647E-03  3.0797E-03  4.1869E-03  5.5333E-03  7.0912E-03  8.8389E-03  1.0915E-02  1.2808E-02  1.5050E-02  1.7175E-02  1.9752E-02  2.8522E-02  3.9933E-02  5.0321E-02  5.6464E-02  5.6459E-02  6.3572E-02  7.1137E-02  7.3564E-02  7.0607E-02  6.4777E-02  5.6086E-02
 3.8920E-06  3.4088E-05  2.3495E-04  4.9641E-04  7.9677E-04  1.2046E-03  1.7365E-03  2.4463E-03  3.3158E-03  4.4443E-03  5.8701E-03  7.6757E-03  9.7502E-03  1.2001E-02  1.4509E-02  1.7135E-02  1.9783E-02  2.6770E-02  4.1523E-02  5.7858E-02  7.1953E-02  7.9370E-02  7.7673E-02  6.8323E-02  7.3773E-02  7.5744E-02  7.4330E-02  7.1109E-02  6.4770E-02
 4.6650E-06  5.3125E-05  2.2971E-04  4.9011E-04  7.9097E-04  1.2220E-03  1.8066E-03  2.5828E-03  3.6019E-03  4.9037E-03  6.4067E-03  8.2342E-03  1.0622E-02  1.3148E-02  1.5883E-02  1.9183E-02  2.2483E-02  3.6494E-02  5.6392E-02  7.7772E-02  9.4868E-02  1.0223E-01  9.8894E-02  9.1976E-02  8.5568E-02  8.2358E-02  7.8415E-02  8.2660E-02  8.1565E-02
 9.0853E-06  4.5735E-05  1.7457E-04  4.4981E-04  7.6381E-04  1.2079E-03  1.8282E-03  2.6508E-03  3.7993E-03  5.2894E-03  7.0853E-03  9.2508E-03  1.1695E-02  1.4364E-02  1.7613E-02  2.3395E-02  3.7178E-02  5.5934E-02  7.0613E-02  9.5801E-02  1.1420E-01  1.2077E-01  1.2270E-01  1.2832E-01  1.4231E-01  1.5629E-01  1.5396E-01  1.2570E-01  1.1744E-01
 1.3449E-05  3.3350E-05  1.3247E-04  3.8047E-04  7.1439E-04  1.1434E-03  1.7988E-03  2.6687E-03  3.9760E-03  5.6715E-03  7.7813E-03  1.0330E-02  1.3286E-02  1.6546E-02  2.0386E-02  3.9753E-02  6.6610E-02  9.6724E-02  1.1874E-01  1.2262E-01  1.2387E-01  1.3454E-01  1.5352E-01  1.9410E-01  2.4461E-01  2.7341E-01  2.5555E-01  1.9458E-01  1.8167E-01
 9.1904E-06  3.4958E-05  1.2374E-04  3.6315E-04  6.3317E-04  1.0605E-03  1.7125E-03  2.6527E-03  4.0605E-03  5.9841E-03  8.3886E-03  1.1404E-02  1.4828E-02  1.8718E-02  3.1919E-02  6.3197E-02  1.0801E-01  1.5739E-01  1.8956E-01  1.8917E-01  1.7163E-01  1.5877E-01  2.0021E-01  2.9638E-01  3.7607E-01  3.9109E-01  3.3307E-01  2.8664E-01  2.6440E-01
 8.9774E-06  3.4562E-05  1.0773E-04  3.2365E-04  5.9564E-04  1.0196E-03  1.7616E-03  2.7062E-03  4.2160E-03  6.1367E-03  8.8644E-03  1.2319E-02  1.6340E-02  2.1849E-02  7.0764E-02  1.3386E-01  1.6927E-01  2.3436E-01  2.7092E-01  2.6892E-01  2.8224E-01  3.1177E-01  2.6870E-01  4.1396E-01  4.7704E-01  4.3365E-01  3.4215E-01  3.7222E-01  3.0938E-01
 6.8659E-06  5.1002E-05  1.2197E-04  3.8943E-04  5.7017E-04  1.3303E-03  2.1220E-03  3.4267E-03  5.1381E-03  7.2195E-03  9.8913E-03  1.4214E-02  2.0024E-02  4.2711E-02  1.3356E-01  3.0441E-01  3.7047E-01  3.2679E-01  3.4590E-01  4.1613E-01  5.5325E-01  5.9729E-01  4.3641E-01  5.2110E-01  4.8346E-01  3.7824E-01  3.6892E-01  3.3164E-01  2.5604E-01
 1.2771E-05  4.1851E-05  1.1775E-04  4.5861E-04  7.4031E-04  1.4358E-03  2.5067E-03  4.0290E-03  5.9426E-03  8.5314E-03  1.1963E-02  1.9169E-02  3.5919E-02  8.4290E-02  4.6568E-01  6.4036E-01  7.0378E-01  5.7180E-01  5.3870E-01  9.0618E-01  1.0517E+00  8.3330E-01  5.9511E-01  5.1895E-01  4.3451E-01  3.9284E-01  3.0318E-01  2.6472E-01  2.1832E-01
 0.0000E+00  2.4555E-05  1.1037E-04  5.6455E-04  8.5457E-04  1.5938E-03  2.9092E-03  4.5201E-03  7.4105E-03  1.1019E-02  1.6595E-02  3.4526E-02  6.1762E-02  1.8492E-01  9.2356E-01  1.2341E+00  1.1151E+00  1.1105E+00  1.2076E+00  1.6925E+00  1.2486E+00  7.2974E-01  5.3236E-01  4.6695E-01  3.7606E-01  3.1632E-01  3.1262E-01  2.7028E-01  2.4246E-01
 0.0000E+00  2.1050E-05  1.4643E-04  6.2914E-04  9.4963E-04  1.7323E-03  3.2432E-03  6.2134E-03  1.0406E-02  1.5140E-02  2.3895E-02  5.0954E-02  1.5081E-01  1.9922E+00  7.9298E-01  4.2464E+00  1.3935E+00  3.7181E+00  2.0623E+00  1.5156E+00  8.9091E-01  7.4117E-01  5.4132E-01  5.5443E-01  4.3044E-01  3.3291E-01  4.5559E-01  5.2177E-01  3.8192E-01
 0.0000E+00  2.7942E-05  1.6173E-04  7.0571E-04  1.1798E-03  1.8767E-03  3.9339E-03  8.3184E-03  1.2510E-02  2.0198E-02  4.3839E-02  9.2987E-02  2.5953E-01  8.6610E+00  6.2197E+00  6.6860E+00  7.8573E+00  3.5901E+00  1.5852E+00  2.1020E+00  2.1212E+00  1.1834E+00  7.1424E-01  6.0548E-01  5.5719E-01  5.2757E-01  5.4074E-01  5.2242E-01  3.5597E-01
 0.0000E+00  3.0847E-05  1.8819E-04  7.4621E-04  1.3718E-03  2.5888E-03  5.4874E-03  9.3994E-03  1.4285E-02  2.7333E-02  5.7756E-02  1.1654E-01  7.6885E-01  3.4297E+00  2.2531E+00  1.0909E+01  4.4473E+00  1.2886E+01  3.8882E+00  2.9077E+00  2.0132E+00  1.1176E+00  9.2729E-01  7.2741E-01  5.7705E-01  4.5938E-01  3.7076E-01  3.0334E-01  2.2049E-01
 0.0000E+00  1.2605E-04  2.5089E-04  9.5202E-04  1.4772E-03  2.7115E-03  5.9853E-03  1.0416E-02  1.6033E-02  3.1423E-02  6.0770E-02  1.2080E-01  8.6492E-01  2.0094E+00  4.3911E+00  1.1424E+01  1.7562E+01  5.3383E+00  2.7597E+00  1.6010E+00  1.7573E+00  1.2155E+00  9.1931E-01  6.7800E-01  5.1261E-01  4.0333E-01  3.6377E-01  3.2836E-01  2.5139E-01
 0.0000E+00  1.6156E-04  3.5877E-04  9.9670E-04  1.5143E-03  3.5621E-03  7.0197E-03  1.0874E-02  1.6392E-02  3.1536E-02  5.6895E-02  1.0836E-01  1.0381E+00  1.2620E+00  1.1916E+00  4.6442E+00  6.8940E+00  8.0981E+00  3.3500E+00  2.1657E+00  1.5787E+00  1.0832E+00  7.2986E-01  5.2636E-01  3.6875E-01  3.3770E-01  3.4180E-01  3.2152E-01  2.6129E-01
 4.7258E-05  1.9043E-04  4.2045E-04  7.8787E-04  1.4634E-03  3.6081E-03  7.3602E-03  1.0785E-02  1.6882E-02  2.8880E-02  5.0819E-02  1.1401E-01  9.9010E-01  8.5282E-01  1.1379E+00  2.6599E+00  4.2333E+00  4.8562E+00  3.6598E+00  1.9153E+00  1.2215E+00  1.2600E+00  1.0174E+00  7.0508E-01  5.3209E-01  5.3925E-01  6.0435E-01  6.0185E-01  4.7513E-01
 5.6205E-05  2.3778E-04  4.8068E-04  8.4275E-04  1.3669E-03  3.6855E-03  6.4851E-03  1.0862E-02  1.5813E-02  2.5756E-02  4.5431E-02  1.1235E-01  8.2222E-01  1.0943E+00  1.2469E+00  2.6507E+00  2.3118E+00  3.3659E+00  3.1313E+00  1.7659E+00  8.1603E-01  6.3353E-01  6.0422E-01  5.2137E-01  3.6190E-01  4.2709E-01  5.9228E-01  6.6879E-01  5.5412E-01
 6.2629E-05  1.9778E-04  5.0494E-04  8.6120E-04  1.3960E-03  3.7352E-03  5.7005E-03  1.0573E-02  1.5221E-02  2.3605E-02  3.7621E-02  1.0593E-01  6.7042E-01  1.1576E+00  9.4279E-01  2.1441E+00  1.5490E+00  2.2060E+00  2.3902E+00  1.7033E+00  9.7458E-01  7.4246E-01  4.8847E-01  4.5047E-01  3.6818E-01  2.5435E-01  3.8572E-01  4.9249E-01  4.3812E-01
 3.2373E-05  2.0489E-04  5.0890E-04  8.3302E-04  1.2870E-03  2.4735E-03  5.4561E-03  9.9152E-03  1.3976E-02  1.9326E-02  3.0635E-02  8.5000E-02  4.3674E-01  7.4839E-01  8.5171E-01  1.4387E+00  1.2432E+00  1.8569E+00  2.2520E+00  1.5156E+00  8.5986E-01  8.9529E-01  6.6317E-01  4.1330E-01  4.1982E-01  3.2200E-01  2.4937E-01  2.7618E-01  2.6629E-01
 4.7129E-05  1.5609E-04  4.3652E-04  7.7467E-04  1.2976E-03  2.0612E-03  5.1317E-03  7.4016E-03  1.2429E-02  1.6777E-02  2.5254E-02  5.8930E-02  2.9484E-01  6.9730E-01  8.5513E-01  1.1794E+00  8.5599E-01  1.3817E+00  1.9059E+00  1.6426E+00  1.0146E+00  5.7987E-01  6.1322E-01  4.5361E-01  3.2063E-01  3.0344E-01  2.1985E-01  2.5718E-01  1.8791E-01
 5.1218E-05  1.2252E-04  3.9714E-04  6.8158E-04  1.0209E-03  2.0131E-03  3.4006E-03  6.5375E-03  1.0651E-02  1.4228E-02  1.8545E-02  4.2591E-02  1.8652E-01  4.8480E-01  7.1422E-01  1.0561E+00  7.2767E-01  1.0217E+00  1.3140E+00  1.3180E+00  9.5974E-01  6.0042E-01  4.5974E-01  4.7367E-01  3.9305E-01  3.1421E-01  2.1875E-01  2.3723E-01  2.1012E-01
 4.0378E-05  1.2599E-04  2.4897E-04  5.2118E-04  8.8089E-04  1.5786E-03  2.7775E-03  5.6129E-03  7.6428E-03  1.1676E-02  1.4996E-02  2.3938E-02  9.2613E-02  2.5509E-01  5.9837E-01  9.2069E-01  7.1944E-01  9.3824E-01  1.2534E+00  1.2538E+00  9.7433E-01  6.5676E-01  4.6241E-01  4.3209E-01  3.8919E-01  3.2911E-01  2.8718E-01  1.8703E-01  1.8383E-01
 4.8513E-05  1.2474E-04  2.2300E-04  3.7650E-04  8.2775E-04  1.4225E-03  2.4443E-03  3.4032E-03  6.2917E-03  8.2329E-03  1.1194E-02  1.4896E-02  3.6522E-02  1.7390E-01  4.9290E-01  7.5172E-01  6.3928E-01  6.8075E-01  1.0800E+00  1.3223E+00  1.0048E+00  7.5802E-01  5.1840E-01  3.9304E-01  3.9334E-01  3.7719E-01  2.7318E-01  2.0878E-01  1.3180E-01
 3.6047E-05  1.0170E-04  2.0916E-04  3.3936E-04  6.4517E-04  1.1292E-03  1.8035E-03  2.5421E-03  3.8429E-03  5.5497E-03  8.1322E-03  1.0197E-02  3.0846E-02  1.4244E-01  3.6332E-01  5.4226E-01  4.8800E-01  3.8822E-01  7.7163E-01  1.0452E+00  9.1709E-01  7.5702E-01  5.0951E-01  3.7063E-01  3.2574E-01  3.4428E-01  3.3371E-01  1.9352E-01  1.2368E-01
 1.9114E-05  8.2234E-05  1.9135E-04  2.9189E-04  4.8302E-04  8.5810E-04  1.3917E-03  2.0693E-03  2.9403E-03  4.3411E-03  5.8867E-03  7.7454E-03  2.5139E-02  9.8444E-02  2.2913E-01  3.3447E-01  3.1285E-01  2.3002E-01  4.6255E-01  6.5796E-01  6.4053E-01  5.9121E-01  4.2217E-01  2.8807E-01  2.6940E-01  2.4935E-01  2.9279E-01  1.8162E-01  1.1357E-01





        1   HOUR AVERAGE  CONCENTRATION VALUES AT EACH RECEPTOR  (ug/m**3)    

         SO2           1

         RANK          1
2

The grid is 3km cells in that LCC coordinate system. So I think the process is as follows:

  • Convert that centre point to the LCC coordinate system.

  • Set the grid extent (using extent(d) <- and not origin).

Working out what the coordinate system is in the standard (PROJ4) format is tricky. Its a Lambert Conformal Conic, that's okay, the parameters are in there, and the projection is based on a spherical earth rather than the standard ellipsoid that GPS uses (The "WGS84" ellispoid).

The relevant bits of the projection info are:

 Map projection
 (PMAP)                     Default: UTM    ! PMAP = LCC  !

 False Easting and Northing (km) at the projection origin
 (Used only if PMAP= TTM, LCC, or LAZA)
 (FEAST)                    Default=0.0     ! FEAST  = 254.763  !
 (FNORTH)                   Default=0.0     ! FNORTH = 6341.324  !

 Latitude and Longitude (decimal degrees) of projection origin
 (Used only if PMAP= TTM, LCC, PS, EM, or LAZA)
 (RLAT0)                    No Default      ! RLAT0 = 33.039S  !
 (RLON0)                    No Default      ! RLON0 = 71.626W  !

     LCC :  RLON0 identifies central (true N/S) meridian of projection
            RLAT0 selected for convenience

 Matching parallel(s) of latitude (decimal degrees) for projection
 (Used only if PMAP= LCC or PS)
 (XLAT1)                    No Default      ! XLAT1 = 33.039S  !
 (XLAT2)                    No Default      ! XLAT2 = 33.039S  !

     LCC :  Projection cone slices through Earth's surface at XLAT1 and XLAT2

 Datum-region for output coordinates
 (DATUM)                    Default: WGS-84    ! DATUM = NWS-84  !
 NWS-84    NWS 6370KM Radius, Sphere

Now if the region is fairly small (and 45km isn't that big on a global scale) you might be able to get away with converting the centre coordinate to any metric coordinate system, and if the distortion is acceptible then job done. Otherwise that info I cut from the control file needs to be converted to PROJ4 - maybe someone has already done this!

The following code makes an attempt, I'll also point out where its wrong:

Read data:

> d = raster("./data.grd")

Centre in lat-long:

> xy = sf::st_point(c(-71.626, - 33.039))
> p = st_sfc(xy)

Now I'm going to mark these as being WGS84 lat-long coordinates - which is wrong, but possibly not by much. Once you've worked out the transform parameters from the description replace this with a string that describes lat-long on that spherical earth:

> st_crs(p)=4326

Now transform that centre to the grid system. Here' I'll transform to UTM in metres for that location, which is wrong, but might be close enough. For precision, work out the LCC parameters and replace the number:

> p2 = st_transform(p, 32619)

Now the upper right and lower left coordinates - multiply by 1000 since the coords are in m but the grid is sized 3km. There are 43.5 km from centre to the corners:

> UR = st_coordinates(p2) + 1000*c(43.5,43.5)
> LL = st_coordinates(p2) - 1000*c(43.5,43.5)

Now set the extent of the raster in projected coordinates:

> extent(d) = c(LL[1], UR[1], LL[2], UR[2])

And tell the raster what the projected coordinate system is:

> projection(d)="+init=epsg:32619"

Save it.

> writeRaster(d,"d2.tiff",overwrite=TRUE)

Load into QGIS and are we even in the right part of the world?

enter image description here

West coast of South America?

With the approximations and incorrect assumptions that might be five metres out (which would be okay I guess) or five kilometres out (which wouldn't be good with 3km cells).....

  • Well, that is very good. A little bit to the north but I will try using the EPSG code this region and will let you know. Thanks – Sergio Dec 16 '18 at 19:02
  • 1
    The proj4 page for LCC is here: proj4.org/operations/projections/lcc.html - but beware of different ways of specifying parameters, the false origin, and having some things in metres and some in km! – Spacedman Dec 16 '18 at 19:07
1

Using QGIS, I get good results with this custom CRS:

+proj=lcc +lat_1=-33.039 +lat_0=-33.039 +lon_0=-71.626 +k_0=1 +x_0=0 +y_0=0 +a=6370000 +b=6370000 +towgs84=0,0,0,0,0,0,0 +units=km +no_defs

getting this output:

enter image description here

The False Easting and Northing seem to have no impact, as the projection center is the same as the image center (the blue dot in my image).

By default, (-42;42) is taken as the origin in the upper left corner (mid-cell). The custom CRS takes it as kilometers from the center.

  • I was looking for a custom CRS like you provided. Thanks! – Sergio Dec 18 '18 at 18:13

Your Answer

By clicking "Post Your Answer", you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.