Changeset 51630


Ignore:
Timestamp:
May 15, 2012 3:56:21 PM (4 years ago)
Author:
huhabla
Message:

Added spatial coordinates reference system metadata to netcdf files
to fulfill CF-1.6 and gdal needs.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • grass/trunk/raster3d/r3.out.netcdf/main.c

    r51625 r51630  
    3939#define TIME_NAME "time"
    4040#define X_NAME "x"
    41 #define X_LONG_NAME "x-axsis coordinates"
     41#define X_STANDARD_NAME "projection_x_coordinate"
     42#define X_LONG_NAME "x coordinate of projection"
    4243#define Y_NAME "y"
    43 #define Y_LONG_NAME "y-axsis coordinates"
     44#define Y_LONG_NAME "y coordinate of projection"
     45#define Y_STANDARD_NAME "projection_y_coordinate"
    4446#define Z_NAME "z"
    45 #define Z_LONG_NAME "z-axsis coordinates"
     47#define Z_LONG_NAME "z coordinate of projection"
    4648#define UNITS "units"
    4749#define DEGREES_EAST "degrees_east"
    4850#define DEGREES_NORTH "degrees_north"
     51#define HISTORY_TEXT "GRASS GIS 7 NetCDF export of r3.out.netcdf"
     52#define CF_SUPPORT "CF-1.5"
    4953
    5054#define ERR(e) {fatalError(nc_strerror(e));}
     
    99103    int lat_dimid = 0, lon_dimid = 0, time_dimid = 0;
    100104    int lat_varid = 0, lon_varid = 0, time_varid = 0;
    101     struct Key_Value *pkv, *ukv;
    102105    int dimids[NDIMS];
    103106    struct Cell_head window;
    104    
     107
     108    /* global attributes */
     109    if ((retval = nc_put_att_text(ncid, NC_GLOBAL, "Conventions", strlen(CF_SUPPORT), CF_SUPPORT))) ERR(retval);
     110    if ((retval = nc_put_att_text(ncid, NC_GLOBAL, "history", strlen(HISTORY_TEXT), HISTORY_TEXT))) ERR(retval);
     111
    105112    G_get_window(&window);
    106    
     113     
     114    /* Projection parameter */
     115    if(window.proj != PROJECTION_XY) {
     116        struct Key_Value *pkv, *ukv;
     117        struct pj_info pjinfo;
     118        char *proj4, *proj4mod;
     119        const char *unfact;
     120        int crs_dimid = 0, crs_varid;
     121
     122        if ((retval = nc_def_var(ncid, "crs", NC_CHAR, 0, &crs_dimid, &crs_varid))) ERR(retval);
     123
     124        pkv = G_get_projinfo();
     125        ukv = G_get_projunits();
     126 
     127        pj_get_kv(&pjinfo, pkv, ukv);
     128        proj4 = pj_get_def(pjinfo.pj, 0);
     129        pj_free(pjinfo.pj);
     130
     131#ifdef HAVE_OGR
     132        /* We support the CF suggestion crs_wkt and the gdal spatil_ref attribute */
     133        if ((retval = nc_put_att_text(ncid, crs_varid, "crs_wkt", strlen(GPJ_grass_to_wkt(pkv, ukv, 0, 0)), GPJ_grass_to_wkt(pkv, ukv, 0, 0)))) ERR(retval);
     134        if ((retval = nc_put_att_text(ncid, crs_varid, "spatial_ref", strlen(GPJ_grass_to_wkt(pkv, ukv, 0, 0)), GPJ_grass_to_wkt(pkv, ukv, 0, 0)))) ERR(retval);
     135#endif
     136        /* Code from g.proj:
     137         * GRASS-style PROJ.4 strings don't include a unit factor as this is
     138         * handled separately in GRASS - must include it here though */
     139        unfact = G_find_key_value("meters", ukv);
     140        if (unfact != NULL && (strcmp(pjinfo.proj, "ll") != 0))
     141            G_asprintf(&proj4mod, "%s +to_meter=%s", proj4, unfact);
     142        else
     143            proj4mod = G_store(proj4);
     144        pj_dalloc(proj4);
     145
     146        if ((retval = nc_put_att_text(ncid, crs_varid, "crs_proj4", strlen(proj4mod), proj4mod))) ERR(retval);
     147
     148        if(pkv)
     149            G_free_key_value(pkv);
     150        if(ukv)
     151            G_free_key_value(ukv);
     152    }
     153
    107154    typeIntern = Rast3d_tile_type_map(map);
    108155   
     
    129176        if ((retval = nc_put_att_text(ncid, lon_varid, UNITS, strlen("meter"), "meter"))) ERR(retval);
    130177        if ((retval = nc_put_att_text(ncid, lon_varid, LONG_NAME, strlen(X_LONG_NAME), X_LONG_NAME))) ERR(retval);
    131         if ((retval = nc_put_att_text(ncid, lon_varid, STANDARD_NAME, strlen(X_NAME), X_NAME))) ERR(retval);
     178        if ((retval = nc_put_att_text(ncid, lon_varid, STANDARD_NAME, strlen(X_STANDARD_NAME), X_STANDARD_NAME))) ERR(retval);
    132179        if ((retval = nc_put_att_text(ncid, lat_varid, UNITS, strlen("meter"), "meter"))) ERR(retval);
    133180        if ((retval = nc_put_att_text(ncid, lat_varid, LONG_NAME, strlen(Y_LONG_NAME), Y_LONG_NAME))) ERR(retval);
    134         if ((retval = nc_put_att_text(ncid, lat_varid, STANDARD_NAME, strlen(Y_NAME), Y_NAME))) ERR(retval);
     181        if ((retval = nc_put_att_text(ncid, lat_varid, STANDARD_NAME, strlen(Y_STANDARD_NAME), Y_STANDARD_NAME))) ERR(retval);
    135182    }
    136183       
     
    176223        if ((retval = nc_put_att_text(ncid, time_varid, LONG_NAME, strlen(Z_LONG_NAME), Z_LONG_NAME))) ERR(retval);
    177224    }
    178     /* z - axsis orientation */
     225    /* z - axis orientation */
    179226    if ((retval = nc_put_att_text(ncid, time_varid, "positive", strlen("up"), "up"))) ERR(retval);
    180227   
     228    /* Axis identifier attributes */
     229    if ((retval = nc_put_att_text(ncid, lon_varid, "axis", 1, "X"))) ERR(retval);
     230    if ((retval = nc_put_att_text(ncid, lat_varid, "axis", 1, "Y"))) ERR(retval);
     231    if(is_time) {
     232        if ((retval = nc_put_att_text(ncid, time_varid, "axis", 1, "T"))) ERR(retval);
     233    } else {
     234        if ((retval = nc_put_att_text(ncid, time_varid, "axis", 1, "Z"))) ERR(retval);
     235    }
     236
    181237    dimids[0] = time_dimid;
    182238    dimids[1] = lat_dimid;
     
    188244        if ((retval = nc_def_var(ncid, param.input->answer, NC_DOUBLE, NDIMS, dimids, varid))) ERR(retval);
    189245    }
     246    if(window.proj != PROJECTION_XY) {
     247        if ((retval = nc_put_att_text(ncid, *varid, "grid_mapping", strlen("crs"), "crs"))) ERR(retval);
     248    }
    190249   
    191250   /* End define mode. */
    192251   if ((retval = nc_enddef(ncid))) ERR(retval); 
    193252   
    194     /* Build coordinates, we need to use the cell center in case of spatial dimensions */
     253    /*
     254     * Build coordinates, we need to use the cell center in case of spatial dimensions
     255     * */
     256
    195257    for(i = 0; i < region->cols; i++) {
    196258        c = region->west + i*region->ew_res + 0.5*region->ew_res;
     
    204266   
    205267    for(i = 0; i < region->depths; i++) {
    206         if(is_time)
     268        if(is_time) {
    207269            c = region->bottom + i*region->tb_res;
    208         else
     270            time = (int)c;
     271            nc_put_var1_int(ncid, time_varid, &i, &time); 
     272        } else {
    209273            c = region->bottom + i*region->tb_res + 0.5*region->tb_res;
    210         time = (int)c;
    211         if(is_time)
    212             nc_put_var1_int(ncid, time_varid, &i, &time); 
    213         else
    214274            nc_put_var1_float(ncid, time_varid, &i, &c);     
    215     }
    216    
    217     /* Lets try to specify the projection information */
    218     pkv = G_get_projinfo();
    219     ukv = G_get_projunits();
    220    
    221     if(pkv && ukv) {
    222         for(i = 0; i < pkv->nitems; i++)
    223             fprintf(stderr, "%s : %s\n", pkv->key[i], pkv->value[i]);
    224         for(i = 0; i < ukv->nitems; i++)
    225             fprintf(stderr, "%s : %s\n", ukv->key[i], ukv->value[i]); 
    226     }
    227    
    228     if(window.proj != PROJECTION_XY)
    229         ; // Do projection
    230    
    231     if(pkv)
    232         G_free_key_value(pkv);
    233     if(ukv)
    234         G_free_key_value(ukv);
     275        }
     276    }
    235277}
    236278
Note: See TracChangeset for help on using the changeset viewer.