main.c (9dd6506) - plDarkSky · Pawel Netzel

main.c

/****************************************************************************
 *
 * PROGREM:      darksky
 *
 * AUTHOR(S):    Pawel Netzel - SIL UC, Henryka Netzel - CAMK
 *
 * PURPOSE:      The program calculates sky lightness in the zenith direction
 *               GHSL data are taken as input (values from 0 to 1).
 *               If DEM is provided, the program takes into account shadowing
 *               effect. The unit of the results is magnitudo.
 *
 *               If the program is used, cite following:
 *
 *               H. Netzel, P. Netzel, 2016: High resolution map of light 
 *               pollution over Poland. J. of Quantitative Spectroscopy and 
 *               Radiative Transfer, Vol. 181, pp. 67-73
 *
 *
 * COPYRIGHT:    (C) 2016 by Pawel Netzel
 *
 *               This program is free software under the GNU General Public
 *               License (>=v2): 
 *                 http://www.gnu.org/licenses/old-licenses/gpl-2.0.en.html
 *
 ***************************************************************************/

#include <stdlib.h>
#include <string.h>
#include <ezgdal.h>
#include <math.h>
#include <omp.h>

#include "argtable3.h"
#include "version.h"
#include "viewshed.h"

/*******************************
 *     Calc matrix with model
 *     weights
 */
double **init_model_matrix(double par_a, double par_U, double par_V, double par_h,double par_k, 
                           int par_N, double resX, double resY) {
  double **mtx;
  int r,c,n;
  double range,px,py,D,t1,t2;
  
  n = par_N*2+1;
  range = resX*par_N;
  
  mtx = malloc(n*sizeof(double *));
  for(r=0; r<n; r++)
    mtx[r] = malloc(n*sizeof(double));

for(r=0; r<n; r++)
    for(c=0; c<n; c++) {
      px = resX*(c-par_N);
      py = resY*(r-par_N);
      D=(px*px+py*py);
      if(sqrt(D)<=range) {
        t1=D+par_h*par_h;
        t2=sqrt(t1);
        mtx[r][c]=par_a*(par_U/t1+par_V/t2)*exp(-par_k*t2);
      } else
        mtx[r][c]=0.0;
    }
    
  return mtx;
}

void free_model_matrix(double **mtx, int par_N) {
  int r, n;
  n = 2*par_N+1;
  for(r=0; r<n; r++)
    free(mtx[r]);
  free(mtx);
}

/**********************************
 *    Calc model for a cell
 *    Multiply model matrix by 
 *    GHSL values in the cell 
 *    perimeter
 */
double calc(EZGDAL_FRAME *f, double **mtx, int **visible) {
  int r, n;
  double sum=0.0;
  EZGDAL_LAYER *l = f->owner.stripe->layer;
  
  n = f->cols;
  
#pragma omp parallel for shared(n,l,f,visible,mtx) reduction( + : sum )
  for(r=0; r<n; r++) {
    double *pm = mtx[r];
    double *pb = f->buffer[r];
    int *vis = visible[r];
    double ss = 0.0;
    int c;
    for(c=0; c<n; c++) {
      if(!ezgdal_is_null(l,*pb) && *vis==1)
        ss += (*pb) * (*pm);
      pm++; pb++; vis++;
    }
    sum = sum + ss;
  }
  
  double v = -2.5*log10(sum+250.0)+27.78;

return v;
}

int quiet_val = 0;

/**********************************
 *     Calculation sky brightness
 */
void calc_model(EZGDAL_LAYER *ghsl_layer, EZGDAL_LAYER *dem_layer, EZGDAL_LAYER *out_layer, 
                double **model_matrix, double **distance, int **visibility, 
                double height0, double radius, double scaleH) {
  int r,c,c1,c2,r1,r2,rN;
  int prc = 0;
  
  rN = ghsl_layer->rows;
  c1 = ghsl_layer->stripe->frame->cols/2;
  c2 = ghsl_layer->cols - c1 + 1;
  r1 = ghsl_layer->stripe->rows/2;
  r2 = ghsl_layer->rows - r1 + 1;
  
  if(quiet_val!=1) ezgdal_show_progress(stderr,0,rN);

for(c=0; c<out_layer->cols; c++)
    out_layer->buffer[c] = out_layer->no_data;
  for(r=0; r<=r1; r++) {
    if(quiet_val!=1) ezgdal_show_progress(stderr,r,rN);
    ezgdal_write_buffer(out_layer,r);
  }

for(r=r1+1; r<r2; r++) {
    if(quiet_val!=1) ezgdal_show_progress(stderr,r,rN);
    
    ezgdal_load_stripe_data(ghsl_layer->stripe, r-r1);
    if(dem_layer!=NULL) ezgdal_load_stripe_data(dem_layer->stripe, r-r1);

for(c=c1; c<c2; c++) {
      ezgdal_shift_frame_pos(ghsl_layer->stripe->frame, c-c1);
      if(dem_layer!=NULL) {
        ezgdal_shift_frame_pos(dem_layer->stripe->frame, c-c1);
        calc_visibility_mask(dem_layer->stripe->frame, distance, visibility, height0, scaleH, radius, dem_layer->stripe->frame->cols);
      }
      out_layer->buffer[c] = calc(ghsl_layer->stripe->frame, model_matrix, visibility);
    }

ezgdal_write_buffer(out_layer,r);

}

for(c=0; c<out_layer->cols; c++)
    out_layer->buffer[c] = out_layer->no_data;
  for(r=r2; r<rN; r++) {
    if(quiet_val!=1) ezgdal_show_progress(stderr,r,rN);
    if(100*r/rN != prc) {prc=100*r/rN; printf("\b\b\b\b%3d%%",prc); fflush(stdout); }
    ezgdal_write_buffer(out_layer,r);
  }
  if(quiet_val!=1) ezgdal_show_progress(stderr,rN,rN);
}

/*************************************
 *     Check the file
 */
int file_exists(const char *fname) {
  FILE *f = fopen(fname,"r");
  if(f) {
    fclose(f);
    return 1;
  }
  return 0;
}

/*************************************
 *     Display program usage and exit
 */
void usage(char *progname, void *argtable) {
      printf("\nSoftware version: %s, build date: %s\n",__RELEASE_VERSION__,__DATE__);
      printf("\nUsage:\n\t%s", progname);
      arg_print_syntax(stdout,argtable,"\n");
      printf("\n");
      arg_print_glossary_gnu(stdout,argtable);
      printf("\n");
      exit(0);
}

int main(int argc, char *argv[])
{
  double par_a  = 2.37583;
  double par_U  = 5.96033;
  double par_V  = 0.02625;
  double par_h  = 1.13661;
  double par_k  = 0.03875;
  int    par_N  = 300;
  double resX   = 0.1;
  double resY   = 0.1;
  double resH   = 0.001;

struct arg_str  *ghsl    = arg_str1("g","ghsl","<file_name>","name of GHSL file (GeoTIFF)");
  struct arg_str  *dem     = arg_str0("d","dem","<file_name>","name of DEM file (GeoTIFF)");
  struct arg_str  *out     = arg_str1("o","output","<file_name>","name of output file (GeoTIFF)");
  struct arg_dbl  *p_par_a = arg_dbl0("a","a","<dbl>","parameter 'a' in Berry's model (default: 2.37583)");
  struct arg_dbl  *p_par_U = arg_dbl0("U","U","<dbl>","parameter 'U' in Berry's model (default: 5.96033)");
  struct arg_dbl  *p_par_V = arg_dbl0("V","V","<dbl>","parameter 'V' in Berry's model (default: 0.02625)");
  struct arg_dbl  *p_par_h = arg_dbl0("h","h","<dbl>","parameter 'h' in Berry's model (default: 1.13661)");
  struct arg_dbl  *p_par_k = arg_dbl0("k","k","<dbl>","parameter 'k' in Berry's model (default: 0.03875)");
  struct arg_int  *p_par_N = arg_int0("N","radius","<n>","model radius in cells (default: 300)");
  struct arg_dbl  *p_resX  = arg_dbl0(NULL,"res_x","<dbl>","scale factor: column to km (default: 0.1)");
  struct arg_dbl  *p_resY  = arg_dbl0(NULL,"res_y","<dbl>","scale factor: row to km (default: 0.1)");
  struct arg_dbl  *p_resH  = arg_dbl0(NULL,"res_h","<dbl>","scale factor: height to km (default: 0.001)");
  struct arg_int  *th      = arg_int0("t","threads","<n>","number of threads (default 1)");
  struct arg_lit  *quiet   = arg_litn("q","quiet",0,2,"print help and exit");
  struct arg_lit  *help    = arg_lit0(NULL,"help","print help and exit");
  struct arg_end  *end     = arg_end(20);
  void* argtable[] = {ghsl,dem,out,p_par_a,p_par_U,p_par_V,p_par_h,p_par_k,p_par_N,p_resX,p_resY,p_resH,th,quiet,help,end};

int nerrors = arg_parse(argc,argv,argtable);

if(help->count > 0 ) 
    usage(argv[0],argtable);

if(nerrors != 0) 
    usage(argv[0],argtable);

quiet_val = quiet->count;

if(!file_exists(ghsl->sval[0])) {
    printf("\nInput file [%s] does not exists!\n\n",ghsl->sval[0]);
    exit(1);
  }

if(dem->count>0 && !file_exists(dem->sval[0])) {
    printf("\nInput file [%s] does not exists!\n\n",dem->sval[0]);
    exit(1);
  }

if(p_par_a->count>0) par_a = p_par_a->dval[0];
  if(p_par_U->count>0) par_U = p_par_U->dval[0];
  if(p_par_V->count>0) par_V = p_par_V->dval[0];
  if(p_par_h->count>0) par_h = p_par_h->dval[0];
  if(p_par_k->count>0) par_k = p_par_k->dval[0];
  if(p_par_N->count>0) par_N = p_par_N->ival[0];

if(p_resX->count>0) resX = p_resX->dval[0];
  if(p_resY->count>0) resY = p_resY->dval[0];

if(th->count>0 && th->ival[0]>0 && th->ival[0]<32) 
    omp_set_num_threads(th->ival[0]);
  else 
    omp_set_num_threads(1);

if(quiet_val==0) {
    printf("\nInput file (GHSL): %s\n",ghsl->sval[0]);
    if(dem->count>0)
      printf("Input file (DEM): %s\n",dem->sval[0]);
    printf("Output file (sky brightness): %s\n",out->sval[0]);
    printf("Berry's model parameters: a=%lf, U=%lf, V=%lf, h=%lf, k=%lf\n",par_a,par_U,par_V,par_h,par_k);
    printf("                  radius: %d\n\n",par_N);
  }

int no_data;
  EZGDAL_LAYER *ghsl_layer = NULL;
  EZGDAL_LAYER *dem_layer = NULL;
  EZGDAL_LAYER *out_layer = NULL;
  double **model_matrix;
  double *at;
  double **distance;
  int **visibility;

if(quiet_val==0) ezgdal_show_message_nolf(stdout, " - preparing environment ... ");
    no_data = -9999;

if(dem->count>0) {
      dem_layer = ezgdal_open_layer((char *)dem->sval[0]);
      if(dem_layer==NULL) exit(1);
      dem_layer->stripe = ezgdal_create_stripe(dem_layer,0,2*par_N+1);
      dem_layer->stripe->frame = ezgdal_create_frame(dem_layer->stripe,0);
    }

ghsl_layer = ezgdal_open_layer((char *)ghsl->sval[0]);
    if(ghsl_layer==NULL) exit(1);
    
    ghsl_layer->stripe = ezgdal_create_stripe(ghsl_layer,0,2*par_N+1);
    ghsl_layer->stripe->frame = ezgdal_create_frame(ghsl_layer->stripe,0);
    at = ezgdal_layer_get_at(ghsl_layer);

out_layer = ezgdal_create_layer((char *)(out->sval[0]),
                                  ezgdal_layer_get_wkt(ghsl_layer),
                                  "Float32",
                                  at,
                                  ghsl_layer->rows,
                                  ghsl_layer->cols,
                                  &no_data,
                                  EZGDAL_COMPRESS_DEFLATE);
    if(out_layer==NULL) exit(1);
    
    if(p_resX->count==0) {
      if(at[2]!=0.0) {
        ezgdal_show_message(stderr,"\n\nYou have to provide 'res_x' parameter\n\n");
        exit(1);
      } else 
        resX = 0.001*fabs(at[1]);
    }
    if(p_resY->count==0) {
      if(at[4]!=0.0) {
        ezgdal_show_message(stderr,"\n\nYou have to provide 'res_y' parameter\n\n");
        exit(1);
      } else 
        resY = 0.001*fabs(at[5]);
    }
    if(p_resH->count>0) 
      resH = p_resH->dval[0];

model_matrix = init_model_matrix(par_a,par_U,par_V,par_h,par_k,par_N,resX,resY);
    distance = create_dist(2*par_N+1,resX);
    visibility = create_mask(2*par_N+1,distance,par_N*resX);
  if(quiet_val==0) ezgdal_show_message(stdout,"OK");

if(quiet_val==0) ezgdal_show_message_nolf(stdout," - calculating ... ");
  
    calc_model(ghsl_layer, dem_layer, out_layer, 
               model_matrix, distance, visibility, 
               par_h, par_N*resX, resH);

if(quiet_val==0) ezgdal_show_message(stdout," ... OK");

if(quiet_val==0) ezgdal_show_message_nolf(stdout," - closing environment ... ");
    ezgdal_close_layer(ghsl_layer);
    if(dem->count>0)
      ezgdal_close_layer(dem_layer);
    ezgdal_close_layer(out_layer);
    free_model_matrix(model_matrix, par_N);
    free_array((void **)distance, 2*par_N+1);
    free_array((void **)visibility, 2*par_N+1);
  if(quiet_val==0) ezgdal_show_message(stdout,"OK");

arg_freetable(argtable, sizeof(argtable) / sizeof(argtable[0]));

return 0;
}

plDarkSky

Changes in progress indicator.