search

cf-convention / cf-conventions

AsciiDoc Source
http://cfconventions.org/cf-conventions/cf-conventions
Creative Commons Zero v1.0 Universal
87 stars 45 forks source link

EASE Grid North 25 km original, in CF grid_mapping #278

Closed twhiteaker closed 4 years ago

twhiteaker commented 4 years ago

As recommended in https://github.com/cf-convention/cf-conventions/issues/222#issuecomment-569789956, I'm inquiring as to whether EASE Grid North 25 km original (EPSG 3408) can be represented using grid_mapping attributes. My team generates hydrologic simulation model data for Arctic regions using that projection. I'd like to store data in CF netCDF, but I don't know if I'm doing it correctly. Like Jim says, there's a good chance I botched it. In particular, there appears to be a number in the WKT related to the degree unit (0.0174532925199433), and I don't know if or how that should be represented in grid_mapping. Here's my attempt at a grid_mapping variable, which includes a crs_wkt attribute. How does it look?

int nsidc_ease_grid_north_25km;
  nsidc_ease_grid_north_25km:grid_mapping_name = "NSIDC_EASE_Grid_North";
  nsidc_ease_grid_north_25km:longitude_of_projection_origin = 0.0;
  nsidc_ease_grid_north_25km:latitude_of_projection_origin = 90.0;
  nsidc_ease_grid_north_25km:false_easting = 0.0;
  nsidc_ease_grid_north_25km:false_northing = 0.0;
  nsidc_ease_grid_north_25km:geographic_crs_name = "GCS_Sphere_International_1924_Authalic";
  nsidc_ease_grid_north_25km:horizontal_datum_name = "Not_specified_based_on_International_1924_Authalic_Sphere";
  nsidc_ease_grid_north_25km:reference_ellipsoid_name = "International 1924 Authalic Sphere";
  nsidc_ease_grid_north_25km:semi_major_axis = 6371228.0;
  nsidc_ease_grid_north_25km:inverse_flattening = 0.0;
  nsidc_ease_grid_north_25km:prime_meridian_name = "Greenwich";
  nsidc_ease_grid_north_25km:longitude_of_prime_meridian = 0.0;
  nsidc_ease_grid_north_25km:projected_crs_name = "Lambert_Azimuthal_Equal_Area";
  nsidc_ease_grid_north_25km:longitude_of_central_meridian  = 0.0;
  nsidc_ease_grid_north_25km:crs_wkt = "PROJCS[\"NSIDC EASE-Grid North\",GEOGCS[\"Unspecified datum based upon the International 1924 Authalic Sphere\",DATUM[\"Not_specified_based_on_International_1924_Authalic_Sphere\",SPHEROID[\"International 1924 Authalic Sphere\",6371228,0,AUTHORITY[\"EPSG\",\"7057\"]],AUTHORITY[\"EPSG\",\"6053\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],   AUTHORITY[\"EPSG\",\"4053\"]],PROJECTION[\"Lambert_Azimuthal_Equal_Area\"],   PARAMETER[\"latitude_of_center\",90],PARAMETER[\"longitude_of_center\",0],   PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,   AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"X\",EAST],AXIS[\"Y\",NORTH],AUTHORITY[\"EPSG\",\"3408\"]]" ;
ethanrd commented 4 years ago

Hi @twhiteaker - Can you include the CDL for the data variable(s) that reference this grid mapping along with the coordinate variables (and auxiliary coordinate variables) used by the data variable(s)? I'm mainly wondering how the units of the coordinate variables relate to the units in the WKT. (I think the 0.01745... is the conversion factor from radians to degrees.)

JonathanGregory commented 4 years ago

Dear Tim @twhiteaker I think the grid_mapping_name should be lambert_azimuthal_equal_area, which is listed in CF Appendix F. The grid_mapping_name is a controlled vocabulary, whereas the projected_crs_name can be anything you choose. Maybe that's where you could record "NSIDC_EASE_Grid_North"? Jonathan

marqh commented 4 years ago

Hello @twhiteaker

regarding:

In particular, there appears to be a number in the WKT related to the degree unit (0.0174532925199433), and I don't know if or how that should be represented in grid_mapping.

The handling of unit of measure within WKT and CF is quite different. There is no expectation in CF that the unit of measure conversion factor to the SI unit is provided explicitly. WKT does provide the conversion from the stated degrees to radians, the SI unit.

So, it is fine to omit this information from the grid mapping attribute metadata, it is implied by other CF metadata.

all the best marqh

marqh commented 4 years ago

Hello @twhiteaker

I am wary of the reference you have provided for this CRS definition. https://epsg.io/3408

The WKT provided appears not conformant to the CRS WKT standard. There are syntax issues suggesting this may have been manufactured a long time ago. I'm not convinced it is valid syntax.

Searching for EASE suggests that http://www.epsg-registry.org/export.htm?wkt=urn:ogc:def:crs:EPSG::6931 may be the definition in this case, but I do not know enough about the EASE grid definitions to understand if this is the case.

The definition in the registry suggests that this is an updated version of 3408:

WGS 84 / NSIDC EASE-Grid 2.0 North | EPSG::6931 | ProjectedCRS | Valid | Northern hemisphere. | Supersedes NSIDC EASE-Grid North, CRS code 3408. This ellipsoidal development is to be preferred for new data.

If you examine the syntax, you will see key differences, such as

I am concerned that the use of https://epsg.io is causing you problems in this definition.

all the best marqh

mjbrodzik commented 4 years ago

Hello, I may be able to help with the issues you are having. You may find this paper helpful, describing issues that we faced in putting EASE-Grid 2.0 data into netCDF-CF compliant files:

Brodzik, M. J., D. G. Long and M. A. Hardman. 2018. Best Practices in Crafting the Calibrated, Enhanced-Resolution Passive-Microwave EASE-Grid 2.0 Brightness Temperature Earth System Data Record. Remote Sensing, 10(11), doi:10.3390/rs10111793.

http://www.mdpi.com/2072-4292/10/11/1793

Some of the things we did will definitely apply to your data as well.

That said, it appears from your initial post that your data are actually in the original EASE-Grid (using the authalic sphere as projection ellipsoid). Note that EPSG:6931 is for EASE-Grid 2.0. There were several early attempts to define EPSG codes for the original EASE-Grids. We would need to figure out which is the correct one for the data you have. There are plenty of data sets still in the original projections, so if there's a problem with the EPSG wkt, we should definitely figure that out. I am glad to help if you have more questions.

twhiteaker commented 4 years ago

Below is an expanded snippet as requested by @ethanrd. I followed @JonathanGregory's suggestion and switched the grid_mapping_name and projected_crs_name values. I also reformatted the crs_wkt to make it easier to read. If you plug in the value from the original post, the file should convert to netCDF. This file just has the bare minimum to cover the crs and x-y locations. It is not the fully compliant CF file, as I haven't gotten that far yet. Some history: the model spits out data in a custom text format, and I was attempting to convert to netCDF to make it easier for others to use. The model uses the old EASE Grid, not the newer 2.0 version. @marqh's comment prompted me to check out the WKT specs, which show to identify WKT built using an older version of the spec. So it looks like the WKT I'm using from EPSG.io is outdated. I'm not well-versed enough in WKT to update it presently.

netcdf from_arc_grid {
  dimensions:
    x = 39;
    y = 26;
  variables:
    double x(x);
      x:long_name = "x coordinate of projection";
      x:standard_name = "projection_x_coordinate";
      x:units = "Meter";

    int nsidc_ease_grid_north_25km;
      nsidc_ease_grid_north_25km:grid_mapping_name = "lambert_azimuthal_equal_area";
      nsidc_ease_grid_north_25km:longitude_of_projection_origin = 0.0;
      nsidc_ease_grid_north_25km:latitude_of_projection_origin = 90.0;
      nsidc_ease_grid_north_25km:false_easting = 0.0;
      nsidc_ease_grid_north_25km:false_northing = 0.0;
      nsidc_ease_grid_north_25km:geographic_crs_name = "GCS_Sphere_International_1924_Authalic";
      nsidc_ease_grid_north_25km:horizontal_datum_name = "Not_specified_based_on_International_1924_Authalic_Sphere";
      nsidc_ease_grid_north_25km:reference_ellipsoid_name = "International 1924 Authalic Sphere";
      nsidc_ease_grid_north_25km:semi_major_axis = 6371228.0;
      nsidc_ease_grid_north_25km:inverse_flattening = 0.0;
      nsidc_ease_grid_north_25km:prime_meridian_name = "Greenwich";
      nsidc_ease_grid_north_25km:longitude_of_prime_meridian = 0.0;
      nsidc_ease_grid_north_25km:projected_crs_name = "NSIDC_EASE_Grid_North";
      nsidc_ease_grid_north_25km:longitude_of_central_meridian  = 0.0;
      nsidc_ease_grid_north_25km:crs_wkt = 
        PROJCS["NSIDC EASE-Grid North",
          GEOGCS["Unspecified datum based upon the International 1924 Authalic Sphere",
            DATUM["Not_specified_based_on_International_1924_Authalic_Sphere",
              SPHEROID["International 1924 Authalic Sphere",6371228,0,AUTHORITY["EPSG","7057"]],
              AUTHORITY["EPSG","6053"]],
            PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
            UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],
            AUTHORITY["EPSG","4053"]],
          PROJECTION["Lambert_Azimuthal_Equal_Area"],
          PARAMETER["latitude_of_center",90],
          PARAMETER["longitude_of_center",0],
          PARAMETER["false_easting",0],
          PARAMETER["false_northing",0],
          UNIT["metre",1,AUTHORITY["EPSG","9001"]],
          AXIS["X",EAST],
          AXIS["Y",NORTH],
          AUTHORITY["EPSG","3408"]] ;

    float fake_var(y, x);
      fake_var:missing_value = -3.4028235E38f;
      fake_var:long_name = "fake_var";
      fake_var:coordinates = "x y";
      fake_var:grid_mapping = "nsidc_ease_grid_north_25km";
      fake_var:units = "m";

    double y(y);
      y:long_name = "y coordinate of projection";
      y:standard_name = "projection_y_coordinate";
      y:units = "Meter";

  // global attributes:
  :Conventions = "CF-1.7";
 data:
x =
  -1729986.943003976, -1704919.417143976, -1679851.8912839761, -1654784.3654239762, -1629716.8395639763, -1604649.3137039763, -1579581.7878439764, -1554514.2619839765, -1529446.7361239765, -1504379.2102639766, -1479311.6844039767, -1454244.1585439767, -1429176.6326839768, -1404109.1068239769, -1379041.580963977, -1353974.055103977, -1328906.529243977, -1303839.003383977, -1278771.4775239772, -1253703.9516639772, -1228636.4258039773, -1203568.8999439774, -1178501.3740839774, -1153433.8482239775, -1128366.3223639776, -1103298.7965039776, -1078231.2706439777, -1053163.7447839777, -1028096.2189239777, -1003028.6930639776, -977961.1672039776, -952893.6413439775, -927826.1154839775, -902758.5896239774, -877691.0637639774, -852623.5379039773, -827556.0120439773, -802488.4861839772, -777420.9603239772 ;
nsidc_ease_grid_north_25km = 1 ;

fake_var =

    1667.6692, 1600.6313, 1383.8802, 1242.0479, 1029.6685, 936.45776, 828.88995, 797.6116, 725.4621, 722.60846, 719.2757, 688.6377, 615.24744, 615.727, 589.234, 561.66626, 534.04517, 530.23157, 500.45145, 466.29483, 387.02902, 341.3668, 257.10144, 239.1657, 205.0, 205.0, 205.0, 205.0, 205.0, 205.0, 205.0, 205.0, 205.0, 205.0, 205.0, 205.0, 205.0, 205.01683, 205.01294,
    1671.1307, 1662.1647, 1563.2664, 1380.5562, 1239.3501, 1030.1027, 936.60925, 834.3774, 796.24316, 725.1207, 721.54004, 717.2908, 686.7519, 615.63, 616.3524, 591.1606, 565.6316, 539.508, 534.6084, 482.53415, 462.78027, 385.90796, 296.25885, 238.92244, 205.0, 205.0, 205.0, 205.0, 205.0, 205.0, 205.0, 205.0, 205.0, 205.0, 205.0, 205.0, 205.04216, 205.03383, 232.59724,
    1699.3593, 1665.6464, 1596.4927, 1543.5199, 1339.8746, 1137.5507, 1043.4258, 936.863, 834.43475, 783.63605, 708.3317, 720.1765, 714.58795, 684.04156, 616.0452, 617.15894, 593.993, 572.10065, 549.24896, 524.8121, 491.95132, 430.44058, 314.72443, 232.49597, 205.0, 205.0, 205.0, 205.0, 205.0, 205.0, 205.0, 205.0, 205.02841, 205.02841, 205.0, 205.12534, 327.90073, 314.9465, 301.0779,
    1703.5214, 1670.0665, 1667.3369, 1556.084, 1477.4391, 1284.2404, 1139.022, 1044.7275, 997.9504, 888.8829, 743.5883, 737.2321, 702.57874, 668.87213, 679.94946, 639.5634, 618.03265, 598.211, 601.6408, 569.72534, 546.22784, 509.46536, 316.63837, 205.06854, 205.0, 205.0, 205.0, 205.0, 205.0, 205.0, 205.02841, 205.08449, 205.19962, 205.20111, 205.18289, 205.47266, 692.5294, 471.48734, 383.9402,
    1708.6459, 1741.8413, 1662.6558, 1570.4044, 1538.1487, 1471.0963, 1288.844, 1194.92, 1094.5859, 998.59424, 889.0004, 797.40015, 737.9417, 685.5752, 665.4116, 673.4193, 635.08606, 637.3434, 603.70435, 597.77313, 617.4955, 813.4537, 205.38863, 205.09479, 205.0, 205.0, 205.0, 205.0, 205.0, 205.08449, 205.19962, 205.36646, 533.6088, 533.6418, 205.61722, 692.518, 692.78107, 582.65753, 513.4857,
    1762.441, 1747.5782, 1735.1123, 1661.5953, 1538.5953, 1530.1918, 1462.7273, 1343.358, 1257.9608, 1096.7577, 999.6101, 849.69305, 800.86365, 739.56274, 684.5412, 633.84485, 637.04254, 626.2405, 624.66675, 614.28253, 589.06635, 813.3728, 205.50923, 205.05258, 205.0, 205.0, 205.0, 205.0, 205.05608, 205.30888, 533.6088, 533.77856, 533.91626, 533.99756, 692.6615, 613.3648, 642.6252, 665.738, 664.3566,
    1834.098, 1712.6528, 1741.3107, 1667.1663, 1580.4484, 1529.6447, 1519.3242, 1329.19, 1344.3364, 1306.4607, 1160.8256, 1001.4501, 859.6937, 806.23724, 743.04, 683.1945, 630.8026, 625.4645, 596.94775, 589.03406, 589.07983, 813.42255, 205.68297, 205.21175, 205.0, 205.0, 205.0, 205.0, 205.16534, 533.6059, 533.8878, 533.97235, 534.0, 534.2034, 692.8236, 692.96436, 693.0, 693.0, 693.0,
    1836.3982, 1768.2614, 1749.5248, 1731.8792, 1585.8538, 1529.7852, 1516.9214, 1403.1202, 1403.422, 1418.7867, 1309.7406, 1174.5691, 1011.39764, 871.79865, 865.1643, 793.7436, 711.93207, 622.8396, 589.0, 589.0, 589.0563, 589.14294, 813.30554, 205.61856, 205.21175, 205.0, 205.0, 205.0, 205.21999, 533.6867, 533.9716, 534.0, 534.1917, 534.59717, 693.156, 693.05316, 693.0, 693.0, 693.0,
    1894.1924, 1840.0977, 1816.8716, 1740.9686, 1593.9713, 1541.0651, 1514.9049, 1447.4691, 1384.5281, 1404.901, 1420.7001, 1316.2999, 1179.4514, 1084.441, 946.4858, 918.61127, 799.68274, 681.6114, 589.0433, 589.0, 589.0, 589.1694, 813.6523, 813.42365, 205.63799, 205.15916, 205.0, 205.0, 205.16391, 533.59924, 533.914, 534.0, 534.56537, 1657.171, 693.88776, 693.25244, 693.0669, 693.0, 759.1413,
    1919.826, 1899.8077, 1850.1716, 1753.5385, 1727.1559, 1548.4854, 1523.2119, 1441.986, 1421.8821, 1384.8119, 1410.5525, 1427.5686, 1431.0393, 1193.6199, 1141.0457, 1024.4269, 963.42004, 829.0665, 589.2228, 589.0, 589.0, 589.1121, 813.72595, 813.7888, 813.3393, 205.37228, 205.05258, 205.0, 205.05464, 205.35693, 533.77045, 533.9716, 534.7522, 1658.0242, 1657.8549, 693.5905, 693.3924, 693.1943, 810.32245,
    1916.8677, 1898.6039, 1882.6769, 1823.1418, 1741.6816, 1578.5299, 1479.2343, 1436.4941, 1423.0979, 1406.2247, 1389.8658, 1521.7025, 1505.3613, 1446.0398, 1293.6143, 1147.3684, 1126.4141, 1212.5907, 589.7096, 589.2666, 589.0703, 589.10767, 589.40594, 813.843, 813.52216, 205.76042, 205.05258, 205.0, 205.54375, 205.21927, 533.6277, 533.914, 535.1601, 1658.9225, 1658.4161, 694.7344, 694.7762, 1293.1703, 1348.4988,
    1852.2799, 1923.5338, 1908.0535, 1900.9149, 1814.9977, 1651.266, 1538.8037, 1446.1113, 1408.2416, 1377.0646, 1400.6017, 1497.011, 1559.9365, 1587.1216, 1437.0786, 1301.3236, 1212.9928, 1212.9364, 1402.0242, 589.9254, 589.7066, 589.81696, 590.0911, 814.0539, 814.0952, 206.9128, 205.9906, 206.0732, 207.20384, 206.24074, 205.86145, 534.37305, 535.7426, 3492.8765, 1660.2894, 697.4469, 3493.158, 2280.8645, 1889.0016,
    1793.0905, 1890.3955, 1907.9243, 2008.6223, 1918.5328, 1822.6978, 1561.8844, 1456.3507, 1393.1827, 1340.623, 1351.8551, 1446.4163, 1603.3896, 1653.6628, 1614.7091, 1378.5635, 1301.9514, 1213.133, 1213.1909, 1402.5244, 1402.3212, 1402.9199, 2460.0369, 590.29816, 815.272, 3492.7468, 208.81946, 208.33824, 6492.856, 1759.349, 535.65643, 535.0456, 1539.8656, 3495.1882, 3019.1829, 6494.647, 3495.7432, 2846.4236, 2293.1934,
    1813.926, 1861.6526, 1903.226, 1973.5342, 2006.0195, 1853.0249, 1728.7468, 1549.7412, 1352.4491, 1266.0322, 1247.9811, 1329.8184, 1594.548, 1811.6636, 1793.6581, 1582.5719, 1302.1743, 1302.2859, 1214.49, 1214.5112, 1403.1327, 1403.5221, 2460.622, 590.78754, 815.7387, 3494.7717, 4389.8647, 2214.3916, 2771.652, 2498.2268, 4389.551, 1920.0123, 1750.7377, 2615.7295, 3145.381, 3749.2434, 3362.282, 2970.332, 2547.1814,
    1777.8163, 1859.0912, 1906.2814, 1974.4171, 2058.386, 1978.9233, 1848.4207, 1653.8483, 1379.5443, 1257.5254, 1139.1401, 1139.3469, 1624.6567, 1882.0488, 1965.9009, 1948.5388, 1302.638, 1303.1387, 3497.7556, 3496.811, 1214.8357, 1403.6735, 2461.0361, 590.55164, 815.0604, 1911.9514, 2621.6926, 2575.646, 2528.9116, 2376.4634, 2611.5413, 2224.9866, 2104.2786, 2577.6401, 3082.2375, 3085.8528, 3098.174, 2879.5906, 2688.4648,
    1764.4244, 1828.2345, 1888.253, 1952.5411, 2042.3264, 2238.2021, 2174.6716, 1799.267, 1543.2142, 1297.1641, 1139.2643, 1139.5691, 1881.4862, 2103.5586, 1966.0468, 1948.68, 1302.7925, 1303.8992, 3498.9954, 3497.8792, 1214.6901, 1646.5715, 1595.8993, 1246.7363, 1369.588, 1702.3994, 2112.3606, 2454.5334, 2478.3987, 2284.9697, 2275.6423, 2102.7283, 2429.367, 2465.4727, 2840.3179, 2825.8884, 2850.8464, 2784.49, 2665.8804,
    1718.6046, 1794.8638, 1877.3363, 1969.0251, 1997.7479, 2144.947, 2138.2874, 2099.7646, 1817.1332, 1139.8975, 1139.5016, 1140.075, 1881.9976, 2103.6313, 2103.7764, 1965.9026, 1949.2183, 1304.675, 2730.1772, 2765.9087, 2253.3982, 1949.9802, 1629.7665, 1353.1683, 1480.873, 1732.5785, 2051.753, 2378.937, 2446.097, 2206.1025, 2148.4534, 2368.8027, 2654.6575, 2429.1726, 2747.0376, 2904.7588, 2711.7297, 2700.2083, 2631.819,
    1685.0023, 1760.5691, 1888.5703, 1987.9873, 2146.2764, 1896.4961, 1896.7046, 2353.752, 2012.221, 2011.868, 1140.4604, 2462.5278, 1882.9994, 1882.6903, 2103.9534, 1966.4017, 1950.3428, 6500.7085, 3488.4758, 2878.9229, 2594.8652, 2204.802, 1776.953, 1542.9617, 1613.0967, 1745.0978, 2001.3356, 2390.0403, 2384.7021, 2126.9602, 2091.6455, 2231.2615, 2688.1504, 2398.241, 2688.7188, 2770.1914, 2788.3667, 2619.468, 2598.233,
    1601.6554, 1685.8271, 2040.5292, 2040.8624, 2278.431, 1897.0819, 1897.9906, 2354.286, 2013.2153, 3499.502, 1141.6252, 2463.6436, 4394.082, 2553.5327, 2185.3118, 2302.386, 2552.0383, 3491.221, 3341.5325, 2998.9543, 2713.1392, 2475.051, 2207.2266, 1585.8168, 1768.311, 1813.2218, 1985.244, 2347.249, 2317.3423, 2112.8645, 2295.0125, 2431.6033, 2706.4424, 2371.9866, 2648.929, 2715.8718, 2740.058, 2804.405, 2558.388,
    1453.1942, 1453.3436, 1453.7155, 2040.5779, 2278.7551, 1897.912, 6502.519, 2355.1963, 2013.2517, 3500.3945, 1141.3937, 2463.6196, 2676.9302, 2550.3787, 2360.5242, 2478.3606, 2636.4678, 2840.4897, 2919.1099, 2861.9917, 2671.157, 2470.295, 2312.865, 1811.0321, 1803.3445, 1842.2676, 1672.5413, 2324.2063, 2272.6497, 2087.9194, 2249.463, 2437.6487, 2716.3655, 2350.0679, 2619.9785, 2665.2917, 2705.6316, 2701.7534, 2762.2524,
    1453.3206, 1453.6749, 1453.901, 1454.434, 2041.37, 2279.8994, 3326.202, 3037.0286, 2558.9536, 2570.1707, 2165.9854, 2301.4426, 2452.9678, 2458.5056, 2324.5159, 2641.0852, 2646.7827, 2633.2043, 2729.1658, 2694.9475, 2661.8206, 2553.631, 2349.6697, 2034.855, 1826.9413, 1849.5292, 1816.9027, 2311.6462, 2378.8582, 2272.8325, 2228.3386, 2482.9895, 2918.7112, 2769.8408, 2597.9568, 2634.0454, 2667.2214, 2676.9048, 2744.7961,
    1893.442, 2465.146, 2465.5337, 2466.6724, 3149.46, 3732.2578, 3709.0935, 3418.0312, 3031.6562, 2534.7842, 2307.0352, 2333.139, 2390.268, 2457.0278, 2273.3552, 2671.3665, 2647.7324, 2599.4822, 2629.165, 2623.177, 2673.639, 2507.4556, 2390.5342, 2093.583, 2012.879, 1859.59, 1877.727, 2141.3372, 2151.4094, 2247.2808, 2278.3223, 2695.0815, 2927.9539, 2769.0503, 2580.6628, 2610.139, 2642.1035, 2657.2253, 2651.9216,
    2144.431, 2404.0876, 2919.3657, 4408.241, 4408.095, 6508.9556, 4535.4814, 3907.9006, 3415.207, 2950.5745, 2830.9404, 2377.2258, 2373.1028, 2435.0918, 2332.1184, 2617.5874, 2674.1086, 2632.4773, 2806.4966, 2573.2188, 2554.642, 2502.0178, 2303.2751, 2335.99, 1975.5845, 2166.781, 2112.0686, 1755.2678, 2134.9758, 2155.5896, 2256.1907, 2691.6484, 2933.6726, 2767.622, 2566.7485, 2591.3108, 2622.185, 2641.215, 2637.106,
    2533.9775, 2885.5159, 3337.9722, 3855.6895, 3505.0732, 3505.219, 4170.714, 4124.779, 3838.8967, 3420.1367, 2816.2097, 2521.3706, 2366.852, 2331.5688, 2321.1252, 2541.9646, 2598.6748, 2683.799, 2768.121, 2575.3772, 2533.7993, 2492.617, 2399.1335, 2360.7708, 2109.1323, 2200.8032, 2118.5393, 1748.1138, 2123.122, 2141.4182, 2239.1445, 2716.2654, 2979.0427, 2765.822, 2555.3345, 2576.1365, 2594.7053, 2613.7024, 2624.6196,
    2906.183, 3221.6003, 3732.7625, 4006.5708, 4415.4165, 4028.8335, 3988.1401, 4113.607, 4089.2922, 3489.6257, 3170.4888, 2841.285, 2542.6426, 2426.273, 2314.1702, 2631.608, 2623.216, 2670.9148, 2751.1199, 2768.3245, 2511.9236, 2508.8943, 2400.0405, 2265.8037, 2133.2334, 2227.4668, 2205.7126, 1985.6766, 2164.1272, 2130.689, 2225.8682, 2682.8875, 2949.8018, 3045.781, 2545.82, 2563.6738, 2579.8035, 2599.585, 2613.9438,
    3168.247, 3626.869, 3843.3628, 4115.968, 4418.7603, 4161.348, 4013.4546, 4105.967, 4108.43, 3583.6033, 3360.78, 3154.0063, 3019.2446, 2495.156, 2278.1326, 2624.5344, 2608.8618, 2587.8606, 2668.6638, 2848.7437, 2710.4495, 2496.5667, 2404.6682, 2329.613, 2149.9753, 2357.1348, 2226.0796, 2462.0925, 2244.6748, 2122.4116, 2387.0222, 2718.6172, 2945.7063, 3047.3606, 2986.2878, 2553.2761, 2567.414, 2587.7053, 2592.2214 ;

y =
  2205614.6174975215, 2180547.0916375215, 2155479.5657775216, 2130412.0399175216, 2105344.5140575217, 2080276.9881975218, 2055209.4623375218, 2030141.936477522, 2005074.410617522, 1980006.884757522, 1954939.358897522, 1929871.8330375222, 1904804.3071775222, 1879736.7813175223, 1854669.2554575223, 1829601.7295975224, 1804534.2037375225, 1779466.6778775225, 1754399.1520175226, 1729331.6261575227, 1704264.1002975227, 1679196.5744375228, 1654129.0485775229, 1629061.522717523, 1603993.996857523, 1578926.470997523 ;

}
twhiteaker commented 4 years ago

Hi @mjbrodzik, I'm think our data are in EPSG 3408. Our data use the original EASE Grid North (yep, authalic sphere-based) at 25 km resolution. At NSIDC's EASE Grid Page, you can scroll down and click the 25 km tab to see where they identify the projection as "Northern Hemisphere, Lambert Azimuthal ​(EPSG: 3408)". Representing this in modern WKT would be a step in the right direction, and I'm also trying to make sure I've represented it correctly using traditional CF attributes.

twhiteaker commented 4 years ago

I plan to follow the pattern in the publication that @mjbrodzik linked and include srid and proj4text attributes to augment the CF attributes and crs_wkt, along with an esri_pe_string (more on that below) to maximize likelihood that end users will be able to interpret the coordinate system. The outdated WKT format is related to WKT and not CF so that would be up to the data producer (my team in this case) to figure out rather than the CF team, though I would happily entertain a WKT guru reading this who wants to offer up the solution. So as long as nothing looks broken in the CF attributes, then I think this issue can be closed.

Note: I was testing by adding the netCDF file to ArcMap 10.6, and the grid cells were not aligning with known point data, which is what prompted me to open this issue. Upon closer inspection, it appears ArcMap is ignoring the inverse_flattening attribute, and assuming the WGS84 datum instead of using the authalic sphere. So I think it's ArcMap that is incorrect and not the netCDF file. I found that by including an esri_pe_string attribute on the data variable that is basically the WKT, the software would read that instead of the CF stuff and project the data correctly. 

fake_var:esri_pe_string = "PROJCS[\"NSIDC_EASE_Grid_North\",GEOGCS[\"GCS_Sphere_International_1924_Authalic\",DATUM[\"D_Sphere_International_1924_Authalic\",SPHEROID[\"Sphere_International_1924_Authalic\",6371228.0,0.0]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Lambert_Azimuthal_Equal_Area\"],PARAMETER[\"False_Easting\",0.0],PARAMETER[\"False_Northing\",0.0],PARAMETER[\"Central_Meridian\",0.0],PARAMETER[\"Latitude_Of_Origin\",90.0],UNIT[\"Meter\",1.0]]";
mjbrodzik commented 4 years ago

Wow, @twhiteaker, that esri_pe_string is a real find. When we first defined EASE-Grid and realized ESRI did that, (assumed WGS84) there was no remedy other than reprojecting to wgs84-based projection, which has always been a pain and was one of the motivators for defining EASE-Grid 2.0.

I do realize that including all that redundant information is fraught with risk for all the reasons people on this thread and related ones mentioned. In an ideal world we should only need to include this kind of technical information once, and not repeat ourselves. I do applaud those who are striving to create a more ideal mapping software world.

marqh commented 4 years ago

Hello @twhiteaker

using http://www.epsg-registry.org/export.htm?wkt=urn:ogc:def:crs:EPSG::6931 as a start point, we can adjust the datum definition and name as follows

(disclaimer, this is a draft, please check the contents)

PROJCRS["NSIDC EASE-Grid North",
  BASEGEODCRS["GCS_Sphere_International_1924_Authalic",
    DATUM["D_Sphere_International_1924_Authalic",
      ELLIPSOID["Sphere_International_1924_Authalic",6371228.0,0.0,LENGTHUNIT["metre",1.0]]],
    PRIMEM["Greenwich",0.0]],
  CONVERSION["US NSIDC EASE-Grid North",
    METHOD["Lambert Azimuthal Equal Area",ID["EPSG",9820]],
    PARAMETER["Latitude of natural origin",90,ANGLEUNIT["degree",0.01745329252]],
    PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.01745329252]],
    PARAMETER["False easting",0,LENGTHUNIT["metre",1.0]],
    PARAMETER["False northing",0,LENGTHUNIT["metre",1.0]]],
  CS[cartesian,2],
    AXIS["easting (X)",south,MERIDIAN[90,ANGLEUNIT["degree",0.01745329252]],ORDER[1]],
    AXIS["northing (Y)",south,MERIDIAN[180,ANGLEUNIT["degree",0.01745329252]],ORDER[2]],
    LENGTHUNIT["metre",1.0]]

note, you can check syntax for errors with a modern version of proj (e.g.proj-6.2.1) 

projinfo 'PROJCRS["NSIDC EASE-Grid North",
>   BASEGEODCRS["GCS_Sphere_International_1924_Authalic",
>     DATUM["D_Sphere_International_1924_Authalic",
>       ELLIPSOID["Sphere_International_1924_Authalic",6371228.0,0.0,LENGTHUNIT["metre",1.0]]],
>     PRIMEM["Greenwich",0.0]],
>   CONVERSION["US NSIDC EASE-Grid North",
>     METHOD["Lambert Azimuthal Equal Area",ID["EPSG",9820]],
>     PARAMETER["Latitude of natural origin",90,ANGLEUNIT["degree",0.01745329252]],
>     PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.01745329252]],
>     PARAMETER["False easting",0,LENGTHUNIT["metre",1.0]],
>     PARAMETER["False northing",0,LENGTHUNIT["metre",1.0]]],
>   CS[cartesian,2],
>     AXIS["easting (X)",south,MERIDIAN[90,ANGLEUNIT["degree",0.01745329252]],ORDER[1]],
>     AXIS["northing (Y)",south,MERIDIAN[180,ANGLEUNIT["degree",0.01745329252]],ORDER[2]],
>     LENGTHUNIT["metre",1.0]]'
PROJ.4 string:
+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +R=6371228 +units=m +no_defs +type=crs

(though also note that projinfo also outputs it's own adapted WKT, which may have ID entries that you do not want)

I hope this helps marqh

mjbrodzik commented 4 years ago

Thanks to @marqh I am learning all kinds of useful tricks, I didn't know about projinfo, that is a very useful tool.

One note to @twhiteaker, please do not start referring to the original EASE-Grids with anything with a version number, like "version 1.0 " or "1.0". When we wrote the definition for EASE-Grid 2.0, we realized our mistake in not giving a version to the original one ("whyever would anyone need to revise such a thing?"). People had been using the original EASE-Grids for almost 20 years and the definitions had never included a version number. We decided that the original definitions should be referred to with no version number, and that the newer version should always be referred to as "EASE-Grid 2.0", and then left room for potential newer revision numbers after that. Respectfully, Mary Jo

marqh commented 4 years ago

please do not start referring to the original EASE-Grids with anything with a version number, like "version 1.0 " or "1.0".

Hello @mjbrodzik

I have edited my previous comment and removed instances of 1.0, to try to limit potential confusion I hope this helps

marqh

mjbrodzik commented 4 years ago

Thanks to @marqh ! Keeping terminology consistent is critical.

twhiteaker commented 4 years ago

Thanks @marqh for the projinfo tip. I found that you could enter

projinfo EPSG:3408

to get this (with some IDs removed). Maybe that's what you were referring to regarding proj's own adapted WKT.

PROJCRS["NSIDC EASE-Grid North",
    BASEGEOGCRS["Unspecified datum based upon the International 1924 Authalic Sphere",
        DATUM["Not specified (based on International 1924 Authalic Sphere)",
            ELLIPSOID["International 1924 Authalic Sphere",6371228,0,LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]]],
    CONVERSION["US NSIDC Equal Area north projection",
        METHOD["Lambert Azimuthal Equal Area (Spherical)"],
        PARAMETER["Latitude of natural origin",90,ANGLEUNIT["degree",0.0174532925199433]],
        PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433]],
        PARAMETER["False easting",0,LENGTHUNIT["metre",1]],
        PARAMETER["False northing",0,LENGTHUNIT["metre",1]]],
    CS[Cartesian,2],
        AXIS["easting (X)",south,MERIDIAN[90,ANGLEUNIT["degree",0.0174532925199433]],ORDER[1],LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",south,MERIDIAN[180,ANGLEUNIT["degree",0.0174532925199433]],ORDER[2],LENGTHUNIT["metre",1]],
    ID["EPSG",3408]]
snowman2 commented 4 years ago

Side note that pyproj may be useful when converting to CF: https://pyproj4.github.io/pyproj/latest/build_crs_cf.html