Skip to content

Commit 37bb856

Browse files
davidfrantzdavidfrantz
committed
Merge branch 's2_angles' into develop
2 parents 6ccf5ba + 5951278 commit 37bb856

File tree

2 files changed

+167
-68
lines changed

2 files changed

+167
-68
lines changed

‎src/lower-level/meta-ll.c‎

Lines changed: 165 additions & 68 deletions
Original file line numberDiff line numberDiff line change
@@ -653,9 +653,9 @@ int i = 0, j = 0, ii, jj;
653653
int left, right;
654654
int top, bottom;
655655
int sv_nx = 0, sv_ny = 0;
656-
float **s_sz = NULL, **s_sa = NULL;
657-
float **s_vz = NULL, **s_va = NULL;
658-
int **k_vz = NULL, **k_va = NULL;
656+
int d = 0, nd;
657+
float *s_sz = NULL, *s_sa = NULL;
658+
float ***s_vz = NULL, ***s_va = NULL;
659659
bool s = false, v = false, z = false, a = false;
660660
char d_img[NPOW_10];
661661
char id_img[NPOW_10];
@@ -677,7 +677,8 @@ int svgrid = 5000;
677677
meta->sat = 65535;
678678

679679

680-
nb = 13;
680+
nb = 13; // number of bands
681+
nd = 12; // number of detectors
681682
DN = allocate_brick(nb, 0, _DT_NONE_);
682683

683684
set_brick_res(DN, INT_MAX);
@@ -961,20 +962,33 @@ int svgrid = 5000;
961962
if (svgrid != atoi(tokenptr)){
962963
printf("SUN_VIEW_GRID is incompatible with Sentinel-2 metadata. "); return FAILURE;}
963964

964-
sv_nx = ceil(get_brick_width(DN)/(float)svgrid);
965-
sv_ny = ceil(get_brick_height(DN)/(float)svgrid);
966-
if (s_sz == NULL) alloc_2D((void***)&s_sz, sv_ny, sv_nx, sizeof(float));
967-
if (s_vz == NULL) alloc_2D((void***)&s_vz, sv_ny, sv_nx, sizeof(float));
968-
if (s_sa == NULL) alloc_2D((void***)&s_sa, sv_ny, sv_nx, sizeof(float));
969-
if (s_va == NULL) alloc_2D((void***)&s_va, sv_ny, sv_nx, sizeof(float));
970-
if (k_vz == NULL) alloc_2D((void***)&k_vz, sv_ny, sv_nx, sizeof(int));
971-
if (k_va == NULL) alloc_2D((void***)&k_va, sv_ny, sv_nx, sizeof(int));
965+
sv_nx = ceil(get_brick_width(DN)/(float)svgrid) + 1;
966+
sv_ny = ceil(get_brick_height(DN)/(float)svgrid) + 1;
967+
if (s_sz == NULL) alloc((void**)&s_sz, sv_ny*sv_nx, sizeof(float));
968+
if (s_sa == NULL) alloc((void**)&s_sa, sv_ny*sv_nx, sizeof(float));
969+
if (s_vz == NULL) alloc_3D((void****)&s_vz, nb, nd, sv_ny*sv_nx, sizeof(float));
970+
if (s_va == NULL) alloc_3D((void****)&s_va, nb, nd, sv_ny*sv_nx, sizeof(float));
971+
//if (k_vz == NULL) alloc_2D((void***)&k_vz, sv_ny, sv_nx, sizeof(int));
972+
//if (k_va == NULL) alloc_2D((void***)&k_va, sv_ny, sv_nx, sizeof(int));
972973
} else if (strcmp(tag, "Sun_Angles_Grid") == 0){
973974
s = true;
974975
} else if (strcmp(tag, "/Sun_Angles_Grid") == 0){
975976
s = false;
976977
} else if (strcmp(tag, "Viewing_Incidence_Angles_Grids") == 0){
977978
v = true;
979+
b = 0;
980+
d = 0;
981+
while (tokenptr != NULL){
982+
if (strcmp(tokenptr, "bandId") == 0){
983+
tokenptr = strtok(NULL, separator);
984+
b = atoi(tokenptr);
985+
}
986+
if (strcmp(tokenptr, "detectorId") == 0){
987+
tokenptr = strtok(NULL, separator);
988+
d = atoi(tokenptr) - 1;
989+
}
990+
tokenptr = strtok(NULL, separator);
991+
}
978992
} else if (strcmp(tag, "/Viewing_Incidence_Angles_Grids") == 0){
979993
v = false;
980994
} else if (strcmp(tag, "Zenith") == 0){
@@ -993,20 +1007,18 @@ int svgrid = 5000;
9931007

9941008
if (s){
9951009
if (z){
996-
if (strcmp(tokenptr, "NaN") !=0) s_sz[i][j] = atof(tokenptr);
1010+
if (strcmp(tokenptr, "NaN") !=0) s_sz[i*sv_ny+j] = atof(tokenptr);
9971011
} else if (a){
998-
if (strcmp(tokenptr, "NaN") !=0) s_sa[i][j] = atof(tokenptr);
1012+
if (strcmp(tokenptr, "NaN") !=0) s_sa[i*sv_ny+j] = atof(tokenptr);
9991013
}
10001014
} else if (v){
10011015
if (z){
10021016
if (strcmp(tokenptr, "NaN") !=0){
1003-
s_vz[i][j] += atof(tokenptr);
1004-
k_vz[i][j]++;
1017+
s_vz[b][d][i*sv_ny+j] += atof(tokenptr);
10051018
}
10061019
} else if (a){
10071020
if (strcmp(tokenptr, "NaN") !=0){
1008-
s_va[i][j] += atof(tokenptr);
1009-
k_va[i][j]++;
1021+
s_va[b][d][i*sv_ny+j] += atof(tokenptr);
10101022
}
10111023
}
10121024
}
@@ -1029,15 +1041,33 @@ int svgrid = 5000;
10291041
fclose(fp);
10301042

10311043

1032-
// get image subset
1044+
// interpolate the grids
1045+
// angles are given at grid intersections, we need one value per cell
1046+
interpolate_sunview_grid(s_sz, sv_nx, sv_ny, meta->s2.nodata);
1047+
interpolate_sunview_grid(s_sa, sv_nx, sv_ny, meta->s2.nodata);
1048+
for (b=0; b<nb; b++){
1049+
for (d=0; d<nd; d++){
1050+
interpolate_sunview_grid(s_vz[b][d], sv_nx, sv_ny, meta->s2.nodata);
1051+
interpolate_sunview_grid(s_va[b][d], sv_nx, sv_ny, meta->s2.nodata);
1052+
}
1053+
}
1054+
sv_nx--;
1055+
sv_ny--;
1056+
1057+
// collapse view grids
1058+
collapse_view_grid(s_vz, nb, nd, sv_nx, sv_ny, meta->s2.nodata);
1059+
collapse_view_grid(s_va, nb, nd, sv_nx, sv_ny, meta->s2.nodata);
10331060

1061+
1062+
// get image subset
10341063
left = sv_nx-1; right = 0;
10351064
top = sv_ny-1; bottom = 0;
10361065

10371066
for (i=0; i<sv_ny; i++){
10381067
for (j=0; j<sv_nx; j++){
10391068

1040-
if (k_vz[i][j] > 0 && k_va[i][j] > 0){
1069+
if (s_vz[0][0][i*sv_ny+j] != meta->s2.nodata &&
1070+
s_va[0][0][i*sv_ny+j] != meta->s2.nodata){
10411071
if (j < left) left = j;
10421072
if (j > right) right = j;
10431073
if (i < top) top = i;
@@ -1081,73 +1111,35 @@ int svgrid = 5000;
10811111
if (meta->s2.vazi == NULL) alloc_2D((void***)&meta->s2.vazi, meta->s2.ny, meta->s2.nx, sizeof(float));
10821112

10831113

1084-
// average of view angles
1114+
// copy values to final view sun grids
10851115
for (i=0; i<meta->s2.ny; i++){
10861116
for (j=0; j<meta->s2.nx; j++){
10871117

10881118
ii = i+top;
10891119
jj = j+left;
10901120

1091-
if (k_vz[ii][jj] > 0){
1092-
meta->s2.vzen[i][j] = s_vz[ii][jj]/k_vz[ii][jj];
1093-
meta->s2.szen[i][j] = s_sz[ii][jj];
1094-
} else {
1095-
meta->s2.vzen[i][j] = meta->s2.nodata;
1096-
meta->s2.szen[i][j] = meta->s2.nodata;
1097-
}
1098-
if (k_va[ii][jj] > 0){
1099-
meta->s2.vazi[i][j] = s_va[ii][jj]/k_va[ii][jj];
1100-
meta->s2.sazi[i][j] = s_sa[ii][jj];
1101-
} else {
1102-
meta->s2.vazi[i][j] = meta->s2.nodata;
1103-
meta->s2.sazi[i][j] = meta->s2.nodata;
1104-
}
1105-
1106-
}
1107-
}
1108-
1121+
meta->s2.szen[i][j] = s_sz[ii*sv_ny+jj];
1122+
meta->s2.sazi[i][j] = s_sa[ii*sv_ny+jj];
1123+
meta->s2.vzen[i][j] = s_vz[0][0][ii*sv_ny+jj];
1124+
meta->s2.vazi[i][j] = s_va[0][0][ii*sv_ny+jj];
11091125

1110-
// try to fill the left edge (duplicate values - silly method, but it will do for now)
1111-
1112-
// average of view angles
1113-
for (i=0; i<meta->s2.ny; i++){
1114-
for (j=0; j<meta->s2.nx; j++){
1115-
1116-
if ((jj = j+1) == meta->s2.nx) continue;
1117-
1118-
if (meta->s2.vzen[i][j] == meta->s2.nodata &&
1119-
meta->s2.vzen[i][jj] != meta->s2.nodata){
1120-
1121-
meta->s2.vzen[i][j] = meta->s2.vzen[i][jj];
1122-
meta->s2.szen[i][j] = meta->s2.szen[i][jj];
1123-
meta->s2.vazi[i][j] = meta->s2.vazi[i][jj];
1124-
meta->s2.sazi[i][j] = meta->s2.sazi[i][jj];
1125-
1126-
}
11271126

11281127
}
11291128
}
11301129

11311130

1132-
1133-
1134-
//right++;
1135-
//bottom++;
1136-
11371131
meta->s2.left = left * svgrid/get_brick_res(DN);
11381132
meta->s2.top = top * svgrid/get_brick_res(DN);
11391133
meta->s2.right = right * svgrid/get_brick_res(DN);
11401134
meta->s2.bottom = bottom * svgrid/get_brick_res(DN);
11411135

1142-
if (meta->s2.right > get_brick_ncols(DN)) meta->s2.right = get_brick_ncols(DN);
1136+
if (meta->s2.right > get_brick_ncols(DN)) meta->s2.right = get_brick_ncols(DN);
11431137
if (meta->s2.bottom > get_brick_nrows(DN)) meta->s2.bottom = get_brick_nrows(DN);
11441138

1145-
free_2D((void**)s_sz, sv_ny);
1146-
free_2D((void**)s_sa, sv_ny);
1147-
free_2D((void**)s_vz, sv_ny);
1148-
free_2D((void**)s_va, sv_ny);
1149-
free_2D((void**)k_vz, sv_ny);
1150-
free_2D((void**)k_va, sv_ny);
1139+
free((void*)s_sz);
1140+
free((void*)s_sa);
1141+
free_3D((void***)s_vz, nb, nd);
1142+
free_3D((void***)s_va, nb, nd);
11511143

11521144
#ifdef FORCE_DEBUG
11531145
printf("active image subset: UL %d/%d to LR %d/%d\n", meta->s2.left, meta->s2.top, meta->s2.right, meta->s2.bottom);
@@ -1333,6 +1325,111 @@ int nchar;
13331325
}
13341326

13351327

1328+
/** This function interpolates the sun and view angle grids given in the
1329+
+++ Sentinel-2 metadata. The interpolation grid will be interpolated to
1330+
+++ a cell-based grid, and then clipped to the image extent.
1331+
+++ int_grid: interpolation grid (modified)
1332+
--- nx: number of columns
1333+
--- ny: number of rows
1334+
--- nodata: nodata value
1335+
+++ Return: void
1336+
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++**/
1337+
void interpolate_sunview_grid(float *int_grid, int int_nx, int int_ny, float nodata){
1338+
int i, j, ii, jj;
1339+
int cell_nx, cell_ny;
1340+
int cell_p, int_p;
1341+
float *cell_grid = NULL;
1342+
float sum, num;
1343+
1344+
1345+
cell_nx = int_nx - 1;
1346+
cell_ny = int_ny - 1;
1347+
1348+
alloc((void**)&cell_grid, cell_ny*cell_nx, sizeof(float));
1349+
1350+
1351+
for (i=0, cell_p=0; i<cell_ny; i++){
1352+
for (j=0; j<cell_nx; j++, cell_p++){
1353+
1354+
sum = num = 0;
1355+
1356+
for (ii=0; ii<=1; ii++){
1357+
for (jj=0; jj<=1; jj++){
1358+
1359+
int_p = (i+ii) * int_ny + (j+jj);
1360+
1361+
if (int_grid[int_p] > 0){
1362+
sum += int_grid[int_p];
1363+
num++;
1364+
}
1365+
}
1366+
}
1367+
1368+
if (num > 0){
1369+
cell_grid[cell_p] = sum/num;
1370+
} else {
1371+
cell_grid[cell_p] = nodata;
1372+
}
1373+
1374+
}
1375+
}
1376+
1377+
memmove(int_grid, cell_grid, cell_nx*cell_ny*sizeof(float));
1378+
re_alloc((void**)&int_grid, int_nx*int_ny, cell_nx*cell_ny, sizeof(float));
1379+
free((void*)cell_grid);
1380+
1381+
return;
1382+
}
1383+
1384+
1385+
/** This function collapses the Sentinel-2 view angle grids to a single
1386+
+++ grid. All band- and detector-based grids will be averaged. There is
1387+
+++ room for improvement at this point. The collapsed grid will be put
1388+
+++ into the first slot.
1389+
+++ grid: 3D grid (modified)
1390+
--- nb: number of bands
1391+
--- nd: number of detectors
1392+
--- nx: number of columns
1393+
--- ny: number of rows
1394+
--- nodata: nodata value
1395+
+++ Return: void
1396+
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++**/
1397+
void collapse_view_grid(float ***grid, int nb, int nd, int nx, int ny, float nodata){
1398+
int b, d, i, j, p;
1399+
float sum, num;
1400+
1401+
1402+
for (i=0, p=0; i<ny; i++){
1403+
for (j=0; j<nx; j++, p++){
1404+
1405+
sum = num = 0;
1406+
1407+
for (b=0; b<nb; b++){
1408+
for (d=0; d<nd; d++){
1409+
1410+
if (grid[b][d][p] != nodata){
1411+
sum += grid[b][d][p];
1412+
num++;
1413+
}
1414+
1415+
}
1416+
}
1417+
1418+
if (num > 0){
1419+
grid[0][0][p] = sum/num;
1420+
} else {
1421+
grid[0][0][p] = nodata;
1422+
}
1423+
1424+
}
1425+
}
1426+
1427+
1428+
return;
1429+
}
1430+
1431+
1432+
13361433
/** This function identifies the satellite mission, i.e. Landsat or Senti-
13371434
+++ nel-2
13381435
--- pl2: L2 parameters

‎src/lower-level/meta-ll.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -83,6 +83,8 @@ int print_metadata(meta_t *meta, int nb);
8383
int parse_metadata_landsat(par_ll_t *pl2, meta_t *meta, brick_t **dn);
8484
int parse_metadata_sentinel2(par_ll_t *pl2, meta_t *meta, brick_t **dn);
8585
void parse_metadata_band(char *d_level1, char *tag, char *value, cal_t *cal, int lid, int type);
86+
void interpolate_sunview_grid(float *int_grid, int int_nx, int int_ny, float nodata);
87+
void collapse_view_grid(float ***grid, int nb, int nd, int nx, int ny, float nodata);
8688
int parse_metadata_mission(par_ll_t *pl2);
8789
int parse_metadata(par_ll_t *pl2, meta_t **metadata, brick_t **DN);
8890

0 commit comments

Comments
 (0)