Commit 7ac7169c authored by Dan Baston's avatar Dan Baston

Use OGC WKT for spatial_ref crs attribute

This change is being made in response to the issue described at
#74, and modifies the
spatial_ref attribute to match the one generated using GDAL. This does
not change the inability of the R raster package to understand this
field, which is a GDAL construct and not incorporated into CF
conventions.
parent aaa3fa64
Pipeline #73783178 passed with stages
in 23 minutes and 59 seconds
......@@ -58,7 +58,7 @@ ncks --no_tmp_fl -h -C -O -x -v time $TEMP_NC2 $TEMP_NC3
# Add a CRS variable
ncap2 --no_tmp_fl -h -O -s 'crs=-9999' $TEMP_NC3 $2
ncatted -h -O \
-a spatial_ref,crs,c,c,'GEOGCS[\"GCS_WGS_1984\",DATUM[\"WGS_1984\",SPHEROID[\"WGS_84\",6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.017453292519943295]]' \
-a spatial_ref,crs,c,c,'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]' \
-a grid_mapping_name,crs,c,c,'latitude_longitude' \
-a longitude_of_prime_meridian,crs,c,d,0 \
-a semi_major_axis,crs,c,d,6378137 \
......
......@@ -426,7 +426,7 @@ write_wgs84_crs_attributes <- function(ncout, var_names) {
ncdf4::ncatt_put(ncout, "crs", "longitude_of_prime_meridian", 0.0)
ncdf4::ncatt_put(ncout, "crs", "semi_major_axis", 6378137.0)
ncdf4::ncatt_put(ncout, "crs", "inverse_flattening", 298.257223563)
ncdf4::ncatt_put(ncout, "crs", "spatial_ref", "GEOGCS[\"GCS_WGS_1984\",DATUM[\"WGS_1984\",SPHEROID[\"WGS_84\",6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.017453292519943295]]")
ncdf4::ncatt_put(ncout, "crs", "spatial_ref", 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]')
for (var in var_names) {
ncdf4::ncatt_put(ncout, var, "grid_mapping", "crs")
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment