diff --git a/xval.d/cmd_prcp_ndays_filter_nc.c b/xval.d/cmd_prcp_ndays_filter_nc.c new file mode 100644 index 0000000..7edba66 --- /dev/null +++ b/xval.d/cmd_prcp_ndays_filter_nc.c @@ -0,0 +1,370 @@ +/* + cmd_prcp_ndays_filter_nc.c + Peter Thornton + 10/12/09 + + Since station-years are not currently being filtered as they are + extracted from the portal for backend processing, a new routine is + introduced here that checks the number of missing days for each station-year + and throws out the stations that don't meet the user-specfied criteria. + + */ + +#include "metsrc2.h" +#include "netcdf.h" + +int main(int argc, char *argv[]) +{ + /* variable declarations */ + + file prcp_nc_f; + + char *good, *good2; + char *listgood; + + short count; + + int ndays,nstns,nregr; + int i,j,k,l; + int day; + int x,y; + int type; + int nrows, ncols; + int width, hw; + int stn_filex, stn_filey; + int test; + int smwidth; + int ngoodpts; + int *keep; + int nkeep,nfiltgood; + + long off; + long int obs_off, obs_offset; + long stn_offset; + + double pred, obs; + double rn; + double dfirad; + double ul_mapx, ul_mapy; + double cellsize; + double stn_mapx, stn_mapy; + double f_x, f_y; + double value,sum_wt; + + double b, ae; + double critp; + + size_t lenp; + int xyloc_flag; + int status, ncid; + int nstnid, descid, yrdid; + int no_dataid; + int in_tileid; + int stnnameid, stnidid, stnelevid, stnlatid, stnlonid; + int stnxid, stnyid, stnzid, goodid; + int stnfxid, stnfyid, stnfzid; + int stntypeid, stnuc1id, stnuc2id; + int prcpid; + int prcp_name; + int ncid_meta; + int year; + int year_day, descriptor_length; + char *stnnames, *stnids; + double *elevs, *lats, *lons, *prcp; + double *stnx, *stnfx, *stny, *stnfy; + double *no_data; + int *in_tile; + char *stnnames_filt, *stnids_filt; + double *elevs_filt, *lats_filt, *lons_filt, *prcp_filt; + double *stnx_filt, *stnfx_filt, *stny_filt, *stnfy_filt; + double *no_data_filt; + int *in_tile_filt; + int *mflag, *mflag_filt; + int mflagid, mflagdims[2]; + int ii, prcpdims[2], stndims[2], gooddims[2]; + int codeid, codedims[2]; + + /* check for command line arguments */ + if (argc != 3) + { + printf("Usage: \n"); + printf("argc = %d\n",argc); + exit(1); + } + /* get initialization data from command line */ + strcpy(prcp_nc_f.name,argv[1]); + critp = atof(argv[2]); +/* ----------------------- netcdf calls for station meta file --------------------------------------- */ +/* */ + if( status = nc_open(prcp_nc_f.name, NC_WRITE, &ncid_meta) ) + { + printf("Error opening %s for netcdf read, exiting\n", prcp_nc_f.name); + exit(1); + } + status = nc_inq_dimid (ncid_meta, "stations", &nstnid); + status = nc_inq_dimlen(ncid_meta, nstnid, &lenp); + nstns = lenp; + status = nc_inq_dimid (ncid_meta, "year_day", &yrdid); + status = nc_inq_dimlen(ncid_meta, yrdid, &lenp); + year_day = lenp; + status = nc_inq_dimid (ncid_meta, "descriptor_length", &descid); + status = nc_inq_dimlen(ncid_meta, descid, &lenp); + descriptor_length = lenp; + status = nc_get_att_int (ncid_meta, NC_GLOBAL, "year", &year); + status = nc_get_att_int (ncid_meta, NC_GLOBAL, "days_per_year", &year_day); + ndays = year_day; + + status = nc_inq_varid (ncid_meta, "station_name", &stnnameid); + if (!(stnnames = (char*) malloc(nstns * descriptor_length * sizeof(char)))) exit(1); + status = nc_get_var_text (ncid_meta, stnnameid, stnnames); + + status = nc_inq_varid (ncid_meta, "station_id", &stnidid); + if (!(stnids = (char*) malloc(nstns * descriptor_length * sizeof(char)))) exit(1); + status = nc_get_var_text (ncid_meta, stnidid, stnids); + + status = nc_inq_varid (ncid_meta, "station_elevation", &stnelevid); + if (!(elevs = (double*) malloc(nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnelevid, elevs); + + status = nc_inq_varid (ncid_meta, "station_latitude", &stnlatid); + if (!(lats = (double*) malloc(nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnlatid, lats); + + status = nc_inq_varid (ncid_meta, "station_longitude", &stnlonid); + if (!(lons = (double*) malloc(nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnlonid, lons); + + status = nc_inq_varid (ncid_meta, "no_data", &no_dataid); + if (!(no_data = (double*) malloc(nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, no_dataid, no_data); + + status = nc_inq_varid (ncid_meta, "in_tile", &in_tileid); + if (!(in_tile = (int*) malloc(nstns * sizeof(int)))) exit(1); + status = nc_get_var_int (ncid_meta, in_tileid, in_tile); + + status = nc_inq_varid (ncid_meta, "stnx", &stnxid); + if (!(stnx = (double*) malloc(nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnxid, stnx); + + status = nc_inq_varid (ncid_meta, "stny", &stnyid); + if (!(stny = (double*) malloc(nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnyid, stny); + + status = nc_get_att_int (ncid_meta, NC_GLOBAL, "xyloc_flag", &xyloc_flag); + if (xyloc_flag == 1) + { + status = nc_inq_varid (ncid_meta, "stnfx", &stnfxid); + if (!(stnfx = (double*) malloc(nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnfxid, stnfx); + + status = nc_inq_varid (ncid_meta, "stnfy", &stnfyid); + if (!(stnfy = (double*) malloc(nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnfyid, stnfy); + } + + status = nc_inq_varid (ncid_meta, "precip", &prcpid); + if (!(prcp = (double*) malloc(nstns * year_day * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, prcpid, prcp); + + nc_close(ncid_meta); +/* */ +/* ------------------- end netcdf calls for station meta file --------------------------------------- */ + /* allocate memory */ + if (!(good = (char*) malloc(nstns * ndays * sizeof(char)))) exit(1); + if (!(mflag = (int*) malloc(nstns * ndays * sizeof(int)))) exit(1); + if (!(keep = (int*) malloc(nstns * sizeof(int)))) exit(1); + /* read station daily data and fill appropriate arrays */ + for (i=0 ; i= critp) + { + /* this station passes the filtering criteria, and will be + retained in the filtered stnmeta and stndata files */ + keep[i] = 1; + nkeep++; + } + } + /* allocate output memory */ + if (!(stnnames_filt = (char*) malloc(nkeep * descriptor_length * sizeof(char)))) exit(1); + if (!(stnids_filt = (char*) malloc(nkeep * descriptor_length * sizeof(char)))) exit(1); + if (!(elevs_filt = (double*) malloc(nkeep * sizeof(double)))) exit(1); + if (!(lats_filt = (double*) malloc(nkeep * sizeof(double)))) exit(1); + if (!(lons_filt = (double*) malloc(nkeep * sizeof(double)))) exit(1); + if (!(no_data_filt = (double*) malloc(nkeep * sizeof(double)))) exit(1); + if (!(in_tile_filt = (int*) malloc(nkeep * sizeof(int)))) exit(1); + if (!(stnx_filt = (double*) malloc(nkeep * sizeof(double)))) exit(1); + if (!(stny_filt = (double*) malloc(nkeep * sizeof(double)))) exit(1); + if (xyloc_flag) + { + if (!(stnfx_filt = (double*) malloc(nkeep * sizeof(double)))) exit(1); + if (!(stnfy_filt = (double*) malloc(nkeep * sizeof(double)))) exit(1); + } + if (!(prcp_filt = (double*) malloc(nkeep * year_day * sizeof(double)))) exit(1); + if (!(mflag_filt = (int*) malloc(nkeep * year_day * sizeof(int)))) exit(1); + + /* read through the meta and daily data, and write out the station records that + are being kept to the new meta and data output files */ + ii = 0; + for (i=0 ; i \n"); + exit(1); + } + + /* scan command line arguments into parameters */ + job = atoi(argv[1]); + strcpy(mask_nc_f.name, argv[2]); + strcpy(dem_nc_f.name, argv[3]); + strcpy(stnmeta_nc_f.name, argv[4]); + strcpy(stndata_nc_f.name, argv[5]); + fzwidth = atof(argv[6]); + dfirad = atof(argv[7]); + init_gsp = atof(argv[8]); + intstr.ans = atof(argv[9]); + prcpstr.f_max = atof(argv[10]); + pop_crit = atof(argv[11]); + smwidth = atoi(argv[12]); + inflag = atoi(argv[13]); + inval = atoi(argv[14]); + strcpy(outprefix,argv[15]); + ans=(int)intstr.ans; + + /* DEBUG SNOTEL CODE ADDED 8/29/00 by PET */ + if (DEBUG_SNOTEL) + { + strcpy(debug_snotel_f.name,"debug_stndat.txt"); + if (file_open(&debug_snotel_f,'w')) + { + printf("Error opening debug snotel file, exiting\n"); + exit(1); + } + } + + /* create netcdf xval output file */ + strcpy(out_f.name, outprefix); + strcat(out_f.name,"_xval.nc"); + if (status = nc_create(out_f.name, 0, &ncid_out)) + { + printf("Error creating %s as netcdf file, exiting\n", out_f.name); + exit(1); + } + /* define scalar variables for xval output file */ + if (status = nc_def_var(ncid_out, "tileid", NC_INT, 0, 0, &outid_tileid)) + { + printf("Error creating tileid variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"Daymet tile ID"); + nc_put_att_text(ncid_out, outid_tileid, "longname", strlen(att), att); + strcpy(att,"index"); + nc_put_att_text(ncid_out, outid_tileid, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "nstns3x3", NC_INT, 0, 0, &outid_nstns3x3)) + { + printf("Error creating nstns3x3 variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"number of stations in 3x3 tile region"); + nc_put_att_text(ncid_out, outid_nstns3x3, "longname", strlen(att), att); + strcpy(att,"stations"); + nc_put_att_text(ncid_out, outid_nstns3x3, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "nstns", NC_INT, 0, 0, &outid_nstns)) + { + printf("Error creating nstns variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"number of stations evaluated (central tile)"); + nc_put_att_text(ncid_out, outid_nstns, "longname", strlen(att), att); + strcpy(att,"stations"); + nc_put_att_text(ncid_out, outid_nstns, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "nstndays", NC_INT, 0, 0, &outid_nstndays)) + { + printf("Error creating nstndays variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"number of station-days evaluated"); + nc_put_att_text(ncid_out, outid_nstndays, "longname", strlen(att), att); + strcpy(att,"days"); + nc_put_att_text(ncid_out, outid_nstndays, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "rad90mean", NC_DOUBLE, 0, 0, &outid_rad90mean)) + { + printf("Error creating rad90mean variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"mean: radius capturing 90% of filter kernel weight"); + nc_put_att_text(ncid_out, outid_rad90mean, "longname", strlen(att), att); + strcpy(att,"m"); + nc_put_att_text(ncid_out, outid_rad90mean, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "rad90stdv", NC_DOUBLE, 0, 0, &outid_rad90stdv)) + { + printf("Error creating rad90stdv variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"std dev: radius capturing 90% of filter kernel weight"); + nc_put_att_text(ncid_out, outid_rad90stdv, "longname", strlen(att), att); + strcpy(att,"m"); + nc_put_att_text(ncid_out, outid_rad90stdv, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "daymae", NC_DOUBLE, 0, 0, &outid_daymae)) + { + printf("Error creating daymae variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"mean absolute error for single-day predictions"); + nc_put_att_text(ncid_out, outid_daymae, "longname", strlen(att), att); + strcpy(att,"cm/day"); + nc_put_att_text(ncid_out, outid_daymae, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "pormae", NC_DOUBLE, 0, 0, &outid_pormae)) + { + printf("Error creating annmae variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"mean absolute error for period-of-record predictions"); + nc_put_att_text(ncid_out, outid_pormae, "longname", strlen(att), att); + strcpy(att,"cm/day"); + nc_put_att_text(ncid_out, outid_pormae, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "pormpae", NC_DOUBLE, 0, 0, &outid_pormpae)) + { + printf("Error creating annmpae variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"mean absolute error as a percentage, for period-of-record predictions"); + nc_put_att_text(ncid_out, outid_pormpae, "longname", strlen(att), att); + strcpy(att,"%"); + nc_put_att_text(ncid_out, outid_pormpae, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "bias", NC_DOUBLE, 0, 0, &outid_bias)) + { + printf("Error creating bias variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"mean prediction bias"); + nc_put_att_text(ncid_out, outid_bias, "longname", strlen(att), att); + strcpy(att,"cm/day"); + nc_put_att_text(ncid_out, outid_bias, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "ppmean", NC_DOUBLE, 0, 0, &outid_ppmean)) + { + printf("Error creating ppmean variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"mean observed precipitation"); + nc_put_att_text(ncid_out, outid_ppmean, "longname", strlen(att), att); + strcpy(att,"cm/day"); + nc_put_att_text(ncid_out, outid_ppmean, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "xlrmean", NC_DOUBLE, 0, 0, &outid_xlrmean)) + { + printf("Error creating xlrmean variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"3-d regression: mean x-component"); + nc_put_att_text(ncid_out, outid_xlrmean, "longname", strlen(att), att); + strcpy(att,"1/m"); + nc_put_att_text(ncid_out, outid_xlrmean, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "xlrstdv", NC_DOUBLE, 0, 0, &outid_xlrstdv)) + { + printf("Error creating xlrstdv variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"3-d regression: among-station std dev of x-component"); + nc_put_att_text(ncid_out, outid_xlrstdv, "longname", strlen(att), att); + strcpy(att,"1/m"); + nc_put_att_text(ncid_out, outid_xlrstdv, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "ylrmean", NC_DOUBLE, 0, 0, &outid_ylrmean)) + { + printf("Error creating ylrmean variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"3-d regression: mean y-component"); + nc_put_att_text(ncid_out, outid_ylrmean, "longname", strlen(att), att); + strcpy(att,"1/m"); + nc_put_att_text(ncid_out, outid_ylrmean, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "ylrstdv", NC_DOUBLE, 0, 0, &outid_ylrstdv)) + { + printf("Error creating ylrstdv variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"3-d regression: among-station std dev of y-component"); + nc_put_att_text(ncid_out, outid_ylrstdv, "longname", strlen(att), att); + strcpy(att,"1/m"); + nc_put_att_text(ncid_out, outid_ylrstdv, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "zlrmean", NC_DOUBLE, 0, 0, &outid_zlrmean)) + { + printf("Error creating zlrmean variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"3-d regression: mean z-component"); + nc_put_att_text(ncid_out, outid_zlrmean, "longname", strlen(att), att); + strcpy(att,"1/m"); + nc_put_att_text(ncid_out, outid_zlrmean, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "zlrstdv", NC_DOUBLE, 0, 0, &outid_zlrstdv)) + { + printf("Error creating zlrstdv variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"3-d regression: among-station std dev of z-component"); + nc_put_att_text(ncid_out, outid_zlrstdv, "longname", strlen(att), att); + strcpy(att,"1/m"); + nc_put_att_text(ncid_out, outid_zlrstdv, "units", strlen(att), att); + + + +/* ----------------------- netcdf calls for mask file --------------------------------------- */ +/* */ + strcat(mask_nc_f.name,".nc"); + if( status = nc_open(mask_nc_f.name, 0, &ncid_mask) ) + { + printf("Error opening %s for netcdf read, exiting\n", mask_nc_f.name); + exit(1); + } + status = nc_inq_dimid (ncid_mask, "ncols", &colsid_mask); + status = nc_inq_dimlen(ncid_mask, colsid_mask, &lenp); + mask_nc_hdr.ncols = lenp; + status = nc_inq_dimid (ncid_mask, "nrows", &rowsid_mask); + status = nc_inq_dimlen(ncid_mask, rowsid_mask, &lenp); + mask_nc_hdr.nrows = lenp; + status = nc_get_att_text (ncid_mask, NC_GLOBAL, "tag", mask_nc_hdr.tag); + status = nc_get_att_double (ncid_mask, NC_GLOBAL, "ulx", &mask_nc_hdr.ulx); + status = nc_get_att_double (ncid_mask, NC_GLOBAL, "uly", &mask_nc_hdr.uly); + status = nc_get_att_double (ncid_mask, NC_GLOBAL, "cellsize", &mask_nc_hdr.cellsize); + status = nc_get_att_long (ncid_mask, NC_GLOBAL, "gctp_outsys", &mask_nc_hdr.gctp.outsys); + status = nc_get_att_long (ncid_mask, NC_GLOBAL, "gctp_outzone", &mask_nc_hdr.gctp.outzone); + status = nc_get_att_long (ncid_mask, NC_GLOBAL, "gctp_outdatum", &mask_nc_hdr.gctp.outdatum); + status = nc_get_att_double (ncid_mask, NC_GLOBAL, "gctp_outparm", mask_nc_hdr.gctp.outparm); + status = nc_get_att_text (ncid_mask, NC_GLOBAL, "gctp_fn27", mask_nc_hdr.gctp.fn27); + status = nc_get_att_text (ncid_mask, NC_GLOBAL, "gctp_fn83", mask_nc_hdr.gctp.fn83); + if (strncmp(mask_nc_hdr.tag,IMGTAG,7)) + { + printf("Error: mask file not a registered image. Use register_image first.\n"); + exit(1); + } + + status = nc_inq_varid (ncid_mask, "image", &imageid_mask); + + nc_close(ncid_mask); +/* */ +/* ----------------------- end netcdf calls for mask file --------------------------------------- */ + +/* ----------------------- netcdf calls for DEM file --------------------------------------- */ +/* */ + strcat(dem_nc_f.name,".nc"); + if( status = nc_open(dem_nc_f.name, 0, &ncid_dem) ) + { + printf("Error opening %s for netcdf read, exiting\n", dem_nc_f.name); + exit(1); + } + status = nc_inq_dimid (ncid_dem, "ncols", &colsid_dem); + status = nc_inq_dimlen(ncid_dem, colsid_dem, &lenp); + dem_nc_hdr.ncols = lenp; + status = nc_inq_dimid (ncid_dem, "nrows", &rowsid_dem); + status = nc_inq_dimlen(ncid_dem, rowsid_dem, &lenp); + dem_nc_hdr.nrows = lenp; + status = nc_get_att_text (ncid_dem, NC_GLOBAL, "tag", dem_nc_hdr.tag); + status = nc_get_att_double (ncid_dem, NC_GLOBAL, "ulx", &dem_nc_hdr.ulx); + status = nc_get_att_double (ncid_dem, NC_GLOBAL, "uly", &dem_nc_hdr.uly); + status = nc_get_att_double (ncid_dem, NC_GLOBAL, "cellsize", &dem_nc_hdr.cellsize); + status = nc_get_att_long (ncid_dem, NC_GLOBAL, "gctp_outsys", &dem_nc_hdr.gctp.outsys); + status = nc_get_att_long (ncid_dem, NC_GLOBAL, "gctp_outzone", &dem_nc_hdr.gctp.outzone); + status = nc_get_att_long (ncid_dem, NC_GLOBAL, "gctp_outdatum", &dem_nc_hdr.gctp.outdatum); + status = nc_get_att_double (ncid_dem, NC_GLOBAL, "gctp_outparm", dem_nc_hdr.gctp.outparm); + status = nc_get_att_text (ncid_dem, NC_GLOBAL, "gctp_fn27", dem_nc_hdr.gctp.fn27); + status = nc_get_att_text (ncid_dem, NC_GLOBAL, "gctp_fn83", dem_nc_hdr.gctp.fn83); + if (strncmp(dem_nc_hdr.tag,IMGTAG,7)) + { + printf("Error: dem file not a registered image. Use register_image first.\n"); + exit(1); + } + + status = nc_inq_varid (ncid_dem, "image", &imageid_dem); + + nc_close(ncid_dem); +/* */ +/* ----------------------- end netcdf calls for DEM file --------------------------------------- */ + +/* ----------------------- netcdf calls for station meta file --------------------------------------- */ +/* */ + if( status = nc_open(stnmeta_nc_f.name, NC_WRITE, &ncid_meta) ) + { + printf("Error opening %s for netcdf read, exiting\n", stnmeta_nc_f.name); + exit(1); + } + status = nc_inq_dimid (ncid_meta, "stations", &nstnid); + status = nc_inq_dimlen(ncid_meta, nstnid, &lenp); + meta_nc_hdr.nstns = lenp; + status = nc_inq_dimid (ncid_meta, "year_day", &yrdid); + status = nc_inq_dimlen(ncid_meta, yrdid, &lenp); + year_day = lenp; + status = nc_inq_dimid (ncid_meta, "descriptor_length", &descid); + status = nc_inq_dimlen(ncid_meta, descid, &lenp); + descriptor_length = lenp; + + status = nc_inq_varid (ncid_meta, "station_name", &stnnameid); + status = nc_inq_varid (ncid_meta, "station_id", &stnidid); + + status = nc_inq_varid (ncid_meta, "station_name", &stnnameid); + if (!(stnnames = (char*) malloc(meta_nc_hdr.nstns * descriptor_length * sizeof(char)))) + { + printf("Error allocating stnnames array.\n"); + exit(1); + } + status = nc_get_var_text (ncid_meta, stnnameid, stnnames); + + status = nc_inq_varid (ncid_meta, "station_id", &stnidid); + if (!(stnids = (char*) malloc(meta_nc_hdr.nstns * descriptor_length * sizeof(char)))) + { + printf("Error allocating stnids array.\n"); + exit(1); + } + status = nc_get_var_text (ncid_meta, stnidid, stnids); + + status = nc_inq_varid (ncid_meta, "in_tile", &in_tileid); + if (!(in_tile = (int*) malloc(meta_nc_hdr.nstns * sizeof(int)))) + { + printf("Error allocating in_tile array.\n"); + exit(1); + } + status = nc_get_var_int (ncid_meta, in_tileid, in_tile); + + status = nc_inq_varid (ncid_meta, "station_elevation", &stnelevid); + if (!(elevs = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) + { + printf("Error allocating elevs array. exiting.\n"); + exit(1); + } + status = nc_get_var_double (ncid_meta, stnelevid, elevs); + + status = nc_inq_varid (ncid_meta, "station_latitude", &stnlatid); + if (!(lats = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) + { + printf("Error allocating lats array. exiting.\n"); + exit(1); + } + status = nc_get_var_double (ncid_meta, stnlatid, lats); + + status = nc_inq_varid (ncid_meta, "station_longitude", &stnlonid); + if (!(lons = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) + { + printf("Error allocating lons array. exiting.\n"); + exit(1); + } + status = nc_get_var_double (ncid_meta, stnlonid, lons); + + status = nc_inq_varid (ncid_meta, "stnx", &stnxid); + if (!(stnx = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) + { + printf("Error allocating stnx array. exiting.\n"); + exit(1); + } + status = nc_get_var_double (ncid_meta, stnxid, stnx); + + status = nc_inq_varid (ncid_meta, "stny", &stnyid); + if (!(stny = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) + { + printf("Error allocating stny array. exiting.\n"); + exit(1); + } + status = nc_get_var_double (ncid_meta, stnyid, stny); + + status = nc_inq_varid (ncid_meta, "stnfx", &stnfxid); + if (!(stnfx = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) + { + printf("Error allocating stnfx array. exiting.\n"); + exit(1); + } + status = nc_get_var_double (ncid_meta, stnfxid, stnfx); + + status = nc_inq_varid (ncid_meta, "stnfy", &stnfyid); + if (!(stnfy = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) + { + printf("Error allocating stnfy array. exiting.\n"); + exit(1); + } + status = nc_get_var_double (ncid_meta, stnfyid, stnfy); + + status = nc_inq_varid (ncid_meta, "mflag", &mflagid); + if (!(mflag = (int*) malloc(meta_nc_hdr.nstns * year_day * sizeof(int)))) + { + printf("Error allocating mflag array.\n"); + exit(1); + } + status = nc_get_var_int (ncid_meta, mflagid, mflag); + + status = nc_inq_varid (ncid_meta, "precip", &prcpid); + if (!(prcp = (double*) malloc(meta_nc_hdr.nstns * year_day * sizeof(double)))) + { + printf("Error allocating prcp array. exiting.\n"); + exit(1); + } + status = nc_get_var_double (ncid_meta, prcpid, prcp); + + status = nc_get_att_int(ncid_meta, NC_GLOBAL, "xyloc_flag", &meta_nc_hdr.xyloc_flag); + + meta_nc_hdr.zfilt_flag = 0; +/* */ +/* ------------------- end netcdf calls for station meta file --------------------------------------- */ + + /* if z-filtering width other than 0.0, open the special ini file + containing filter type and mask and DEM filenames */ + if (fzwidth) + { + strcpy(ini.name, "zfilter.ini"); + if (file_open(&ini, 'i')) + { + printf("Error opening %s as special zfilter ini file\n",ini.name); + exit(1); + } + if (scan_value(ini, &type, 'i')) + { + printf("Error scanning type\n"); + exit(1); + } + if (scan_open(ini, &mask_f, 'r')) + { + printf("Error opening mask file\n"); + exit(1); + } + if (scan_open(ini, &dem_f, 'r')) + { + printf("Error opening dem file\n"); + exit(1); + } + fclose(ini.ptr); + } + + /* perform error checking on input file headers */ + if (fzwidth) + { + /* read mask file header and verify */ + if (strncmp(mask_nc_hdr.tag,IMGTAG,6)) + { + printf("Error: %s not a registered image. Use register_image first.\n",mask_f.name); + exit(1); + } + + /* read dem file header and verify */ + if (strncmp(dem_nc_hdr.tag,IMGTAG,6)) + { + printf("Error: %s not a registered image. Use register_image first.\n",dem_f.name); + exit(1); + } + + /* compare headers for DEM and mask */ + test = (mask_nc_hdr.ncols == dem_nc_hdr.ncols ) * (mask_nc_hdr.nrows == dem_nc_hdr.nrows) * + (mask_nc_hdr.ulx == dem_nc_hdr.ulx ) * (mask_nc_hdr.uly == dem_nc_hdr.uly ) * + (mask_nc_hdr.cellsize == dem_nc_hdr.cellsize); + if (!test) + { + printf("Error: Header information differs between mask and DEM...exiting\n"); + exit(1); + } + } + + /* set internal variables */ + if (fzwidth) + { + ncols = mask_nc_hdr.ncols; + nrows = mask_nc_hdr.nrows; + cellsize = dem_nc_hdr.cellsize; + ul_mapx = dem_nc_hdr.ulx; + ul_mapy = dem_nc_hdr.uly; + } + nstns = prcpstr.nstns = intstr.nstns = meta_nc_hdr.nstns; + ndays = prcpstr.ndays = year_day; + prcpstr.noz = 0; + /* the next two lines hardwire the xval period to the number of days + in the data files */ + start_day = 0; + nsimdays = ndays; + + if (fzwidth) + { + /* allocate zfilter memory */ + if (!(dem_array = (int*) malloc(nrows*ncols*sizeof(int)))) + { + printf("Error allocating dem array...exiting\n"); + exit(1); + } + if (!(mask_array = (long int*) malloc(nrows*ncols*sizeof(long int)))) + { + printf("Error allocating mask array...exiting\n"); + exit(1); + } + + /* read and close dem and mask files */ + fread(dem_array,sizeof(int),nrows*ncols,dem_f.ptr); + fread(mask_array,sizeof(long int),nrows*ncols,mask_f.ptr); + fclose(dem_f.ptr); + fclose(mask_f.ptr); + + /* set size of elevation filter window */ + width = (int) (fzwidth / cellsize); + if (!(width % 2)) width++; + hw = width/2; + + /* generate filter kernel */ + if (!(kernel = (double*) malloc(width * width * sizeof(double)))) + { + printf("Error allocating for kernel... exiting\n"); + exit(1); + } + if (make_kernel(kernel, width, type)) + { + printf("Error in call to make_kernel() ... exiting\n"); + exit(1); + } + } + + /* allocate memory */ + if (!(good = (char*) malloc(nstns * ndays * sizeof(char)))) + { + printf("Error allocating good array...exiting\n"); + exit(1); + } + if (!(prcpstr.stnobs = (double*) malloc(nstns * ndays * sizeof(double)))) + { + printf("Error allocating stnobs array...exiting\n"); + exit(1); + } + if (!(prcpstr.stnsmobs = (double*) malloc(nstns * ndays * sizeof(double)))) + { + printf("Error allocating stnmobs array...exiting\n"); + exit(1); + } + if (!(prcpstr.stnx = (double*) malloc(nstns * sizeof(double)))) + { + printf("Error allocating stnx array...exiting\n"); + exit(1); + } + if (!(prcpstr.stny = (double*) malloc(nstns * sizeof(double)))) + { + printf("Error allocating stny array...exiting\n"); + exit(1); + } + if (!(prcpstr.stnz = (double*) malloc(nstns * sizeof(double)))) + { + printf("Error allocating stnx array...exiting\n"); + exit(1); + } + if (!(intstr.stnx = (double*) malloc(nstns * sizeof(double)))) + { + printf("Error allocating stnx array...exiting\n"); + exit(1); + } + if (!(intstr.stny = (double*) malloc(nstns * sizeof(double)))) + { + printf("Error allocating stny array...exiting\n"); + exit(1); + } + if (!(intstr.sqdist = (double*) malloc(nstns * sizeof(double)))) + { + printf("Error allocating sqdist array...exiting\n"); + exit(1); + } + if (!(intstr.wt = (double*) malloc(nstns * sizeof(double)))) + { + printf("Error allocating wt array...exiting\n"); + exit(1); + } + if (!(intstr.id = (short*) malloc(nstns * sizeof(short)))) + { + printf("Error allocating id array...exiting\n"); + exit(1); + } + /* index array needed for bubble sort algorithm */ + if (!(bsi = (short*) malloc(nstns * sizeof(short)))) {printf("Couild not malloc bsi\n"); exit(1);} + /* DEBUG SNOTEL CODE ADDED 8/29/00 by PET */ + if (DEBUG_SNOTEL) + { + if (!(debug_metastr = (stnmeta_struct*) malloc(nstns * sizeof(stnmeta_struct)))) + { + printf("Error allocating debug_metastr...exiting\n"); + exit(1); + } + } + + /* istr and tstr share space for wt and id arrays */ + prcpstr.wt = intstr.wt; + prcpstr.id = intstr.id; + if (inflag) + { + if (!(in = (int*) malloc(nstns * sizeof(int)))) + { + printf("Error allocating in array...exiting\n"); + exit(1); + } + } + if (!(snot = (int*) malloc(nstns * sizeof(int)))) + { + printf("Error allocating snot array...exiting\n"); + exit(1); + } + + /* allocate memory for output arrays */ + if (!(ngooddays = (double*) malloc(nstns * sizeof(double)))) + { + printf("Error allocating ngooddays array...exiting\n"); + exit(1); + } + if (!(ngoodstns = (double*) malloc(ndays * sizeof(double)))) + { + printf("Error allocating ngoodstns array...exiting\n"); + exit(1); + } + if (!(xlrdaysum = (double*) malloc(ndays * sizeof(double)))) + { + printf("Error allocating xlrdaysum array...exiting\n"); + exit(1); + } + if (!(xlrdaysumsq = (double*) malloc(ndays * sizeof(double)))) + { + printf("Error allocating xlrdaysumsq array...exiting\n"); + exit(1); + } + if (!(ylrdaysum = (double*) malloc(ndays * sizeof(double)))) + { + printf("Error allocating ylrdaysum array...exiting\n"); + exit(1); + } + if (!(ylrdaysumsq = (double*) malloc(ndays * sizeof(double)))) + { + printf("Error allocating ylrdaysumsq array...exiting\n"); + exit(1); + } + if (!(zlrdaysum = (double*) malloc(ndays * sizeof(double)))) + { + printf("Error allocating zlrdaysum array...exiting\n"); + exit(1); + } + if (!(zlrdaysumsq = (double*) malloc(ndays * sizeof(double)))) + { + printf("Error allocating zlrdaysumsq array...exiting\n"); + exit(1); + } + + /* read stnmeta data and fill appropriate arrays */ + nstns_nws = nstns_snot = 0; + nstns_intile = 0; + for (i=0 ; i= ncols) || + (stn_filey < 0) || (stn_filey >= nrows) || + (!mask_array[stn_offset])) + { + prcpstr.stnz[i] = metastr.z = elevs[i]; + continue; + } + + /* station is in mask, and inside DEM data region */ + /* begin loops through kernel to calculate filtered elevation */ + value = 0.0; + sum_wt = 0.0; + for (k=0 ; k= nrows)) + continue; + + koff = k * width; + img_off = y * ncols; + + for (l=0 ; l= ncols)) + continue; + + koffset = koff + l; + img_offset = img_off + x; + + /* finally, check mask for convolution overlap */ + if (!mask_array[img_offset]) + continue; + + /* accumulate convolved values */ + value += dem_array[img_offset] * kernel[koffset]; + sum_wt += kernel[koffset]; + } + } + + /* check that at least some good values were found */ + if (!sum_wt) + { + prcpstr.stnz[i] = metastr.z = elevs[i]; + continue; + } + prcpstr.stnz[i] = value/sum_wt; + } + else /* fzwidth == 0 */ + { + prcpstr.stnz[i] = metastr.z = elevs[i]; + } + + if (inflag) + { + in[i] = in_tile[i]; + if (in[i] == inval) nstns_intile++; + } + + } /* end of nstns loop */ + + if (nstns_intile) + { + /* allocate space for xval ouput */ + if (!(out_stnnames = (char*)malloc(nstns_intile * descriptor_length * sizeof(char)))) + { + printf("Error allocating for out_stnnames array.\n"); + exit(1); + } + if (!(out_stnids = (char*)malloc(nstns_intile * descriptor_length * sizeof(char)))) + { + printf("Error allocating for out_stnids array.\n"); + exit(1); + } + if (!(out_stnx = (double*)malloc(nstns_intile * sizeof(double)))) + { + printf("Error allocating for out_stnx array... exiting\n"); + exit(1); + } + if (!(out_stny = (double*)malloc(nstns_intile * sizeof(double)))) + { + printf("Error allocating for out_stny array... exiting\n"); + exit(1); + } + if (!(out_stnz = (double*)malloc(nstns_intile * sizeof(double)))) + { + printf("Error allocating for out_stnz array... exiting\n"); + exit(1); + } + if (!(out_obs = (double*)malloc(nstns_intile*ndays* sizeof(double)))) + { + printf("Error allocating for out_obs array... exiting\n"); + exit(1); + } + if (!(out_pred = (double*)malloc(nstns_intile*ndays* sizeof(double)))) + { + printf("Error allocating for out_pred array... exiting\n"); + exit(1); + } + if (!(out_good = (int*)malloc(nstns_intile*ndays* sizeof(int)))) + { + printf("Error allocating for out_good array... exiting\n"); + exit(1); + } + + /* netcdf define dimensions and array variables for output arrays */ + len_stns = nstns_intile; + len_days = ndays; + nc_def_dim(ncid_out, "stns", len_stns, &outdid_stns); + nc_def_dim(ncid_out, "descriptor_length", descriptor_length, &outid_desc); + nc_def_dim(ncid_out, "days", len_days, &outdid_days); + + stndims[0] = outdid_stns; + stndims[1] = outid_desc; + out_did[0] = outdid_stns; + out_did[1] = outdid_days; + out_did1[0] = outdid_stns; + + nc_def_var (ncid_out, "station_name", NC_CHAR, 2, stndims, &outid_stnname); + nc_def_var (ncid_out, "station_id", NC_CHAR, 2, stndims, &outid_stnid); + + nc_def_var(ncid_out, "stnx", NC_DOUBLE, 1, out_did1, &outid_stnx); + strcpy(att,"station projected x coordinate"); + nc_put_att_text(ncid_out, outid_stnx, "longname", strlen(att), att); + strcpy(att,"m"); + nc_put_att_text(ncid_out, outid_stnx, "units", strlen(att), att); + + nc_def_var(ncid_out, "stny", NC_DOUBLE, 1, out_did1, &outid_stny); + strcpy(att,"station projected y coordinate"); + nc_put_att_text(ncid_out, outid_stny, "longname", strlen(att), att); + strcpy(att,"m"); + nc_put_att_text(ncid_out, outid_stny, "units", strlen(att), att); + + nc_def_var(ncid_out, "stnz", NC_DOUBLE, 1, out_did1, &outid_stnz); + strcpy(att,"station elevation"); + nc_put_att_text(ncid_out, outid_stnz, "longname", strlen(att), att); + strcpy(att,"m"); + nc_put_att_text(ncid_out, outid_stnz, "units", strlen(att), att); + + nc_def_var(ncid_out, "obs", NC_DOUBLE, 2, out_did, &outid_obs); + strcpy(att,"observed precipitation"); + nc_put_att_text(ncid_out, outid_obs, "longname", strlen(att), att); + strcpy(att,"cm/day"); + nc_put_att_text(ncid_out, outid_obs, "units", strlen(att), att); + + nc_def_var(ncid_out, "pred", NC_DOUBLE, 2, out_did, &outid_pred); + strcpy(att,"predicted precipitation"); + nc_put_att_text(ncid_out, outid_pred, "longname", strlen(att), att); + strcpy(att,"cm/day"); + nc_put_att_text(ncid_out, outid_pred, "units", strlen(att), att); + + nc_def_var(ncid_out, "good", NC_INT, 2, out_did, &outid_good); + strcpy(att,"good value flag (0=missing obs, 1=good obs)"); + nc_put_att_text(ncid_out, outid_good, "longname", strlen(att), att); + strcpy(att,"flag"); + nc_put_att_text(ncid_out, outid_good, "units", strlen(att), att); + } /* end if nstns_intile */ + + /* take output file out of define mode, ready for write */ + nc_enddef(ncid_out); + + /* read station daily data and fill appropriate arrays */ + for (i=0 ; i 49) i = 49; + frad = f90[i]; + + /* calculate initial density filter parameters */ + intstr.dfsqrad = dfirad * dfirad; + intstr.inv_dfsqrad = 1.0 / intstr.dfsqrad; + intstr.dfarea = PI * intstr.dfsqrad; + intstr.trunc = exp(-intstr.gsp); + intstr.dftrunc = exp(-DFGSP); + intstr.dfavgwt = ((1.0 - intstr.dftrunc)/DFGSP) - intstr.dftrunc; + + /* start loop through stations */ + mae1 = mpae1 = bias1 = 0.0; + mae2 = mpae2 = bias2 = 0.0; + prcpstr.switch_init = 1; + n_in = 0; + nper = 0; + ngood = ngoodper = 0.0; + nlr = 0.0; + ngoodlr = 0; + for (i=0 ; ij ; k--) { + if (intstr.sqdist[bsi[k]] < intstr.sqdist[bsi[k-1]]) { + temp_short = bsi[k]; + bsi[k]=bsi[k-1]; + bsi[k-1]=temp_short; + } + } + } + + /* set the weighting kernel radius to 1m greater than the furthest included station */ + /* the extra 1.0m ensures that none of the stations in the list will have wt=0 */ + intstr.sqrad=intstr.sqdist[bsi[ans]]+1.0; + + /* fill the list of included stations, being sure to leave out the prediction point + station */ + k=0; + wtsum=0.0; + for (j=0 ; j 0.0)); + nwetstns += wet; + } + } + if (sumpop) pop /= sumpop; + else + { + printf("no good data for the day\n"); + continue; + } + + /* if precipitation occurrence threshold exceeded, predict + precipitation amount */ + if (pop > pop_crit) + { + + /* only generate regression stats when nwetstns > 3 */ + if (nwetstns > 3) + { + /* generate regression arrays that change by day */ + if (xv_prcp_regr_y_2switch(listgood, &prcpstr)) + { + printf("Error in xv_prcp_regr_y_switch() ... exiting\n"); + exit(1); + } + + /* weighted multiple least squares fit */ + if (prcp_mlr(&prcpstr)) + { + printf("Error in wt_regr() ... exiting\n"); + exit(1); + } + + } /* end if nwet > 3 */ + + else + { + prcpstr.coef[0] = prcpstr.coef[1] = prcpstr.coef[2] = + prcpstr.coef[3] = 0.0; + } + + /* update mean and stdv for slope (lapse rate) */ + if (goodlr) + { + ngoodstns[day] ++; + xlrdaysum[day] += prcpstr.coef[0]; + xlrdaysumsq[day] += prcpstr.coef[0] * prcpstr.coef[0]; + ylrdaysum[day] += prcpstr.coef[1]; + ylrdaysumsq[day] += prcpstr.coef[1] * prcpstr.coef[1]; + zlrdaysum[day] += prcpstr.coef[2]; + zlrdaysumsq[day] += prcpstr.coef[2] * prcpstr.coef[2]; + } + + /* predict prcp for the day */ + if (xv_predict_prcp_noint(listgood, &prcpstr)) + { + printf("Error in xv_predict_prcp_noint() ... exiting\n"); + exit(1); + } + + } /* end if pop > pop_crit */ + + /* calculate the error and bias for this station-day */ + pred = prcpstr.pp; + obs = prcpstr.stnobs[obs_off+day]; + + out_obs[(n_in-1)*nsimdays+day]=obs; + out_pred[(n_in-1)*nsimdays+day]=pred; + out_good[(n_in-1)*nsimdays+day]=1; + + b = pred - obs; + ae = fabs(b); + + /* this day has good observation and good predictor data */ + ngooddays[i]++; + ngood++; + + bias1 += b; + mae1 += ae; + + pred_tot += pred; + obs_tot += obs; + + } /* end of nsimdays loop */ + + /* free memory depending on count */ + free(listgood); + free(prcpstr.listobs); + free(prcpstr.listsmobs); + free(prcpstr.listx); + free(prcpstr.listy); + free(prcpstr.listz); + for (j=0 ; j<3 ; j++) + { + free(prcpstr.regx[j]); + } + free(prcpstr.regx); + free(prcpstr.regy); + free(prcpstr.regwtall); + free(prcpstr.regwt); + + /* calculate period-of-record averages */ + if (ngooddays[i]) + { + /* DEBUG SNOTEL CODE ADDED 8/29/00 by PET */ + if (DEBUG_SNOTEL) + { + fprintf(debug_snotel_f.ptr,"%3d%10.4lf%10.4lf%10.1lf%10.1lf%10.1lf%10s %s\n", + debug_metastr[i].type,debug_metastr[i].lon,debug_metastr[i].lat, + debug_metastr[i].z,obs_tot,pred_tot,debug_metastr[i].code,debug_metastr[i].name); + } + rn = ngooddays[i]; + /* calculate the errors over the entire ndays period, expressed + as an average daily error */ + pred_tot /= rn; + obs_tot /= rn; + b = pred_tot - obs_tot; + ae = fabs(b); + /* only do the percent calculation if there is observed precip + for the period, to avoid division by zero error */ + if (obs_tot) + { + pae = 100.0 * ae / obs_tot; + /* weight the multi-station average by the number of good + days at each station */ + mpae2 += rn * pae; + ngoodper += rn; + } + mae2 += rn * ae; + bias2 += rn * b; + } + + /* switch initial sign for next station */ + prcpstr.switch_init = -prcpstr.switch_init; + + } /* end of nstations loop */ + + if (nstns_intile) + { + /* write netcdf output arrays */ + nc_put_var_text(ncid_out, outid_stnname, out_stnnames); + nc_put_var_text(ncid_out, outid_stnid, out_stnids); + nc_put_var_double(ncid_out, outid_stnx, out_stnx); + nc_put_var_double(ncid_out, outid_stny, out_stny); + nc_put_var_double(ncid_out, outid_stnz, out_stnz); + nc_put_var_int(ncid_out, outid_good, out_good); + nc_put_var_double(ncid_out, outid_obs, out_obs); + nc_put_var_double(ncid_out, outid_pred, out_pred); + } + + /* calculate final xval stats */ + xlrstdv = xlrmean = 0.0; + ylrstdv = ylrmean = 0.0; + zlrstdv = zlrmean = 0.0; + if (ngood) + { + rn = ngood; + mae1 /= rn; + bias1 /= rn; + mae2 /= rn; + bias2 /= rn; + + if (ngoodper) mpae2 /= ngoodper; + + /* calculate means and stdevs for lr */ + for (day=0 ; day \n"); + exit(1); + } + + /* scan command line arguments into parameters */ + job = atoi(argv[1]); + strcpy(mask_nc_f.name, argv[2]); + strcpy(dem_nc_f.name, argv[3]); + strcpy(stnmeta_nc_f.name, argv[4]); + strcpy(stndata_nc_f.name, argv[5]); + fzwidth = atof(argv[6]); + dfirad = atof(argv[7]); + init_gsp = atof(argv[8]); + intstr.ans = atof(argv[9]); + smwidth = atoi(argv[10]); + inflag = atoi(argv[11]); + inval = atoi(argv[12]); + strcpy(outprefix,argv[13]); + ans=(int)intstr.ans; + + /* create netcdf xval output file */ + strcpy(out_f.name, outprefix); + strcat(out_f.name,"_xval.nc"); + if (status = nc_create(out_f.name, 0, &ncid_out)) + { + printf("Error creating %s as netcdf file, exiting\n", out_f.name); + exit(1); + } + /* define scalar variables for xval output file */ + if (status = nc_def_var(ncid_out, "tileid", NC_INT, 0, 0, &outid_tileid)) + { + printf("Error creating tileid variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"Daymet tile ID"); + nc_put_att_text(ncid_out, outid_tileid, "longname", strlen(att), att); + strcpy(att,"index"); + nc_put_att_text(ncid_out, outid_tileid, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "nstns3x3", NC_INT, 0, 0, &outid_nstns3x3)) + { + printf("Error creating nstns3x3 variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"number of stations in 3x3 tile region"); + nc_put_att_text(ncid_out, outid_nstns3x3, "longname", strlen(att), att); + strcpy(att,"stations"); + nc_put_att_text(ncid_out, outid_nstns3x3, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "nstns", NC_INT, 0, 0, &outid_nstns)) + { + printf("Error creating nstns variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"number of stations evaluated (central tile)"); + nc_put_att_text(ncid_out, outid_nstns, "longname", strlen(att), att); + strcpy(att,"stations"); + nc_put_att_text(ncid_out, outid_nstns, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "nstndays", NC_INT, 0, 0, &outid_nstndays)) + { + printf("Error creating nstndays variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"number of station-days evaluated"); + nc_put_att_text(ncid_out, outid_nstndays, "longname", strlen(att), att); + strcpy(att,"days"); + nc_put_att_text(ncid_out, outid_nstndays, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "rad90mean", NC_DOUBLE, 0, 0, &outid_rad90mean)) + { + printf("Error creating rad90mean variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"mean: radius capturing 90% of filter kernel weight"); + nc_put_att_text(ncid_out, outid_rad90mean, "longname", strlen(att), att); + strcpy(att,"m"); + nc_put_att_text(ncid_out, outid_rad90mean, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "rad90stdv", NC_DOUBLE, 0, 0, &outid_rad90stdv)) + { + printf("Error creating rad90stdv variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"std dev: radius capturing 90% of filter kernel weight"); + nc_put_att_text(ncid_out, outid_rad90stdv, "longname", strlen(att), att); + strcpy(att,"m"); + nc_put_att_text(ncid_out, outid_rad90stdv, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "daymae", NC_DOUBLE, 0, 0, &outid_daymae)) + { + printf("Error creating daymae variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"mean absolute error for single-day predictions"); + nc_put_att_text(ncid_out, outid_daymae, "longname", strlen(att), att); + strcpy(att,"degrees C"); + nc_put_att_text(ncid_out, outid_daymae, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "pormae", NC_DOUBLE, 0, 0, &outid_pormae)) + { + printf("Error creating annmae variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"mean absolute error for period-of-record predictions"); + nc_put_att_text(ncid_out, outid_pormae, "longname", strlen(att), att); + strcpy(att,"degrees C"); + nc_put_att_text(ncid_out, outid_pormae, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "bias", NC_DOUBLE, 0, 0, &outid_bias)) + { + printf("Error creating bias variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"mean prediction bias"); + nc_put_att_text(ncid_out, outid_bias, "longname", strlen(att), att); + strcpy(att,"degrees C"); + nc_put_att_text(ncid_out, outid_bias, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "tamean", NC_DOUBLE, 0, 0, &outid_tamean)) + { + printf("Error creating ppmean variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"mean observed temperature"); + nc_put_att_text(ncid_out, outid_tamean, "longname", strlen(att), att); + strcpy(att,"degrees C"); + nc_put_att_text(ncid_out, outid_tamean, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "xlrmean", NC_DOUBLE, 0, 0, &outid_xlrmean)) + { + printf("Error creating xlrmean variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"3-d regression: mean x-component"); + nc_put_att_text(ncid_out, outid_xlrmean, "longname", strlen(att), att); + strcpy(att,"degrees C/m"); + nc_put_att_text(ncid_out, outid_xlrmean, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "xlrstdv", NC_DOUBLE, 0, 0, &outid_xlrstdv)) + { + printf("Error creating xlrstdv variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"3-d regression: among-station std dev of x-component"); + nc_put_att_text(ncid_out, outid_xlrstdv, "longname", strlen(att), att); + strcpy(att,"degrees C/m"); + nc_put_att_text(ncid_out, outid_xlrstdv, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "ylrmean", NC_DOUBLE, 0, 0, &outid_ylrmean)) + { + printf("Error creating ylrmean variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"3-d regression: mean y-component"); + nc_put_att_text(ncid_out, outid_ylrmean, "longname", strlen(att), att); + strcpy(att,"degrees C/m"); + nc_put_att_text(ncid_out, outid_ylrmean, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "ylrstdv", NC_DOUBLE, 0, 0, &outid_ylrstdv)) + { + printf("Error creating ylrstdv variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"3-d regression: among-station std dev of y-component"); + nc_put_att_text(ncid_out, outid_ylrstdv, "longname", strlen(att), att); + strcpy(att,"degrees C/m"); + nc_put_att_text(ncid_out, outid_ylrstdv, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "zlrmean", NC_DOUBLE, 0, 0, &outid_zlrmean)) + { + printf("Error creating zlrmean variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"3-d regression: mean z-component"); + nc_put_att_text(ncid_out, outid_zlrmean, "longname", strlen(att), att); + strcpy(att,"derees C/m"); + nc_put_att_text(ncid_out, outid_zlrmean, "units", strlen(att), att); + + if (status = nc_def_var(ncid_out, "zlrstdv", NC_DOUBLE, 0, 0, &outid_zlrstdv)) + { + printf("Error creating zlrstdv variable for %s, exiting\n", out_f.name); + exit(1); + } + strcpy(att,"3-d regression: among-station std dev of z-component"); + nc_put_att_text(ncid_out, outid_zlrstdv, "longname", strlen(att), att); + strcpy(att,"degrees C/m"); + nc_put_att_text(ncid_out, outid_zlrstdv, "units", strlen(att), att); + + +/* ----------------------- netcdf calls for mask file --------------------------------------- */ +/* */ + strcat(mask_nc_f.name,".nc"); + if( status = nc_open(mask_nc_f.name, 0, &ncid_mask) ) + { + printf("Error opening %s for netcdf read, exiting\n", mask_nc_f.name); + exit(1); + } + status = nc_inq_dimid (ncid_mask, "ncols", &colsid_mask); + status = nc_inq_dimlen(ncid_mask, colsid_mask, &lenp); + mask_nc_hdr.ncols = lenp; + status = nc_inq_dimid (ncid_mask, "nrows", &rowsid_mask); + status = nc_inq_dimlen(ncid_mask, rowsid_mask, &lenp); + mask_nc_hdr.nrows = lenp; + status = nc_get_att_text (ncid_mask, NC_GLOBAL, "tag", mask_nc_hdr.tag); + status = nc_get_att_double (ncid_mask, NC_GLOBAL, "ulx", &mask_nc_hdr.ulx); + status = nc_get_att_double (ncid_mask, NC_GLOBAL, "uly", &mask_nc_hdr.uly); + status = nc_get_att_double (ncid_mask, NC_GLOBAL, "cellsize", &mask_nc_hdr.cellsize); + status = nc_get_att_long (ncid_mask, NC_GLOBAL, "gctp_outsys", &mask_nc_hdr.gctp.outsys); + status = nc_get_att_long (ncid_mask, NC_GLOBAL, "gctp_outzone", &mask_nc_hdr.gctp.outzone); + status = nc_get_att_long (ncid_mask, NC_GLOBAL, "gctp_outdatum", &mask_nc_hdr.gctp.outdatum); + status = nc_get_att_double (ncid_mask, NC_GLOBAL, "gctp_outparm", mask_nc_hdr.gctp.outparm); + status = nc_get_att_text (ncid_mask, NC_GLOBAL, "gctp_fn27", mask_nc_hdr.gctp.fn27); + status = nc_get_att_text (ncid_mask, NC_GLOBAL, "gctp_fn83", mask_nc_hdr.gctp.fn83); + if (strncmp(mask_nc_hdr.tag,IMGTAG,7)) + { + printf("Error: mask file not a registered image. Use register_image first.\n"); + exit(1); + } + + status = nc_inq_varid (ncid_mask, "image", &imageid_mask); + + nc_close(ncid_mask); +/* */ +/* ----------------------- end netcdf calls for mask file --------------------------------------- */ + +/* ----------------------- netcdf calls for DEM file --------------------------------------- */ +/* */ + strcat(dem_nc_f.name,".nc"); + if( status = nc_open(dem_nc_f.name, 0, &ncid_dem) ) + { + printf("Error opening %s for netcdf read, exiting\n", dem_nc_f.name); + exit(1); + } + status = nc_inq_dimid (ncid_dem, "ncols", &colsid_dem); + status = nc_inq_dimlen(ncid_dem, colsid_dem, &lenp); + dem_nc_hdr.ncols = lenp; + status = nc_inq_dimid (ncid_dem, "nrows", &rowsid_dem); + status = nc_inq_dimlen(ncid_dem, rowsid_dem, &lenp); + dem_nc_hdr.nrows = lenp; + status = nc_get_att_text (ncid_dem, NC_GLOBAL, "tag", dem_nc_hdr.tag); + status = nc_get_att_double (ncid_dem, NC_GLOBAL, "ulx", &dem_nc_hdr.ulx); + status = nc_get_att_double (ncid_dem, NC_GLOBAL, "uly", &dem_nc_hdr.uly); + status = nc_get_att_double (ncid_dem, NC_GLOBAL, "cellsize", &dem_nc_hdr.cellsize); + status = nc_get_att_long (ncid_dem, NC_GLOBAL, "gctp_outsys", &dem_nc_hdr.gctp.outsys); + status = nc_get_att_long (ncid_dem, NC_GLOBAL, "gctp_outzone", &dem_nc_hdr.gctp.outzone); + status = nc_get_att_long (ncid_dem, NC_GLOBAL, "gctp_outdatum", &dem_nc_hdr.gctp.outdatum); + status = nc_get_att_double (ncid_dem, NC_GLOBAL, "gctp_outparm", dem_nc_hdr.gctp.outparm); + status = nc_get_att_text (ncid_dem, NC_GLOBAL, "gctp_fn27", dem_nc_hdr.gctp.fn27); + status = nc_get_att_text (ncid_dem, NC_GLOBAL, "gctp_fn83", dem_nc_hdr.gctp.fn83); + if (strncmp(dem_nc_hdr.tag,IMGTAG,7)) + { + printf("dem file not a registered image. Use register_image first.\n"); + exit(1); + } + + status = nc_inq_varid (ncid_dem, "image", &imageid_dem); + + nc_close(ncid_dem); +/* */ +/* ----------------------- end netcdf calls for DEM file --------------------------------------- */ + +/* ----------------------- netcdf calls for station meta file --------------------------------------- */ +/* */ + if( status = nc_open(stnmeta_nc_f.name, NC_WRITE, &ncid_meta) ) + { + printf("Error opening %s for netcdf read, exiting\n", stnmeta_nc_f.name); + exit(1); + } + status = nc_inq_dimid (ncid_meta, "stations", &nstnid); + status = nc_inq_dimlen(ncid_meta, nstnid, &lenp); + meta_nc_hdr.nstns = lenp; + status = nc_inq_dimid (ncid_meta, "year_day", &yrdid); + status = nc_inq_dimlen(ncid_meta, yrdid, &lenp); + year_day = lenp; + status = nc_inq_dimid (ncid_meta, "descriptor_length", &descid); + status = nc_inq_dimlen(ncid_meta, descid, &lenp); + descriptor_length = lenp; + + status = nc_inq_varid (ncid_meta, "station_name", &stnnameid); + status = nc_inq_varid (ncid_meta, "station_id", &stnidid); + + status = nc_inq_varid (ncid_meta, "station_name", &stnnameid); + if (!(stnnames = (char*) malloc(meta_nc_hdr.nstns * descriptor_length * sizeof(char)))) + { + printf("Error allocating stnnames array.\n"); + exit(1); + } + status = nc_get_var_text (ncid_meta, stnnameid, stnnames); + + status = nc_inq_varid (ncid_meta, "station_id", &stnidid); + if (!(stnids = (char*) malloc(meta_nc_hdr.nstns * descriptor_length * sizeof(char)))) + { + printf("Error allocating stnids array.\n"); + exit(1); + } + status = nc_get_var_text (ncid_meta, stnidid, stnids); + + status = nc_inq_varid (ncid_meta, "in_tile", &in_tileid); + if (!(in_tile = (int*) malloc(meta_nc_hdr.nstns * sizeof(int)))) + { + printf("Error allocating in_tile array.\n"); + exit(1); + } + status = nc_get_var_int (ncid_meta, in_tileid, in_tile); + + status = nc_inq_varid (ncid_meta, "station_elevation", &stnelevid); + if (!(elevs = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) + { + printf("Error allocating elevs array.\n"); + exit(1); + } + status = nc_get_var_double (ncid_meta, stnelevid, elevs); + + status = nc_inq_varid (ncid_meta, "station_latitude", &stnlatid); + if (!(lats = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) + { + printf("Error allocating lats array.\n"); + exit(1); + } + status = nc_get_var_double (ncid_meta, stnlatid, lats); + + status = nc_inq_varid (ncid_meta, "station_longitude", &stnlonid); + if (!(lons = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) + { + printf("Error allocating lons array.\n"); + exit(1); + } + status = nc_get_var_double (ncid_meta, stnlonid, lons); + + status = nc_inq_varid (ncid_meta, "stnx", &stnxid); + if (!(stnx = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) + { + printf("Error allocating stnx array.\n"); + exit(1); + } + status = nc_get_var_double (ncid_meta, stnxid, stnx); + + status = nc_inq_varid (ncid_meta, "stny", &stnyid); + if (!(stny = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) + { + printf("Error allocating stny array.\n"); + exit(1); + } + status = nc_get_var_double (ncid_meta, stnyid, stny); + + status = nc_inq_varid (ncid_meta, "stnfx", &stnfxid); + if (!(stnfx = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) + { + printf("Error allocating stnfx array.\n"); + exit(1); + } + status = nc_get_var_double (ncid_meta, stnfxid, stnfx); + + status = nc_inq_varid (ncid_meta, "stnfy", &stnfyid); + if (!(stnfy = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) + { + printf("Error allocating stnfy array.\n"); + exit(1); + } + status = nc_get_var_double (ncid_meta, stnfyid, stnfy); + + status = nc_inq_varid (ncid_meta, "mflag", &mflagid); + if (!(mflag = (int*) malloc(meta_nc_hdr.nstns * year_day * sizeof(int)))) + { + printf("Error allocating mflag array.\n"); + exit(1); + } + status = nc_get_var_int (ncid_meta, mflagid, mflag); + + status = nc_inq_varid (ncid_meta, "code", &codeid); + if (!(code = (int*) malloc(meta_nc_hdr.nstns * year_day * sizeof(int)))) + { + printf("Error allocating code array.\n"); + exit(1); + } + status = nc_get_var_int (ncid_meta, codeid, code); + + status = nc_inq_varid (ncid_meta, "tmax", &tairid); + if (status) + { + status = nc_inq_varid (ncid_meta, "tmin", &tairid); + } + if (status) + { + status = nc_inq_varid (ncid_meta, "tair", &tairid); + } + if (!(tair = (double*) malloc(meta_nc_hdr.nstns * year_day * sizeof(double)))) + { + printf("Error allocating tair array.\n"); + exit(1); + } + status = nc_get_var_double (ncid_meta, tairid, tair); + + + status = nc_get_att_int(ncid_meta, NC_GLOBAL, "xyloc_flag", &meta_nc_hdr.xyloc_flag); +/* */ +/* ------------------- end netcdf calls for station meta file --------------------------------------- */ + + /* if z-filtering width other than 0.0, open the special ini file + containing filter type and mask and DEM filenames */ + if (fzwidth) + { + strcpy(ini.name, "zfilter.ini"); + if (file_open(&ini, 'i')) + { + printf("Error opening %s as special zfilter ini file\n",ini.name); + exit(1); + } + if (scan_value(ini, &type, 'i')) + { + printf("Error scanning type\n"); + exit(1); + } + if (scan_open(ini, &mask_f, 'r')) + { + printf("Error opening mask file\n"); + exit(1); + } + if (scan_open(ini, &dem_f, 'r')) + { + printf("Error opening dem file\n"); + exit(1); + } + fclose(ini.ptr); + } + + /* perform error checking on input file headers */ + if (fzwidth) + { + /* read mask file header and verify */ + if (strncmp(mask_nc_hdr.tag,IMGTAG,6)) + { + printf("Error: %s not a registered image. Use register_image first.\n",mask_f.name); + exit(1); + } + + /* read dem file header and verify */ + if (strncmp(dem_nc_hdr.tag,IMGTAG,6)) + { + printf("Error: %s not a registered image. Use register_image first.\n",dem_f.name); + exit(1); + } + + /* compare headers for DEM and mask */ + test = (mask_nc_hdr.ncols == dem_nc_hdr.ncols ) * (mask_nc_hdr.nrows == dem_nc_hdr.nrows) * + (mask_nc_hdr.ulx == dem_nc_hdr.ulx ) * (mask_nc_hdr.uly == dem_nc_hdr.uly ) * + (mask_nc_hdr.cellsize == dem_nc_hdr.cellsize); + if (!test) + { + printf("Error: Header information differs between mask and DEM...exiting\n"); + exit(1); + } + } + + /* set internal variables */ + if (fzwidth) + { + ncols = mask_nc_hdr.ncols; + nrows = mask_nc_hdr.nrows; + cellsize = dem_nc_hdr.cellsize; + ul_mapx = dem_nc_hdr.ulx; + ul_mapy = dem_nc_hdr.uly; + } + nstns = tairstr.nstns = intstr.nstns = meta_nc_hdr.nstns; + ndays = tairstr.ndays = year_day; + tairstr.noz = 0; + /* the next two lines hardwire the xval period to the number of days + in the data files */ + start_day = 0; + nsimdays = ndays; + + if (fzwidth) + { + /* allocate zfilter memory */ + if (!(dem_array = (int*) malloc(nrows*ncols*sizeof(int)))) + { + printf("Error allocating dem array.\n"); + exit(1); + } + if (!(mask_array = (long int*) malloc(nrows*ncols*sizeof(long int)))) + { + printf("Error allocating mask array.\n"); + exit(1); + } + + /* read and close dem and mask files */ + fread(dem_array,sizeof(int),nrows*ncols,dem_f.ptr); + fread(mask_array,sizeof(long int),nrows*ncols,mask_f.ptr); + fclose(dem_f.ptr); + fclose(mask_f.ptr); + + /* set size of elevation filter window */ + width = (int) (fzwidth / cellsize); + if (!(width % 2)) width++; + hw = width/2; + + /* generate filter kernel */ + if (!(kernel = (double*) malloc(width * width * sizeof(double)))) + { + printf("Error allocating for kernel... exiting\n"); + exit(1); + } + if (make_kernel(kernel, width, type)) + { + printf("Error in call to make_kernel() ... exiting\n"); + exit(1); + } + } + + /* allocate memory */ + if (!(good = (char*) malloc(nstns * ndays * sizeof(char)))) + { + printf("Error allocating good array... exiting\n"); + exit(1); + } + if (!(tairstr.stnobs = (double*) malloc(nstns * ndays * sizeof(double)))) + { + printf("Error allocating stnsobs array... exiting\n"); + exit(1); + } + if (!(tairstr.stnsmobs = (double*) malloc(nstns * ndays * sizeof(double)))) + { + printf("Error allocating stnsmobs array... exiting\n"); + exit(1); + } + if (!(tairstr.stnx = (double*) malloc(nstns * sizeof(double)))) + { + printf("Error allocating stnx array... exiting\n"); + exit(1); + } + if (!(tairstr.stny = (double*) malloc(nstns * sizeof(double)))) + { + printf("Error allocating stny array... exiting\n"); + exit(1); + } + if (!(tairstr.stnz = (double*) malloc(nstns * sizeof(double)))) + { + printf("Error allocating stnz array... exiting\n"); + exit(1); + } + if (!(intstr.stnx = (double*) malloc(nstns * sizeof(double)))) + { + printf("Error allocating stnx array... exiting\n"); + exit(1); + } + if (!(intstr.stny = (double*) malloc(nstns * sizeof(double)))) + { + printf("Error allocating stny array... exiting\n"); + exit(1); + } + if (!(intstr.sqdist = (double*) malloc(nstns * sizeof(double)))) + { + printf("Error allocating sqdist array... exiting\n"); + exit(1); + } + if (!(intstr.wt = (double*) malloc(nstns * sizeof(double)))) + { + printf("Error allocating wt array... exiting\n"); + exit(1); + } + if (!(intstr.id = (short*) malloc(nstns * sizeof(short)))) + { + printf("Error allocating id array... exiting\n"); + exit(1); + } + /* index array needed for bubble sort algorithm */ + if (!(bsi = (short*) malloc(nstns * sizeof(short)))) {printf("Couild not malloc bsi\n"); exit(1);} + /* istr and tstr share space for wt and id arrays */ + tairstr.wt = intstr.wt; + tairstr.id = intstr.id; + if (inflag) + { + if (!(in = (int*) malloc(nstns * sizeof(int)))) + { + printf("Error allocating in array... exiting\n"); + exit(1); + } + } + if (!(snot = (int*) malloc(nstns * sizeof(int)))) + { + printf("Error allocating snot array... exiting\n"); + exit(1); + } + + /* allocate memory for output arrays */ + if (!(ngooddays = (double*) malloc(nstns * sizeof(double)))) + { + printf("Error allocating ngooddays array... exiting\n"); + exit(1); + } + if (!(ngoodstns = (double*) malloc(ndays * sizeof(double)))) + { + printf("Error allocating ngoodstns array... exiting\n"); + exit(1); + } + if (!(xlrdaysum = (double*) malloc(ndays * sizeof(double)))) + { + printf("Error allocating xlrdaysum array... exiting\n"); + exit(1); + } + if (!(xlrdaysumsq = (double*) malloc(ndays * sizeof(double)))) + { + printf("Error allocating xlrdaysumsq array... exiting\n"); + exit(1); + } + if (!(ylrdaysum = (double*) malloc(ndays * sizeof(double)))) + { + printf("Error allocating ylrdaysum array... exiting\n"); + exit(1); + } + if (!(ylrdaysumsq = (double*) malloc(ndays * sizeof(double)))) + { + printf("Error allocating ylrdaysumsq array... exiting\n"); + exit(1); + } + if (!(zlrdaysum = (double*) malloc(ndays * sizeof(double)))) + { + printf("Error allocating zlrdaysum array... exiting\n"); + exit(1); + } + if (!(zlrdaysumsq = (double*) malloc(ndays * sizeof(double)))) + { + printf("Error allocating zlrdaysumsq array... exiting\n"); + exit(1); + } + + /* read stnmeta data and fill appropriate arrays */ + nstns_nws = nstns_snot = 0; + nstns_intile = 0; + for (i=0 ; i= ncols) || + (stn_filey < 0) || (stn_filey >= nrows) || + (!mask_array[stn_offset])) + { + tairstr.stnz[i] = metastr.z = elevs[i]; + continue; + } + + /* station is in mask, and inside DEM data region */ + /* begin loops through kernel to calculate filtered elevation */ + value = 0.0; + sum_wt = 0.0; + for (k=0 ; k= nrows)) + continue; + + koff = k * width; + img_off = y * ncols; + + for (l=0 ; l= ncols)) + continue; + + koffset = koff + l; + img_offset = img_off + x; + + /* finally, check mask for convolution overlap */ + if (!mask_array[img_offset]) + continue; + + /* accumulate convolved values */ + value += dem_array[img_offset] * kernel[koffset]; + sum_wt += kernel[koffset]; + } + } + + /* check that at least some good values were found */ + if (!sum_wt) + { + tairstr.stnz[i] = metastr.z = elevs[i]; + continue; + } + tairstr.stnz[i] = value/sum_wt; + } + else /* fzwidth == 0 */ + { + tairstr.stnz[i] = metastr.z = elevs[i]; + } + + } /* end of nstns loop */ + + if (nstns_intile) + { + /* allocate space for xval ouput */ + if (!(out_stnnames = (char*)malloc(nstns_intile * descriptor_length * sizeof(char)))) + { + printf("Error allocating for out_stnnames array.\n"); + exit(1); + } + if (!(out_stnids = (char*)malloc(nstns_intile * descriptor_length * sizeof(char)))) + { + printf("Error allocating for out_stnids array.\n"); + exit(1); + } + if (!(out_stnx = (double*)malloc(nstns_intile * sizeof(double)))) + { + printf("Error allocating for out_stnx array... exiting\n"); + exit(1); + } + if (!(out_stny = (double*)malloc(nstns_intile * sizeof(double)))) + { + printf("Error allocating for out_stny array... exiting\n"); + exit(1); + } + if (!(out_stnz = (double*)malloc(nstns_intile * sizeof(double)))) + { + printf("Error allocating for out_stnz array... exiting\n"); + exit(1); + } + if (!(out_obs = (double*)malloc(nstns_intile*ndays* sizeof(double)))) + { + printf("Error allocating for out_obs array... exiting\n"); + exit(1); + } + if (!(out_pred = (double*)malloc(nstns_intile*ndays* sizeof(double)))) + { + printf("Error allocating for out_pred array... exiting\n"); + exit(1); + } + if (!(out_good = (int*)malloc(nstns_intile*ndays* sizeof(int)))) + { + printf("Error allocating for out_good array... exiting\n"); + exit(1); + } + + /* netcdf define dimensions and array variables for output arrays */ + len_stns = nstns_intile; + len_days = ndays; + nc_def_dim(ncid_out, "stns", len_stns, &outdid_stns); + nc_def_dim(ncid_out, "descriptor_length", descriptor_length, &outid_desc); + nc_def_dim(ncid_out, "days", len_days, &outdid_days); + + stndims[0] = outdid_stns; + stndims[1] = outid_desc; + out_did[0] = outdid_stns; + out_did[1] = outdid_days; + out_did1[0] = outdid_stns; + + nc_def_var (ncid_out, "station_name", NC_CHAR, 2, stndims, &outid_stnname); + nc_def_var (ncid_out, "station_id", NC_CHAR, 2, stndims, &outid_stnid); + + nc_def_var(ncid_out, "stnx", NC_DOUBLE, 1, out_did1, &outid_stnx); + strcpy(att,"station projected x coordinate"); + nc_put_att_text(ncid_out, outid_stnx, "longname", strlen(att), att); + strcpy(att,"m"); + nc_put_att_text(ncid_out, outid_stnx, "units", strlen(att), att); + + nc_def_var(ncid_out, "stny", NC_DOUBLE, 1, out_did1, &outid_stny); + strcpy(att,"station projected y coordinate"); + nc_put_att_text(ncid_out, outid_stny, "longname", strlen(att), att); + strcpy(att,"m"); + nc_put_att_text(ncid_out, outid_stny, "units", strlen(att), att); + + nc_def_var(ncid_out, "stnz", NC_DOUBLE, 1, out_did1, &outid_stnz); + strcpy(att,"station elevation"); + nc_put_att_text(ncid_out, outid_stnz, "longname", strlen(att), att); + strcpy(att,"m"); + nc_put_att_text(ncid_out, outid_stnz, "units", strlen(att), att); + + nc_def_var(ncid_out, "obs", NC_DOUBLE, 2, out_did, &outid_obs); + strcpy(att,"observed temperature"); + nc_put_att_text(ncid_out, outid_obs, "longname", strlen(att), att); + strcpy(att,"degrees C"); + nc_put_att_text(ncid_out, outid_obs, "units", strlen(att), att); + + nc_def_var(ncid_out, "pred", NC_DOUBLE, 2, out_did, &outid_pred); + strcpy(att,"predicted temperature"); + nc_put_att_text(ncid_out, outid_pred, "longname", strlen(att), att); + strcpy(att,"degrees C"); + nc_put_att_text(ncid_out, outid_pred, "units", strlen(att), att); + + nc_def_var(ncid_out, "good", NC_INT, 2, out_did, &outid_good); + strcpy(att,"good value flag (0=missing obs, 1=good obs)"); + nc_put_att_text(ncid_out, outid_good, "longname", strlen(att), att); + strcpy(att,"flag"); + nc_put_att_text(ncid_out, outid_good, "units", strlen(att), att); + } /* end if nstns_intile */ + + /* take output file out of define mode, ready for write */ + nc_enddef(ncid_out); + + /* read station daily data and fill appropriate arrays */ + for (i=0 ; i 49) i = 49; + frad = f90[i]; + + /* calculate initial density filter parameters */ + intstr.dfsqrad = dfirad * dfirad; + intstr.inv_dfsqrad = 1.0 / intstr.dfsqrad; + intstr.dfarea = PI * intstr.dfsqrad; + intstr.trunc = exp(-intstr.gsp); + intstr.dftrunc = exp(-DFGSP); + intstr.dfavgwt = ((1.0 - intstr.dftrunc)/DFGSP) - intstr.dftrunc; + + /* start loop through stations */ + mae1 = bias1 = 0.0; + mae2 = bias2 = 0.0; + tairstr.switch_init = 1; + n_in = 0; + ngood = 0.0; + nlr = 0.0; + ngoodlr = 0; + for (i=0 ; ij ; k--) { + if (intstr.sqdist[bsi[k]] < intstr.sqdist[bsi[k-1]]) { + temp_short = bsi[k]; + bsi[k]=bsi[k-1]; + bsi[k-1]=temp_short; + } + } + } + + /* set the weighting kernel radius to 1m greater than the furthest included station */ + /* the extra 1.0m ensures that none of the stations in the list will have wt=0 */ + intstr.sqrad=intstr.sqdist[bsi[ans]]+1.0; + + /* fill the list of included stations, being sure to leave out the prediction point + station */ + k=0; + wtsum=0.0; + for (j=0 ; j \n"); + exit(1); + } + /* get initialization data from command line */ + strcpy(stnmeta_nc_f.name,argv[1]); + strcpy(stndata_nc_f.name,argv[2]); + critp = atof(argv[3]); + dfirad = atof(argv[4]); + intstr.gsp = atof(argv[5]); + intstr.ans = atof(argv[6]); + smwidth = atoi(argv[7]); + strcpy(outprefix,argv[8]); + + +/* ----------------------- netcdf calls for station meta file --------------------------------------- */ +/* */ + if( status = nc_open(stnmeta_nc_f.name, NC_SHARE, &ncid_meta) ) + { + printf("Error opening %s for netcdf read, exiting\n", stnmeta_nc_f.name); + exit(1); + } + status = nc_inq_dimid (ncid_meta, "stations", &nstnid); + status = nc_inq_dimlen(ncid_meta, nstnid, &lenp); + meta_nc_hdr.nstns = lenp; + nstns = meta_nc_hdr.nstns; + status = nc_inq_dimid (ncid_meta, "year_day", &yrdid); + status = nc_inq_dimlen(ncid_meta, yrdid, &lenp); + year_day = lenp; + status = nc_inq_dimid (ncid_meta, "descriptor_length", &descid); + status = nc_inq_dimlen(ncid_meta, descid, &lenp); + descriptor_length = lenp; + + status = nc_get_att_text (ncid_meta, NC_GLOBAL, "year", ayear); + year = atoi(ayear); + status = nc_get_att_int (ncid_meta, NC_GLOBAL, "days_per_year", &year_day); + + status = nc_inq_varid (ncid_meta, "station_name", &stnnameid); + if (!(stnnames = (char*) malloc(nstns * descriptor_length * sizeof(char)))) exit(1); + status = nc_get_var_text (ncid_meta, stnnameid, stnnames); + + status = nc_inq_varid (ncid_meta, "station_id", &stnidid); + if (!(stnids = (char*) malloc(nstns * descriptor_length * sizeof(char)))) exit(1); + status = nc_get_var_text (ncid_meta, stnidid, stnids); + + status = nc_inq_varid (ncid_meta, "station_elevation", &stnelevid); + if (!(elevs = (double*) malloc(nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnelevid, elevs); + + status = nc_inq_varid (ncid_meta, "station_latitude", &stnlatid); + if (!(lats = (double*) malloc(nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnlatid, lats); + + status = nc_inq_varid (ncid_meta, "station_longitude", &stnlonid); + if (!(lons = (double*) malloc(nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnlonid, lons); + + status = nc_inq_varid (ncid_meta, "no_data", &no_dataid); + if (!(no_data = (double*) malloc(nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, no_dataid, no_data); + + status = nc_inq_varid (ncid_meta, "in_tile", &in_tileid); + if (!(in_tile = (int*) malloc(nstns * sizeof(int)))) exit(1); + status = nc_get_var_int (ncid_meta, in_tileid, in_tile); + + status = nc_inq_varid (ncid_meta, "stnx", &stnxid); + if (!(stnx = (double*) malloc(nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnxid, stnx); + + status = nc_inq_varid (ncid_meta, "stny", &stnyid); + if (!(stny = (double*) malloc(nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnyid, stny); + + status = nc_inq_varid (ncid_meta, "stnfx", &stnfxid); + if (!(stnfx = (double*) malloc(nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnfxid, stnfx); + + status = nc_inq_varid (ncid_meta, "stnfy", &stnfyid); + if (!(stnfy = (double*) malloc(nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnfyid, stnfy); + + tair_name=0; + status = nc_inq_varid (ncid_meta, "tmax", &tairid); + if (status) + { + tair_name=1; + status = nc_inq_varid (ncid_meta, "tmin", &tairid); + } + if (!(tair = (double*) malloc(nstns * year_day * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, tairid, tair); + + status = nc_get_att_int(ncid_meta, NC_GLOBAL, "xyloc_flag", &meta_nc_hdr.xyloc_flag); + + printf("xyloc_flag = %d\n",meta_nc_hdr.xyloc_flag); + + nc_close(ncid_meta); +/* */ +/* ------------------- end netcdf calls for station meta file --------------------------------------- */ + + /* set internal variables */ + nstns = tairstr.nstns = intstr.nstns = meta_nc_hdr.nstns; + ndays = tairstr.ndays = year_day; + tairstr.noz = 0; + + //BWM + //printf("nstns: %d ndays: %d\n", nstns, ndays); + /* allocate memory */ + if (!(good = (char*) malloc(nstns * ndays * sizeof(char)))) exit(1); + if (!(good2 = (char*) malloc(nstns * ndays * sizeof(char)))) exit(1); + if (!(mflag = (int*) malloc(nstns * ndays * sizeof(int)))) exit(1); + if (!(code = (int*) malloc(nstns * ndays * sizeof(int)))) exit(1); + if (!(keep = (int*) malloc(nstns * sizeof(int)))) exit(1); + if (!(tairstr.stnobs = (double*) malloc(nstns * ndays * sizeof(double)))) exit(1); + if (!(tairstr.stnsmobs = (double*) malloc(nstns * ndays * sizeof(double)))) exit(1); + if (!(tairstr.stnx = (double*) malloc(nstns * sizeof(double)))) exit(1); + if (!(tairstr.stny = (double*) malloc(nstns * sizeof(double)))) exit(1); + if (!(tairstr.stnz = (double*) malloc(nstns * sizeof(double)))) exit(1); + if (!(intstr.stnx = (double*) malloc(nstns * sizeof(double)))) exit(1); + if (!(intstr.stny = (double*) malloc(nstns * sizeof(double)))) exit(1); + if (!(intstr.sqdist = (double*) malloc(nstns * sizeof(double)))) exit(1); + if (!(intstr.wt = (double*) malloc(nstns * sizeof(double)))) exit(1); + if (!(intstr.id = (short*) malloc(nstns * sizeof(short)))) exit(1); + /* istr and tstr share space for wt and id arrays */ + tairstr.wt = intstr.wt; + tairstr.id = intstr.id; + + /* read stnmeta data and fill appropriate arrays */ + for (i=0 ; i= 10.0) good2[obs_off + day] = 0; + + } /* end of ndays loop */ + + /* free memory depending on count */ + free(listgood); + free(tairstr.listobs); + free(tairstr.listsmobs); + free(tairstr.listx); + free(tairstr.listy); + free(tairstr.listz); + for (j=0 ; j<3 ; j++) + { + free(tairstr.regx[j]); + } + free(tairstr.regx); + free(tairstr.regy); + free(tairstr.regwtall); + free(tairstr.regwt); + + /* switch initial sign for next station */ + tairstr.switch_init = -tairstr.switch_init; + + } /* end of nstations loop */ + + /* loop through stations and cast out all those which exceed the + missing days criteria after the xval filtering */ + nkeep = 0; + for (i=0 ; i= critp) + { + /* this station passes the filtering criteria, and will be + retained in the filtered stnmeta and stndata files */ + keep[i] = 1; + nkeep++; + } + } + printf("nkeep = %d\n",nkeep); + error_check=2; + + /* allocate output memory */ + if (!(stnnames_filt = (char*) malloc(nkeep * descriptor_length * sizeof(char)))) exit(1); + if (!(stnids_filt = (char*) malloc(nkeep * descriptor_length * sizeof(char)))) exit(1); + if (!(elevs_filt = (double*) malloc(nkeep * sizeof(double)))) exit(1); + if (!(lats_filt = (double*) malloc(nkeep * sizeof(double)))) exit(1); + if (!(lons_filt = (double*) malloc(nkeep * sizeof(double)))) exit(1); + if (!(no_data_filt = (double*) malloc(nkeep * sizeof(double)))) exit(1); + if (!(in_tile_filt = (int*) malloc(nkeep * sizeof(int)))) exit(1); + if (!(stnx_filt = (double*) malloc(nkeep * sizeof(double)))) exit(1); + if (!(stny_filt = (double*) malloc(nkeep * sizeof(double)))) exit(1); + if (!(stnfx_filt = (double*) malloc(nkeep * sizeof(double)))) exit(1); + if (!(stnfy_filt = (double*) malloc(nkeep * sizeof(double)))) exit(1); + if (!(tair_filt = (double*) malloc(nkeep * year_day * sizeof(double)))) exit(1); + if (!(mflag_filt = (int*) malloc(nkeep * year_day * sizeof(int)))) exit(1); + if (!(code_filt = (int*) malloc(nkeep * year_day * sizeof(int)))) exit(1); + + ii = 0; + + /* read through the meta and daily data, and write out the station records that + are being kept to the new meta and data output files */ + for (i=0 ; i \n"); + printf("argc = %d\n",argc); + exit(1); + } + /* get initialization data from command line */ + strcpy(tair_nc_f.name,argv[1]); + critp = atof(argv[2]); +/* ----------------------- netcdf calls for station meta file --------------------------------------- */ +/* */ + if( status = nc_open(tair_nc_f.name, NC_WRITE, &ncid_meta) ) + { + printf("Error opening %s for netcdf read, exiting\n", tair_nc_f.name); + exit(1); + } + status = nc_inq_dimid (ncid_meta, "stations", &nstnid); + status = nc_inq_dimlen(ncid_meta, nstnid, &lenp); + nstns = lenp; + status = nc_inq_dimid (ncid_meta, "year_day", &yrdid); + status = nc_inq_dimlen(ncid_meta, yrdid, &lenp); + year_day = lenp; + status = nc_inq_dimid (ncid_meta, "descriptor_length", &descid); + status = nc_inq_dimlen(ncid_meta, descid, &lenp); + descriptor_length = lenp; + status = nc_get_att_int (ncid_meta, NC_GLOBAL, "year", &year); + status = nc_get_att_int (ncid_meta, NC_GLOBAL, "days_per_year", &year_day); + ndays = year_day; + + status = nc_inq_varid (ncid_meta, "station_name", &stnnameid); + if (!(stnnames = (char*) malloc(nstns * descriptor_length * sizeof(char)))) exit(1); + status = nc_get_var_text (ncid_meta, stnnameid, stnnames); + + status = nc_inq_varid (ncid_meta, "station_id", &stnidid); + if (!(stnids = (char*) malloc(nstns * descriptor_length * sizeof(char)))) exit(1); + status = nc_get_var_text (ncid_meta, stnidid, stnids); + + status = nc_inq_varid (ncid_meta, "station_elevation", &stnelevid); + if (!(elevs = (double*) malloc(nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnelevid, elevs); + + status = nc_inq_varid (ncid_meta, "station_latitude", &stnlatid); + if (!(lats = (double*) malloc(nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnlatid, lats); + + status = nc_inq_varid (ncid_meta, "station_longitude", &stnlonid); + if (!(lons = (double*) malloc(nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnlonid, lons); + + status = nc_inq_varid (ncid_meta, "no_data", &no_dataid); + if (!(no_data = (double*) malloc(nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, no_dataid, no_data); + + status = nc_inq_varid (ncid_meta, "in_tile", &in_tileid); + if (!(in_tile = (int*) malloc(nstns * sizeof(int)))) exit(1); + status = nc_get_var_int (ncid_meta, in_tileid, in_tile); + + status = nc_inq_varid (ncid_meta, "stnx", &stnxid); + if (!(stnx = (double*) malloc(nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnxid, stnx); + + status = nc_inq_varid (ncid_meta, "stny", &stnyid); + if (!(stny = (double*) malloc(nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnyid, stny); + + status = nc_get_att_int (ncid_meta, NC_GLOBAL, "xyloc_flag", &xyloc_flag); + if (xyloc_flag == 1) + { + status = nc_inq_varid (ncid_meta, "stnfx", &stnfxid); + if (!(stnfx = (double*) malloc(nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnfxid, stnfx); + + status = nc_inq_varid (ncid_meta, "stnfy", &stnfyid); + if (!(stnfy = (double*) malloc(nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnfyid, stnfy); + } + + tair_name=0; + status = nc_inq_varid (ncid_meta, "tmax", &tairid); + if (status) + { + tair_name=1; + status = nc_inq_varid (ncid_meta, "tmin", &tairid); + } + if (!(tair = (double*) malloc(nstns * year_day * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, tairid, tair); + + nc_close(ncid_meta); +/* */ +/* ------------------- end netcdf calls for station meta file --------------------------------------- */ + /* allocate memory */ + if (!(good = (char*) malloc(nstns * ndays * sizeof(char)))) exit(1); + if (!(mflag = (int*) malloc(nstns * ndays * sizeof(int)))) exit(1); + if (!(keep = (int*) malloc(nstns * sizeof(int)))) exit(1); + /* read station daily data and fill appropriate arrays */ + for (i=0 ; i= critp) + { + /* this station passes the filtering criteria, and will be + retained in the filtered stnmeta and stndata files */ + keep[i] = 1; + nkeep++; + } + } + /* allocate output memory */ + if (!(stnnames_filt = (char*) malloc(nkeep * descriptor_length * sizeof(char)))) exit(1); + if (!(stnids_filt = (char*) malloc(nkeep * descriptor_length * sizeof(char)))) exit(1); + if (!(elevs_filt = (double*) malloc(nkeep * sizeof(double)))) exit(1); + if (!(lats_filt = (double*) malloc(nkeep * sizeof(double)))) exit(1); + if (!(lons_filt = (double*) malloc(nkeep * sizeof(double)))) exit(1); + if (!(no_data_filt = (double*) malloc(nkeep * sizeof(double)))) exit(1); + if (!(in_tile_filt = (int*) malloc(nkeep * sizeof(int)))) exit(1); + if (!(stnx_filt = (double*) malloc(nkeep * sizeof(double)))) exit(1); + if (!(stny_filt = (double*) malloc(nkeep * sizeof(double)))) exit(1); + if (xyloc_flag) + { + if (!(stnfx_filt = (double*) malloc(nkeep * sizeof(double)))) exit(1); + if (!(stnfy_filt = (double*) malloc(nkeep * sizeof(double)))) exit(1); + } + if (!(tair_filt = (double*) malloc(nkeep * year_day * sizeof(double)))) exit(1); + if (!(mflag_filt = (int*) malloc(nkeep * year_day * sizeof(int)))) exit(1); + + /* read through the meta and daily data, and write out the station records that + are being kept to the new meta and data output files */ + ii = 0; + for (i=0 ; i \n"); + exit(1); + } + strcpy(mask_nc_f.name,argv[1]); + strcpy(stnmeta_nc_f.name, argv[2]); + strcpy(stndata_nc_f.name, argv[3]); + fzwidth = atof(argv[4]); + dfirad = atof(argv[5]); + intstr.gsp = atof(argv[6]); + intstr.ans = atof(argv[7]); + prcpstr.f_max = atof(argv[8]); + pop_crit = atof(argv[9]); + smwidth = atoi(argv[10]); + +/* ----------------------- netcdf calls for mask file --------------------------------------- */ +/* */ + strcat(mask_nc_f.name,".nc"); + if( status = nc_open(mask_nc_f.name, 0, &ncid_mask) ) + { + printf("Error opening %s for netcdf read, exiting\n", mask_nc_f.name); + exit(1); + } + status = nc_inq_dimid (ncid_mask, "ncols", &colsid_mask); + status = nc_inq_dimlen(ncid_mask, colsid_mask, &lenp); + mask_nc_hdr.ncols = lenp; + status = nc_inq_dimid (ncid_mask, "nrows", &rowsid_mask); + status = nc_inq_dimlen(ncid_mask, rowsid_mask, &lenp); + mask_nc_hdr.nrows = lenp; + status = nc_get_att_text (ncid_mask, NC_GLOBAL, "tag", mask_nc_hdr.tag); + status = nc_get_att_double (ncid_mask, NC_GLOBAL, "ulx", &mask_nc_hdr.ulx); + status = nc_get_att_double (ncid_mask, NC_GLOBAL, "uly", &mask_nc_hdr.uly); + status = nc_get_att_double (ncid_mask, NC_GLOBAL, "cellsize", &mask_nc_hdr.cellsize); + status = nc_get_att_long (ncid_mask, NC_GLOBAL, "gctp_outsys", &mask_nc_hdr.gctp.outsys); + status = nc_get_att_long (ncid_mask, NC_GLOBAL, "gctp_outzone", &mask_nc_hdr.gctp.outzone); + status = nc_get_att_long (ncid_mask, NC_GLOBAL, "gctp_outdatum", &mask_nc_hdr.gctp.outdatum); + status = nc_get_att_double (ncid_mask, NC_GLOBAL, "gctp_outparm", mask_nc_hdr.gctp.outparm); + status = nc_get_att_text (ncid_mask, NC_GLOBAL, "gctp_fn27", mask_nc_hdr.gctp.fn27); + status = nc_get_att_text (ncid_mask, NC_GLOBAL, "gctp_fn83", mask_nc_hdr.gctp.fn83); + if (strncmp(mask_nc_hdr.tag,IMGTAG,7)) + { + printf("mask file not a registered image. Use register_image first.\n"); + exit(1); + } + + status = nc_inq_varid (ncid_mask, "image", &imageid_mask); + if (!(mask_array = (long int*) malloc(mask_nc_hdr.nrows * mask_nc_hdr.ncols * sizeof(long int)))) exit(1); + status = nc_get_var_long (ncid_mask, imageid_mask, mask_array); + + nc_close(ncid_mask); +/* */ +/* ----------------------- end netcdf calls for mask file --------------------------------------- */ + +/* ----------------------- netcdf calls for station meta file --------------------------------------- */ +/* */ + if( status = nc_open(stnmeta_nc_f.name, NC_WRITE, &ncid_meta) ) + { + printf("Error opening %s for netcdf read, exiting\n", stnmeta_nc_f.name); + exit(1); + } + status = nc_inq_dimid (ncid_meta, "stations", &nstnid); + status = nc_inq_dimlen(ncid_meta, nstnid, &lenp); + meta_nc_hdr.nstns = lenp; + status = nc_inq_dimid (ncid_meta, "year_day", &yrdid); + status = nc_inq_dimlen(ncid_meta, yrdid, &lenp); + year_day = lenp; + status = nc_inq_dimid (ncid_meta, "descriptor_length", &descid); + status = nc_inq_dimlen(ncid_meta, descid, &lenp); + descriptor_length = lenp; + + status = nc_inq_varid (ncid_meta, "station_name", &stnnameid); + status = nc_inq_varid (ncid_meta, "station_id", &stnidid); + + status = nc_inq_varid (ncid_meta, "station_elevation", &stnelevid); + if (!(elevs = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnelevid, elevs); + + status = nc_inq_varid (ncid_meta, "station_latitude", &stnlatid); + if (!(lats = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnlatid, lats); + + status = nc_inq_varid (ncid_meta, "station_longitude", &stnlonid); + if (!(lons = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnlonid, lons); + + status = nc_inq_varid (ncid_meta, "stnx", &stnxid); + if (!(stnx = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnxid, stnx); + + status = nc_inq_varid (ncid_meta, "stny", &stnyid); + if (!(stny = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnyid, stny); + + status = nc_get_att_int (ncid_meta, NC_GLOBAL, "xyloc_flag", &xyloc_flag); + if (xyloc_flag == 1) + { + status = nc_inq_varid (ncid_meta, "stnfx", &stnfxid); + if (!(stnfx = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnfxid, stnfx); + + status = nc_inq_varid (ncid_meta, "stnfy", &stnfyid); + if (!(stnfy = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnfyid, stnfy); + } + + status = nc_inq_varid (ncid_meta, "mflag", &mflagid); + if (!(mflag = (int*) malloc(meta_nc_hdr.nstns * year_day * sizeof(int)))) + { + printf("Error allocating mflag array.\n"); + exit(1); + } + status = nc_get_var_int (ncid_meta, mflagid, mflag); + + status = nc_inq_varid (ncid_meta, "precip", &precipid); + if (!(precip = (double*) malloc(meta_nc_hdr.nstns * year_day * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, precipid, precip); + + meta_nc_hdr.zfilt_flag = 0; +/* */ +/* ------------------- end netcdf calls for station meta file --------------------------------------- */ + + /* if z-filtering width other than 0.0, open the special ini file + containing filter type and mask and DEM filenames */ + if (fzwidth) + { + strcpy(ini.name, "zfilter.ini"); + if (file_open(&ini, 'i')) + { + printf("Error opening %s as special zfilter ini file\n",ini.name); + exit(1); + } + if (scan_value(ini, &type, 'i')) exit(1); + if (scan_open(ini, &mask_f, 'r')) exit(1); + if (scan_open(ini, &dem_f, 'r')) exit(1); + fclose(ini.ptr); + } + + /* perform error checking on input file headers */ + if (fzwidth) + { + /* read mask file header and verify */ + fread(&maskhdr,sizeof(imghdr_struct), 1, mask_f.ptr); + if (strcmp(maskhdr.tag,IMGTAG)) + { + printf("%s not a registered image. Use register_image first.\n",mask_f.name); + exit(1); + } + + /* read dem file header and verify */ + fread(&demhdr,sizeof(imghdr_struct), 1, dem_f.ptr); + if (strcmp(demhdr.tag,IMGTAG)) + { + printf("%s not a registered image. Use register_image first.\n",dem_f.name); + exit(1); + } + + /* compare headers for DEM and mask */ + test = (maskhdr.ncols == demhdr.ncols)*(maskhdr.nrows == demhdr.nrows)* + (maskhdr.ulx == demhdr.ulx)*(maskhdr.uly == demhdr.uly)* + (maskhdr.cellsize == demhdr.cellsize); + if ((!test) || memcmp(&maskhdr.gctp,&demhdr.gctp,sizeof(gctp_struct))) + { + printf("Header information differs between mask and DEM...exiting\n"); + exit(1); + } + } + + if (fzwidth) + { + /* compare gctp structures between images and stnmeta */ + if (memcmp(&maskhdr.gctp,&metahdr.gctp,sizeof(gctp_struct))) + { + printf("Different projection parameters in images and stnmeta.\n"); + exit(1); + } + } + + /* set internal variables */ + if (fzwidth) + { + ncols = maskhdr.ncols; + nrows = maskhdr.nrows; + cellsize = demhdr.cellsize; + ul_mapx = demhdr.ulx; + ul_mapy = demhdr.uly; + } + nstns = prcpstr.nstns = intstr.nstns = meta_nc_hdr.nstns; + ndays = prcpstr.ndays = year_day; + prcpstr.noz = 0; + + if (fzwidth) + { + /* allocate zfilter memory */ + if (!(dem_array = (short*) malloc(nrows*ncols*sizeof(short)))) exit(1); + if (!(mask_array = (long int*) malloc(nrows*ncols*sizeof(long int)))) exit(1); + + /* read and close dem and mask files */ + fread(dem_array,sizeof(short),nrows*ncols,dem_f.ptr); + fread(mask_array,sizeof(long int),nrows*ncols,mask_f.ptr); + fclose(dem_f.ptr); + fclose(mask_f.ptr); + + /* set size of elevation filter window */ + width = (int) (fzwidth / cellsize); + if (!(width % 2)) width++; + hw = width/2; + + /* generate filter kernel */ + if (!(kernel = (double*) malloc(width * width * sizeof(double)))) + { + printf("Error allocating for kernel... exiting\n"); + exit(1); + } + if (make_kernel(kernel, width, type)) + { + printf("Error in call to make_kernel() ... exiting\n"); + exit(1); + } + } + + /* allocate memory */ + if (!(good = (char*) malloc(nstns * ndays * sizeof(char)))) exit(1); + if (!(prcpstr.stnobs = (double*) malloc(nstns * ndays * sizeof(double)))) exit(1); + if (!(prcpstr.stnsmobs = (double*) malloc(nstns * ndays * sizeof(double)))) exit(1); + if (!(prcpstr.stnx = (double*) malloc(nstns * sizeof(double)))) exit(1); + if (!(prcpstr.stny = (double*) malloc(nstns * sizeof(double)))) exit(1); + if (!(prcpstr.stnz = (double*) malloc(nstns * sizeof(double)))) exit(1); + if (!(intstr.stnx = (double*) malloc(nstns * sizeof(double)))) exit(1); + if (!(intstr.stny = (double*) malloc(nstns * sizeof(double)))) exit(1); + if (!(intstr.sqdist = (double*) malloc(nstns * sizeof(double)))) exit(1); + if (!(intstr.wt = (double*) malloc(nstns * sizeof(double)))) exit(1); + if (!(intstr.id = (short*) malloc(nstns * sizeof(short)))) exit(1); + + if (!(bsi = (short*) malloc(nstns * sizeof (short)))) exit(1); + + /* istr and tstr share space for wt and id arrays */ + prcpstr.wt = intstr.wt; + prcpstr.id = intstr.id; + + /* read stnmeta data and fill appropriate arrays */ + for (i=0 ; i= ncols) || + (stn_filey < 0) || (stn_filey >= nrows) || + (!mask_array[stn_offset])) + { + prcpstr.stnz[i] = elevs[i]; + continue; + } + + /* station is in mask, and inside DEM data region */ + /* begin loops through kernel to calculate filtered elevation */ + value = 0.0; + sum_wt = 0.0; + for (k=0 ; k= nrows)) + continue; + + koff = k * width; + img_off = y * ncols; + + for (l=0 ; l= ncols)) + continue; + + koffset = koff + l; + img_offset = img_off + x; + + /* finally, check mask for convolution overlap */ + if (!mask_array[img_offset]) + continue; + + /* accumulate convolved values */ + value += dem_array[img_offset] * kernel[koffset]; + sum_wt += kernel[koffset]; + } + } + + /* check that at least some good values were found */ + if (!sum_wt) + { + prcpstr.stnz[i] = elevs[i]; + continue; + } + prcpstr.stnz[i] = value/sum_wt; + } + else /* fzwidth == 0 */ + { + prcpstr.stnz[i] = elevs[i]; + } + + } /* end of nstns loop */ + + /* read station daily data and fill appropriate arrays */ + for (i=0 ; ij ; k--) + { + if (good[bsi[k]*ndays+day] && ((intstr.sqdist[bsi[k]] < intstr.sqdist[bsi[k-1]]) || !good[bsi[k-1]*ndays+day])) + { + /* switch these two stations in bsi */ + temp_short = bsi[k]; + bsi[k]=bsi[k-1]; + bsi[k-1]=temp_short; + } + } + } + + /* PET (4/12/2020): Count the number of good stations for the day */ + /* up to a max of ANS */ + ngood_today = 0; + for (j=0 ; j 0.0)); + nwetstns += wet; + } + } + + if (sumpop) pop /= sumpop; + /* otherwise, no good data for the day, pop = 0.0, and assume + prcp = 0.0 */ + + /* if precipitation occurrence threshold exceeded, predict + precipitation amount */ + if (pop > pop_crit) + { + // /* only generate regression stats when nwetstns > 3 */ + // if (nwetstns > 3) + // { + // /* generate regression arrays that change by day */ + // if (xv_prcp_regr_y_2switch(listgood, &prcpstr)) + // { + // printf("Error in xv_prcp_regr_y_switch() ... exiting\n"); + // exit(1); + // } + + // /* weighted least squares slope and intercept */ + // if (prcp_mlr(&prcpstr)) + // { + // printf("Error in wt_regr() ... exiting\n"); + // exit(1); + // } + + // } /* end if nwet > 3 */ + + // else + // { + prcpstr.coef[0] = prcpstr.coef[1] = prcpstr.coef[2] = + prcpstr.coef[3] = 0.0; + // } + + /* predict prcp for the day */ + if (xv_predict_prcp_noint(listgood, &prcpstr)) + { + printf("Error in xv_predict_prcp_noint() ... exiting\n"); + exit(1); + } + + } /* end if pop > pop_crit */ + + precip[obs_off+day] = prcpstr.pp; + + /* free memory depending on count */ + free(listgood); + free(prcpstr.listobs); + free(prcpstr.listsmobs); + free(prcpstr.listx); + free(prcpstr.listy); + free(prcpstr.listz); + for (j=0 ; j<3 ; j++) + { + free(prcpstr.regx[j]); + } + free(prcpstr.regx); + free(prcpstr.regy); + free(prcpstr.regwtall); + free(prcpstr.regwt); + + } /* end of nsimdays loop */ + + + /* switch initial sign for next station */ + prcpstr.switch_init = -prcpstr.switch_init; + + } /* end of nstations loop */ + +/* ----------------------- netcdf calls for station meta file update --------------------------------------- */ +/* */ + status = nc_redef(ncid_meta); + + gooddims[0] = nstnid; + gooddims[1] = yrdid; + status = nc_def_var(ncid_meta, "good", NC_BYTE, 2, gooddims, &goodid); + + status = nc_put_att_int(ncid_meta, NC_GLOBAL, "fill_flag", NC_INT, 1, &fill_flag); + + status = nc_enddef(ncid_meta); + + status = nc_put_var_double(ncid_meta, precipid, precip); + status = nc_put_var_uchar(ncid_meta, goodid, good); + +/* */ +/* ----------------------- end netcdf calls for station meta file update ----------------------------------- */ + + /* close files and exit */ + nc_close(ncid_meta); + + /* free remaining memory */ + if (fzwidth) + { + free(dem_array); + free(mask_array); + free(kernel); + } + free(good); + free(prcpstr.stnobs); + free(prcpstr.stnsmobs); + free(prcpstr.stnz); + free(intstr.stnx); + free(intstr.stny); + free(intstr.sqdist); + free(intstr.wt); + free(intstr.id); + + return(0); + +} /* end of main */ diff --git a/xval.d/fill_daily_tair3_nc.c b/xval.d/fill_daily_tair3_nc.c new file mode 100644 index 0000000..61286b9 --- /dev/null +++ b/xval.d/fill_daily_tair3_nc.c @@ -0,0 +1,759 @@ +/* +fill_daily_tair3_nc.c +Peter Thornton +11/11/96 + +Fills in the missing data values in a station data file using a specified +set of interpolation and prediction parameters. +The missing data flags remain in the data file, and so these missing +values can be filled in again (overwritten) using other interpolation and/or +prediction parameters at a later time if desired. + +Updated: +9/5/00, PET: from fill_daily_tair.c, added 3-variable regression +July 2005 : netcdf version - CSH +4/12/2020, PET: reduce number of stations when not enough good data (short station lists like PR) +*/ + +#include "metsrc2.h" +#include "dailywx.h" +#include "cproj.h" +#include "proj.h" +#include "netcdf.h" + +int main(int argc, char *argv[]) +{ + /* variable declarations */ + interpolate_struct intstr; + stnmetahdr_struct metahdr; + stnmeta_struct metastr; + daydathdr_struct datahdr; + dat_struct datastr; + daymet_tair_struct tairstr; + imghdr_struct maskhdr, demhdr; + gctp_struct gctp; + long iflg; + + char *good; + int gooddims[2]; + char *listgood; + + file ini; + file stnmeta_f, stndata_f; + file mask_f, dem_f; + file tmp_f; + + char endkey[8]; + char round[32]; + char outprefix[80]; + char jobstr[80]; + + short *dem_array; + short count; + + int ndays,nstns,nregr; + int i,j,k,l; + int day; + int x,y; + int type; + int nrows, ncols; + int width, hw; + int stn_filex, stn_filey; + int test; + int smwidth; + int ngoodpts; + int sstn; + int n_fill; + int job; + + long off; + long int obs_off, obs_offset; + long img_off, img_offset, stn_offset, koff, koffset; + + double dfirad; + double ul_mapx, ul_mapy; + double cellsize; + double fzwidth; + double *kernel; + double stn_mapx, stn_mapy; + double f_x, f_y; + double value,sum_wt; + double t1; + + size_t lenp; + int xyloc_flag; + stnmetahdr_struct meta_nc_hdr; + imghdr_struct mask_nc_hdr, dem_nc_hdr; + long int *mask_array; + file mask_nc_f, stnmeta_nc_f; + file stndata_nc_f; + int status, ncid; + int nstnid, descid, yrdid; + int stnnameid, stnidid, stnelevid, stnlatid, stnlonid; + int stnxid, stnyid, stnzid, goodid; + int stnfxid, stnfyid, stnfzid; + int stntypeid, stnuc1id, stnuc2id; + int tairid; + int ncid_mask, colsid_mask, rowsid_mask, imageid_mask; + int ncid_meta; + int *mflag, *code; + int mflagid, codeid; + int year_day, descriptor_length; + double *elevs, *lats, *lons, *tair; + double *stnx, *stny, *stnfx, *stnfy; + int fill_flag = 1; + + short *bsi; + short temp_short; + double wtsum; + int ngood_today; + + /* check for appropriate number of command line arguments */ + if (argc != 9) + { + printf("Usage: \n"); + printf(" \n"); + exit(1); + } + /* assuming command-line parameter passing */ + /* scan command line arguments into parameters */ + /* scan next parameter as the mask filename */ + strcpy(mask_nc_f.name,argv[1]); + strcpy(stnmeta_nc_f.name, argv[2]); + strcpy(stndata_nc_f.name, argv[3]); + fzwidth = atof(argv[4]); + dfirad = atof(argv[5]); + intstr.gsp = atof(argv[6]); + intstr.ans = atof(argv[7]); + smwidth = atoi(argv[8]); + +/* ----------------------- netcdf calls for mask file --------------------------------------- */ +/* */ + strcat(mask_nc_f.name,".nc"); + if( status = nc_open(mask_nc_f.name, 0, &ncid_mask) ) + { + printf("Error opening %s for netcdf read, exiting\n", mask_nc_f.name); + exit(1); + } + status = nc_inq_dimid (ncid_mask, "ncols", &colsid_mask); + status = nc_inq_dimlen(ncid_mask, colsid_mask, &lenp); + mask_nc_hdr.ncols = lenp; + status = nc_inq_dimid (ncid_mask, "nrows", &rowsid_mask); + status = nc_inq_dimlen(ncid_mask, rowsid_mask, &lenp); + mask_nc_hdr.nrows = lenp; + status = nc_get_att_text (ncid_mask, NC_GLOBAL, "tag", mask_nc_hdr.tag); + status = nc_get_att_double (ncid_mask, NC_GLOBAL, "ulx", &mask_nc_hdr.ulx); + status = nc_get_att_double (ncid_mask, NC_GLOBAL, "uly", &mask_nc_hdr.uly); + status = nc_get_att_double (ncid_mask, NC_GLOBAL, "cellsize", &mask_nc_hdr.cellsize); + status = nc_get_att_long (ncid_mask, NC_GLOBAL, "gctp_outsys", &mask_nc_hdr.gctp.outsys); + status = nc_get_att_long (ncid_mask, NC_GLOBAL, "gctp_outzone", &mask_nc_hdr.gctp.outzone); + status = nc_get_att_long (ncid_mask, NC_GLOBAL, "gctp_outdatum", &mask_nc_hdr.gctp.outdatum); + status = nc_get_att_double (ncid_mask, NC_GLOBAL, "gctp_outparm", mask_nc_hdr.gctp.outparm); + status = nc_get_att_text (ncid_mask, NC_GLOBAL, "gctp_fn27", mask_nc_hdr.gctp.fn27); + status = nc_get_att_text (ncid_mask, NC_GLOBAL, "gctp_fn83", mask_nc_hdr.gctp.fn83); + if (strncmp(mask_nc_hdr.tag,IMGTAG,7)) + { + printf("mask file not a registered image. Use register_image first.\n"); + exit(1); + } + + status = nc_inq_varid (ncid_mask, "image", &imageid_mask); + if (!(mask_array = (long int*) malloc(mask_nc_hdr.nrows * mask_nc_hdr.ncols * sizeof(long int)))) exit(1); + status = nc_get_var_long (ncid_mask, imageid_mask, mask_array); + + nc_close(ncid_mask); +/* */ +/* ----------------------- end netcdf calls for mask file --------------------------------------- */ + +/* ----------------------- netcdf calls for station meta file --------------------------------------- */ +/* */ + if( status = nc_open(stnmeta_nc_f.name, NC_WRITE, &ncid_meta) ) + { + printf("Error opening %s for netcdf read, exiting\n", stnmeta_nc_f.name); + exit(1); + } + status = nc_inq_dimid (ncid_meta, "stations", &nstnid); + status = nc_inq_dimlen(ncid_meta, nstnid, &lenp); + meta_nc_hdr.nstns = lenp; + status = nc_inq_dimid (ncid_meta, "year_day", &yrdid); + status = nc_inq_dimlen(ncid_meta, yrdid, &lenp); + year_day = lenp; + status = nc_inq_dimid (ncid_meta, "descriptor_length", &descid); + status = nc_inq_dimlen(ncid_meta, descid, &lenp); + descriptor_length = lenp; + + status = nc_inq_varid (ncid_meta, "station_name", &stnnameid); + status = nc_inq_varid (ncid_meta, "station_id", &stnidid); + + status = nc_inq_varid (ncid_meta, "station_elevation", &stnelevid); + if (!(elevs = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnelevid, elevs); + + status = nc_inq_varid (ncid_meta, "station_latitude", &stnlatid); + if (!(lats = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnlatid, lats); + + status = nc_inq_varid (ncid_meta, "station_longitude", &stnlonid); + if (!(lons = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnlonid, lons); + + status = nc_inq_varid (ncid_meta, "stnx", &stnxid); + if (!(stnx = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnxid, stnx); + + status = nc_inq_varid (ncid_meta, "stny", &stnyid); + if (!(stny = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnyid, stny); + + status = nc_get_att_int (ncid_meta, NC_GLOBAL, "xyloc_flag", &xyloc_flag); + if (xyloc_flag == 1) + { + status = nc_inq_varid (ncid_meta, "stnfx", &stnfxid); + if (!(stnfx = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnfxid, stnfx); + + status = nc_inq_varid (ncid_meta, "stnfy", &stnfyid); + if (!(stnfy = (double*) malloc(meta_nc_hdr.nstns * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, stnfyid, stnfy); + } + + status = nc_inq_varid (ncid_meta, "mflag", &mflagid); + if (!(mflag = (int*) malloc(meta_nc_hdr.nstns * year_day * sizeof(int)))) + { + printf("Error allocating mflag array.\n"); + exit(1); + } + status = nc_get_var_int (ncid_meta, mflagid, mflag); + + status = nc_inq_varid (ncid_meta, "tmax", &tairid); + if (status) + { + status = nc_inq_varid (ncid_meta, "tmin", &tairid); + } + if (!(tair = (double*) malloc(meta_nc_hdr.nstns * year_day * sizeof(double)))) exit(1); + status = nc_get_var_double (ncid_meta, tairid, tair); + + meta_nc_hdr.zfilt_flag = 0; +/* */ +/* ------------------- end netcdf calls for station meta file --------------------------------------- */ + + /* if z-filtering width other than 0.0, open the special ini file + containing filter type and mask and DEM filenames */ + if (fzwidth) + { + strcpy(ini.name, "zfilter.ini"); + if (file_open(&ini, 'i')) + { + printf("Error opening %s as special zfilter ini file\n",ini.name); + exit(1); + } + if (scan_value(ini, &type, 'i')) exit(1); + if (scan_open(ini, &mask_f, 'r')) exit(1); + if (scan_open(ini, &dem_f, 'r')) exit(1); + fclose(ini.ptr); + } + + /* perform error checking on input file headers */ + if (fzwidth) + { + /* read mask file header and verify */ + fread(&maskhdr,sizeof(imghdr_struct), 1, mask_f.ptr); + if (strcmp(maskhdr.tag,IMGTAG)) + { + printf("%s not a registered image. Use register_image first.\n",mask_f.name); + exit(1); + } + + /* read dem file header and verify */ + fread(&demhdr,sizeof(imghdr_struct), 1, dem_f.ptr); + if (strcmp(demhdr.tag,IMGTAG)) + { + printf("%s not a registered image. Use register_image first.\n",dem_f.name); + exit(1); + } + /* check that dem is sizeof(short) */ + + /* compare headers for DEM and mask */ + test = (maskhdr.ncols == demhdr.ncols)*(maskhdr.nrows == demhdr.nrows)* + (maskhdr.ulx == demhdr.ulx)*(maskhdr.uly == demhdr.uly)* + (maskhdr.cellsize == demhdr.cellsize); + if ((!test) || memcmp(&maskhdr.gctp,&demhdr.gctp,sizeof(gctp_struct))) + { + printf("Header information differs between mask and DEM...exiting\n"); + exit(1); + } + } + + if (fzwidth) + { + /* compare gctp structures between images and stnmeta */ + if (memcmp(&maskhdr.gctp,&metahdr.gctp,sizeof(gctp_struct))) + { + printf("Different projection parameters in images and stnmeta.\n"); + exit(1); + } + } + + /* set internal variables */ + if (fzwidth) + { + ncols = maskhdr.ncols; + nrows = maskhdr.nrows; + cellsize = demhdr.cellsize; + ul_mapx = demhdr.ulx; + ul_mapy = demhdr.uly; + } + nstns = tairstr.nstns = intstr.nstns = meta_nc_hdr.nstns; + printf("Number stations: %d\n", nstns); + ndays = tairstr.ndays = year_day; + tairstr.noz = 0; + + if (fzwidth) + { + /* allocate zfilter memory */ + if (!(dem_array = (short*) malloc(nrows*ncols*sizeof(short)))) exit(1); + if (!(mask_array = (long int*) malloc(nrows*ncols*sizeof(long int)))) exit(1); + + /* read and close dem and mask files */ + fread(dem_array,sizeof(short),nrows*ncols,dem_f.ptr); + fread(mask_array,sizeof(char),nrows*ncols,mask_f.ptr); + fclose(dem_f.ptr); + fclose(mask_f.ptr); + + /* set size of elevation filter window */ + width = (int) (fzwidth / cellsize); + if (!(width % 2)) width++; + hw = width/2; + + /* generate filter kernel */ + if (!(kernel = (double*) malloc(width * width * sizeof(double)))) + { + printf("Error allocating for kernel... exiting\n"); + exit(1); + } + if (make_kernel(kernel, width, type)) + { + printf("Error in call to make_kernel() ... exiting\n"); + exit(1); + } + } + + /* allocate memory */ + if (!(good = (char*) malloc(nstns * ndays * sizeof(char)))) exit(1); + if (!(tairstr.stnobs = (double*) malloc(nstns * ndays * sizeof(double)))) exit(1); + if (!(tairstr.stnsmobs = (double*) malloc(nstns * ndays * sizeof(double)))) exit(1); + if (!(tairstr.stnx = (double*) malloc(nstns * sizeof(double)))) exit(1); + if (!(tairstr.stny = (double*) malloc(nstns * sizeof(double)))) exit(1); + if (!(tairstr.stnz = (double*) malloc(nstns * sizeof(double)))) exit(1); + if (!(intstr.stnx = (double*) malloc(nstns * sizeof(double)))) exit(1); + if (!(intstr.stny = (double*) malloc(nstns * sizeof(double)))) exit(1); + if (!(intstr.sqdist = (double*) malloc(nstns * sizeof(double)))) exit(1); + if (!(intstr.wt = (double*) malloc(nstns * sizeof(double)))) exit(1); + if (!(intstr.id = (short*) malloc(nstns * sizeof(short)))) exit(1); + + if (!(bsi = (short*) malloc(nstns * sizeof(short)))) exit(1); + + /* istr and tstr share space for wt and id arrays */ + tairstr.wt = intstr.wt; + tairstr.id = intstr.id; + + /* read stnmeta data and fill appropriate arrays */ + for (i=0 ; i= ncols) || + (stn_filey < 0) || (stn_filey >= nrows) || + (!mask_array[stn_offset])) + { + tairstr.stnz[i] = elevs[i]; + continue; + } + + /* station is in mask, and inside DEM data region */ + /* begin loops through kernel to calculate filtered elevation */ + value = 0.0; + sum_wt = 0.0; + for (k=0 ; k= nrows)) + continue; + + koff = k * width; + img_off = y * ncols; + + for (l=0 ; l= ncols)) + continue; + + koffset = koff + l; + img_offset = img_off + x; + + /* finally, check mask for convolution overlap */ + if (!mask_array[img_offset]) + continue; + + /* accumulate convolved values */ + value += dem_array[img_offset] * kernel[koffset]; + sum_wt += kernel[koffset]; + } + } + + /* check that at least some good values were found */ + if (!sum_wt) + { + tairstr.stnz[i] = elevs[i]; + continue; + } + tairstr.stnz[i] = value/sum_wt; + } + else /* fzwidth == 0 */ + { + tairstr.stnz[i] = elevs[i]; + } + + } /* end of nstns loop */ + + /* read station daily data and fill appropriate arrays */ + for (i=0 ; ij ; k--) + { + if (good[bsi[k]*ndays+day] && ((intstr.sqdist[bsi[k]] < intstr.sqdist[bsi[k-1]]) || !good[bsi[k-1]*ndays+day])) + { + /* switch these two stations in bsi */ + temp_short = bsi[k]; + bsi[k]=bsi[k-1]; + bsi[k-1]=temp_short; + } + } + } + + /* PET (4/12/2020): Count the number of good stations for the day */ + /* up to a max of ANS */ + ngood_today = 0; + for (j=0 ; j= 3) + { + // /* generate regression arrays that change by day */ + // if (xv_tair_regr_y_2switch(listgood, &tairstr)) + // { + // printf("Error in xv_tair_regr_y_switch() ... exiting\n"); + // exit(1); + // } + + // /* weighted least squares slope and intercept */ + // if (xv_tair_mlr(&tairstr)) + // { + // printf("Error in wt_regr() ... exiting\n"); + // exit(1); + // } + + tairstr.coef[0] = tairstr.coef[1] = tairstr.coef[3] = 0.0; + tairstr.coef[2] = DEFAULT_TVSZ; + + /* predict tair for the day */ + if (xv_predict_tair_noint(listgood, &tairstr)) + { + printf("Error in xv_predict_tair_noint() ... exiting\n"); + exit(1); + } + } /* end if ngoodpts >= 3 */ + else if (ngoodpts) /* 1 or 2 good stations */ + { + /* find the closest of these stations */ + t1=0.0; + for (j=0 ; j t1) + { + sstn = j; + t1 = tairstr.wt[j]; + } + } + /* if a good station is found (should always be true) */ + if (t1) + { + tairstr.ta = tairstr.listobs[tairstr.dam_off + sstn] + + (DEFAULT_TVSZ * (tairstr.elev - tairstr.listz[sstn])); + } + else + { + printf("Error: ngoodpts > 0, but no good stations found\n"); + printf("This should never happen!\n"); + printf("exiting...\n"); + exit(1); + } + } + else /* no good satations */ + { + printf("Error: No good stations found\n"); + printf("Suggest using a larger value for average number of stations\n"); + printf("exiting...\n"); + for (j=0 ; jnstns; + sqdist = str->sqdist; + sqrad = str->sqrad; + gsp = str->gsp; + trunc = str->trunc; + inv_sqrad = 1.0/sqrad; + + count = 0; + wtsum = 0.0; + + for (i=0 ; iid[count] = i; + str->wt[count] = exp(-gsp * sqd * inv_sqrad) - trunc; + wtsum += str->wt[count]; + count++; + } + } + + /* normalize weights */ + if (count > 1) + { + for (i=0 ; iwt[i] /= wtsum; + } + str->count = count; + } + else + { + printf("fewer than 2 stations in list\n"); + ok = 0; + } + + return (!ok); +} + +int xv_fill_list(char* stmgood, char* listgood, daymet_prcp_struct* str); +int xv_fill_list(char* stmgood, char* listgood, daymet_prcp_struct* str) +/* fills the station-major list arrays (station-major -> day-major) +includes an array of good value flags for the cross-validation routine */ +{ + int ok=1; + int i,j; + int ndays,nstns; + long int stm_off; + + ndays = str->ndays; + nstns = str->count; + + for (i=0 ; iid[i] * ndays; + for (j=0 ; jlistobs[j*nstns + i] = str->stnobs[stm_off + j]; + str->listsmobs[j*nstns + i] = str->stnsmobs[stm_off + j]; + listgood[j*nstns + i] = stmgood[stm_off + j]; + } + str->listx[i] = str->stnx[str->id[i]]; + str->listy[i] = str->stny[str->id[i]]; + str->listz[i] = str->stnz[str->id[i]]; + } + + return (!ok); +} + + +int xv_prcp_regr_y_2switch(char* listgood, daymet_prcp_struct* str); +int xv_prcp_regr_y_2switch(char* listgood, daymet_prcp_struct* str) +/* make time variant regression array */ +/* generates the y variable for daily regressions */ +/* this version employs a switching mechanism to keep differencing unbiased */ +/* this version uses a double switching mechanism, alternates initial sign */ +/* PET (10/18/2015: adding always_good option since this is now being called with only good data */ +{ + int ok=1; + int nstns; + int sign_init, sign; + int stni, stnj, n=0; + double* prcp; + double prcp1, prcp2; + char* good; + int always_good = 1; + + nstns = str->count; + sign_init = str->switch_init; + prcp = str->listsmobs+str->dam_off; + good = listgood+str->dam_off; + + for (stni=0 ; stniregy[n] = sign * (prcp1-prcp2)/(prcp1+prcp2); */ + str->regy[n] = sign * (prcp1-prcp2); + str->regwt[n] = str->regwtall[n]; + } + else /* !prcp2 or missing data */ + { + str->regy[n] = str->regwt[n] = 0.0; + } + n++; + sign = -sign; + + } /* end for stnj */ + + } /* end if (prcp1) */ + + else /* !prcp1 or missing data */ + { + for (stnj=stni+1 ; stnjregy[n] = str->regwt[n] = 0.0; + n++; + } + } + + } /* end for stni */ + + return (!ok); +} + + +int xv_predict_prcp_noint(char* listgood, daymet_prcp_struct* str) +/* generates predicted prcp, with and without elevation correction +assumes the intercept from regression is 0.0 +Uses a list of good data indicators for cross-validation treatment of +missing data */ +/* PET (10/18/2015): added always_good option since this now is called with good data always */ +{ + int ok=1; + int nstns; + int stni; + double p1, pc, p, sumwt; + double* prcp; + double mx,my,mz,x,y,z; + double f, f_max; + double max_stn_day_prcp; + double max_factor = 2.0; + char *good; + int always_good = 1; + + mx = str->coef[0]; + my = str->coef[1]; + mz = str->coef[2]; + x = str->mapx; + y = str->mapy; + z = str->elev; + f_max = str->f_max; + nstns = str->count; + prcp = str->listobs+str->dam_off; + good = listgood+str->dam_off; + pc = p = sumwt = 0.0; + + max_stn_day_prcp = 0.0; + for (stni=0 ; stniwt[stni]*(p1 + mx*(x-str->listx[stni]) + + my*(y-str->listy[stni]) + mz*(z-str->listz[stni])); + sumwt += str->wt[stni]; + /* + f = mx*(x-str->listx[stni]) + my*(y-str->listy[stni]) + mz*(z - str->listz[stni]); + if (f > f_max) f = f_max; + if (f < -f_max) f = -f_max; + pc += str->wt[stni] * p1 * (1.0 + f) / (1.0 - f); + sumwt += str->wt[stni]; + */ + + if (prcp[stni] > max_stn_day_prcp) max_stn_day_prcp = prcp[stni]; + } + } + + if (sumwt) + { + pc /= sumwt; + } + else + { + printf("predict_prcp() sumwt = 0.0 ... exiting\n"); + ok=0; + } + + /* set some constraints on daily precip */ + /* PET (10/18/2015): reset the upper limit from 20 cm/day to 40 cm/day */ + if (pc < 0.0) pc = 0.0; +// if (pc > 40.0) pc = 40.0; + + // constrain to range of obs values, times a fixed factor + if (pc > max_stn_day_prcp * max_factor) pc = max_stn_day_prcp * max_factor; + + str->pp = pc; + + return (!ok); +} + + +/* xv_prcp_boxcar_smooth() generates smoothed values +from non-zero observations in the smoothing window, and returns a value which +reflects the average value per non-zero element. +For x-validation code, this function includes a check for missing data, and +requires an additional array as input, consisting of character data, +1 for good values, 0 for missing values, corresponding to the actual input +data value array. +10/25/2015, PET: Switching to always_good treatment, since for V3 code the +arrays are all filled with good values before xval is called. +*/ + +int xv_prcp_boxcar_smooth(char* good, double* input,double* output,int n,int w,int w_flag); +int xv_prcp_boxcar_smooth(char* good, double* input,double* output,int n,int w,int w_flag) +{ + int ok=1; + int tail,i,j; + int *wt; + double p,total,sum; + int always_good = 1; + + if (ok && (w > n/2)) + { + printf("Boxcar window longer than 1/2 array length...\n"); + printf("Resize window and try again\n"); + ok=0; + } + + /* establish the lengths of the boxcar tails */ + if (ok) + { + if (!(w % 2)) + w += 1; + tail = w/2; + } + + if (ok && (!(wt = (int*) malloc(w * sizeof(int))))) + { + printf("Allocation error in prcp_boxcar... Exiting\n"); + ok=0; + } + + /* when w_flag != 0, use linear ramp to weight tails, + otherwise use constant weight */ + if (ok) + { + if (w_flag) + { + for (i=0 ; i= tail) && (i < n-tail)) + { + for (j=0 ; j= n-tail) + { + for (j=0 ; j n/2)) + { + printf("Boxcar window longer than 1/2 array length...\n"); + printf("Resize window and try again\n"); + ok=0; + } + + /* establish the lengths of the boxcar tails */ + if (ok) + { + if (!(w % 2)) + w += 1; + tail = w/2; + } + + if (ok && (!(wt = (int*) malloc(w * sizeof(int))))) + { + printf("Allocation error in prcp_boxcar... Exiting\n"); + ok=0; + } + + /* when w_flag != 0, use linear ramp to weight tails, + otherwise use constant weight */ + if (ok) + { + if (w_flag) + { + for (i=0 ; i= tail) && (i < n-tail)) + { + for (j=0 ; j= n-tail) + { + for (j=0 ; j +#endif + +#include "metsrc2.h" + +/* xv_boxcar_smooth() performs a windowed smoothing on the input array, returns +result in output array. Both arrays must be doubles. n=array length, +w = windowing width, w_flag (0=flat boxcar, 1=ramped boxcar, e.g. [1 2 3 2 1]) +For x-validation code, this function includes a check for missing data, and +requires an additional array as input, consisting of character data, +1 for good values, 0 for missing values, corresponding to the actual input +data value array. + +Updated: +9/2/00, PET: added support for three-dimensional regression. +*/ +int xv_boxcar_smooth(char* good, double* input, double* output, int n, int w, int w_flag); +int xv_boxcar_smooth(char* good, double* input, double* output, int n, int w, int w_flag) +{ + int ok=1; + int tail,i,j; + int* wt; + double total,sum; + int always_good = 1; + + if (ok && (w > n/2)) + { + printf("Boxcar window longer than 1/2 array length...\n"); + printf("Resize window and try again\n"); + ok=0; + } + + /* establish the lengths of the boxcar tails */ + if (ok) + { + if (!(w % 2)) + w += 1; + tail = w/2; + } + + if (ok && (!(wt = (int*) malloc(w * sizeof(int))))) + { + printf("Allocation error in tair_boxcar... Exiting\n"); + ok=0; + } + + /* when w_flag != 0, use linear ramp to weight tails, + otherwise use constant weight */ + if (ok) + { + if (w_flag) + { + for (i=0 ; i= tail) && (i < n-tail)) + { + for (j=0 ; j= n-tail) + { + for (j=0 ; j n/2)) + { + printf("Boxcar window longer than 1/2 array length...\n"); + printf("Resize window and try again\n"); + ok=0; + } + + /* establish the lengths of the boxcar tails */ + if (ok) + { + if (!(w % 2)) + w += 1; + tail = w/2; + } + + if (ok && (!(wt = (int*) malloc(w * sizeof(int))))) + { + printf("Allocation error in tair_boxcar... Exiting\n"); + ok=0; + } + + /* when w_flag != 0, use linear ramp to weight tails, + otherwise use constant weight */ + if (ok) + { + if (w_flag) + { + for (i=0 ; i= tail) && (i < n-tail)) + { + for (j=0 ; j= n-tail) + { + for (j=0 ; jnstns; + sqdist = str->sqdist; + sqrad = str->sqrad; + gsp = str->gsp; + trunc = str->trunc; + inv_sqrad = 1.0/sqrad; + + count = 0; + wtsum = 0.0; + + for (i=0 ; iid[count] = i; + str->wt[count] = exp(-gsp * sqd * inv_sqrad) - trunc; + wtsum += str->wt[count]; + count++; + } + } + + /* normalize weights */ + if (count > 1) + { + for (i=0 ; iwt[i] /= wtsum; + } + str->count = count; + } + else + { + printf("less than 2 stations in list\n"); + ok = 0; + } + + return (!ok); +} + +int xv_fill_list(char* stmgood, char* listgood, daymet_tair_struct* str); +int xv_fill_list(char* stmgood, char* listgood, daymet_tair_struct* str) +/* fills the station-major list arrays (station-major -> day-major) +includes an array of good value flags for the cross-validation routine */ +{ + int ok=1; + int i,j; + int ndays,nstns; + long int stm_off; + + ndays = str->ndays; + nstns = str->count; + + for (i=0 ; iid[i] * ndays; + for (j=0 ; jlistobs[j*nstns + i] = str->stnobs[stm_off + j]; + str->listsmobs[j*nstns + i] = str->stnsmobs[stm_off + j]; + listgood[j*nstns + i] = stmgood[stm_off + j]; + } + str->listx[i] = str->stnx[str->id[i]]; + str->listy[i] = str->stny[str->id[i]]; + str->listz[i] = str->stnz[str->id[i]]; + } + + return (!ok); +} + +int xv_tair_regr_xwt_2switch(daymet_tair_struct* str); +int xv_tair_regr_xwt_2switch(daymet_tair_struct* str) +/* make time invariant regression arrays */ +/* generates the x and wt variables once that are constant between days */ +/* this version employs a switching mechanism to keep differencing unbiased */ +/* this version uses a double switching mechanism, alternates initial sign */ +{ + int ok=1; + int nstns; + int sign_init, sign; + int stni, stnj, n=0; + double xi,yi,elevi,wti; + + nstns = str->count; + sign_init = str->switch_init; + + for (stni=0 ; stnilistx[stni]; + yi = str->listy[stni]; + elevi = str->listz[stni]; + wti = str->wt[stni]; + for (stnj=stni+1 ; stnjregx[0][n] = sign * (xi - str->listx[stnj]); + str->regx[1][n] = sign * (yi - str->listy[stnj]); + str->regx[2][n] = sign * (elevi - str->listz[stnj]); + str->regwtall[n] = wti * str->wt[stnj]; + n++; + sign = -sign; + } + } + + return (!ok); +} + +int xv_tair_regr_y_2switch(char* listgood, daymet_tair_struct* str); +int xv_tair_regr_y_2switch(char* listgood, daymet_tair_struct* str) +/* make time variant regression array */ +/* generates the y variable for daily regressions */ +/* this version employs a switching mechanism to keep differencing unbiased */ +/* this version uses a double switching mechanism, alternates initial sign */ +/* 9/20/2015, PET: replacing the check of good with always good, since this is called + after the fill_missing_values call. +*/ +{ + int ok=1; + int nstns; + int sign_init, sign; + int stni, stnj, n=0; + double* tair; + double tair1, tair2; + char* good; + int always_good = 1; + + nstns = str->count; + sign_init = str->switch_init; + tair = str->listsmobs+str->dam_off; + good = listgood+str->dam_off; + for (stni=0 ; stniregy[n] = sign * (tair1-tair2); + str->regwt[n] = str->regwtall[n]; + } + else /* missing data */ + { + str->regy[n] = str->regwt[n] = 0.0; + } + n++; + sign = -sign; + + } /* end for stnj */ + + } /* end if (tair1) */ + + else /* missing data */ + { + for (stnj=stni+1 ; stnjregy[n] = str->regwt[n] = 0.0; + n++; + } + } + + } /* end for stni */ + + return (!ok); +} + +int xv_tair_mlr(daymet_tair_struct* str); +int xv_tair_mlr(daymet_tair_struct* str) +{ + int ok=1; + + double **x; + double *y; + double *w; + int npts; + double *coef; + double x1m,x2m,x3m,ym,wm; + double sx1,sx2,sx3,ssx1,ssx2,ssx3,ssy; + double sx1x2,sx1x3,sx2x3; + double sx1y,sx2y,sx3y; + + double dx1,dx2,dx3,dy; + double n,wt; + + double r[3]; + double a,b,c,aa,bb,cc,det; + double inv[3][3]; + double sigx[3], sigy; + double t1; + + int i,j; + + npts = str->nregr; + x=str->regx; + y=str->regy; + w=str->regwt; + coef=str->coef; + + n = (double)npts; + x1m=x2m=x3m=ym=wm=0.0; + sx1=sx2=sx3=ssx1=ssx2=ssx3=ssy=0.0; + sx1x2=sx1x3=sx2x3=0.0; + sx1y=sx2y=sx3y=0.0; + + /* step 1 : calculate means for x1,x2,x3,y,w */ + for (i=0 ; i 0.001) coef[2] = 0.001; + if (coef[2] < -0.012) coef[2] = -0.012; + + return(!ok); +} + + +int xv_predict_tair_noint(char* listgood, daymet_tair_struct* str); +int xv_predict_tair_noint(char* listgood, daymet_tair_struct* str) +/* generates predicted tair, with and without elevation correction +no regression intercept (assumed = 0.0) +Uses a list of good data indicators for cross-validation treatment of +missing data */ +/* 9/20/2015, PET: Modified to consider all station data to be good, since it has been filled. */ +{ + int ok=1; + int nstns; + int stni; + double tc, t, sumwt; + double* tair; + double mx,my,mz,x,y,z; + char *good; + int always_good = 1; + double max_stn_day_tair; + double tair_factor = 10.0; + + mx = str->coef[0]; + my = str->coef[1]; + mz = str->coef[2]; + x = str->mapx; + y = str->mapy; + z = str->elev; + nstns = str->count; + tair = str->listobs+str->dam_off; + good = listgood+str->dam_off; + tc = t = sumwt = 0.0; + + max_stn_day_tair = 0.0; + for (stni=0 ; stniwt[stni]*(tair[stni] + mx*(x-str->listx[stni]) + + my*(y-str->listy[stni]) + mz*(z-str->listz[stni])); + sumwt += str->wt[stni]; + + if (tair[stni] > max_stn_day_tair) max_stn_day_tair = tair[stni]; + } + } + + tc /= sumwt; + if (tc > max_stn_day_tair + tair_factor) tc = max_stn_day_tair + tair_factor; + + str->ta = tc; + if (str->noz) + { + t /= sumwt; + if (t > max_stn_day_tair + tair_factor) t = max_stn_day_tair + tair_factor; + str->tanoz = t; + } + + return (!ok); +}