TheBoussinesqModel  3.2.1
 All Data Structures Files Functions Variables Typedefs Macros Pages
g_raster2plvector.c
Go to the documentation of this file.
1 
23 #include "turtle.h"
24 #include "t_alloc.h"
25 
26 
36 #include "linear_span.h"
37 #include "geometry.h"
38 #include "geometry_io.h"
39 #include "geometry_freememory.h"
40 #include "rw_maps.h"
41 //#include "gridded.element.input.geotop.h"
42 
43 #include "g_raster2plvector.h"
44 
45 #define NO_ELEVATION 0.0
46 #define LX_CENTER -0.5
47 #define LY_CENTER 0.5
48 #define TOP 1.0
49 #define BOTTOM 0.0
50 #define LEFT -1.0
51 #define RIGHT 0.0
52 
53 #define IVERTEX 10000
54 
55 #define INIT_VALUE -999
56 
57 #define FLOATING_POINT_TYPE 0
58 #define MAP_FORMAT 3 // output esri asccii according to Stefano Endrizzi's 2007 functions
59 //#define MAP_FORMAT 2
60 
61 
62 
63 //LONGBIN *addresses(LONGMATRIX *lmask,DOUBLEVECTOR *V, int (* check_novalue2)(double x,double y,DOUBLEVECTOR *V),int reply,long row_shift,long col_shift){
64 /*!* \param V (DOUBLEVECTOR *) - vector regarding the query appied (searching novalue)
65  * \param (* check_novalue)(double x,DOUBLEVECTOR *V) (int) - query expressed as a function returning a long integer
66  * \param reply - expected reply from the query to save the element addresses
67  */
68 
69 
70 //int create_raster_index_matrices(LONGMATRIX *lmask,long novalue, LONGMATRIX *i_pixels,LONGMATRIX *i_horizontal_lines, LONGMATRIX *i_vertical_lines, LONGMATRIX *i_vertex){
78 LONGBIN *addresses(LONGMATRIX *lmask,long row_shift,long col_shift){
92  LONGVECTOR *ndata;
93 
94  LONGBIN *LB;
95  long r,c,j,data_counter;
96 
97  printf("ENTRANCE: %ld %ld %ld %ld ",row_shift,col_shift,lmask->nrh,lmask->nch);
98 
99  if ((row_shift<-1)|| (row_shift>1) || (col_shift<-1)|| (col_shift>1) ) {
100  printf ("Warning in addresses function source file g_raster2plvector_novalue_manegement: parameters row_shift or col_shift exceed 1 (in absolute value) \n");
101  printf ("This could cause bus errors!!\n ");
102  }
103 
104 
105  for (r=lmask->nrl;r<=lmask->nrh;r++) {
106  if ((lmask->element[r][lmask->ncl]==1) || (lmask->element[r][lmask->nch]==1)) printf ("Warning in addresses function source file g_raster2plvector_novalue_manegement: at row %ld the mask of the basin reaches the border !\n",r);
107  }
108  for (c=lmask->ncl;c<=lmask->nch;c++) {
109  if ((lmask->element[lmask->nrl][c]==1) || (lmask->element[lmask->nrh][c]==1)) printf ("Warning in addresses function source file g_raster2plvector_novalue_manegement: at column %ld the mask of the basin reaches the border!\n",c);
110  }
111  ndata=new_longvector(lmask->nrh);
112  ndata->element[ndata->nl]=0;
113 // ndata->element[ndata->nh]=0;
114  for (r=lmask->nrl+1;r<=lmask->nrh;r++){
115  data_counter=0;
116  for (c=lmask->ncl+1;c<=lmask->nch;c++){
117  if ((lmask->element[r][c]==1) || (lmask->element[r][c+col_shift]==1) || (lmask->element[r+row_shift][c]==1) || (lmask->element[r+row_shift][c+col_shift]==1))
118  data_counter++;
119 
120  }
121  ndata->element[r]=data_counter;
122  }
123 // print_longvector_elements(ndata,PRINT);
124 // stop_execution();
125  LB=new_longbin(ndata);
126 // printf("OKK pl");
127 // stop_execution();
128  for (r=lmask->nrl+1;r<=lmask->nrh;r++){
129  j=0;
130  for (c=lmask->ncl+1;c<=lmask->nch;c++){
131  if ((lmask->element[r][c]==1) || (lmask->element[r][c+col_shift]==1) || (lmask->element[r+row_shift][c]==1) || (lmask->element[r+row_shift][c+col_shift]==1)) {
132  j++;
133  LB->element[r][j]=c;
134  }
135  }
136 
137  }
138  free_longvector(ndata);
139  return LB;
140 
141 }
142 //COME FARE CON I VERTICI???
143 
144 
145 
146 long index_pixel_from_a_bin(long r, long c, LONGVECTOR *s_index){
162  long l;
163 
164  if (r%2==0) {
165  l=s_index->element[r]-c+1;
166  }else {
167  if (r>1) {
168  l=s_index->element[r-1]+c;
169  } else {
170  l=c;
171  }
172 
173 
174  }
175  return l;
176 
177 }
178 
179 LONGMATRIX *m_indices_with_novalues(LONGBIN *laddresses, long nch, long novalue, long IBASE, long (*index_pixel_from_a_bin)(long r, long c, LONGVECTOR *s_index)){
193  LONGMATRIX *m_ind;
194  LONGVECTOR *s_index;
195  long r,c,j,l;
196 
197  s_index=new_longvector(laddresses->index->nh);
198  s_index->element[s_index->nl]=laddresses->index->element[s_index->nl];
199  for (l=s_index->nl+1;l<=s_index->nh;l++){
200  if (laddresses->index->element[l]>nch) printf("Error: in m_indices_with_novalues at row %ld , %ld exceeds number of columns %ld \n",l,laddresses->index->element[l],nch);
201  s_index->element[l]=s_index->element[l-1]+laddresses->index->element[l];
202  }
203 // stop_execution();
204 // printf("22novalue %ld :\n",novalue);
205 // stop_execution();
206  m_ind=new_longmatrix(laddresses->index->nh,nch);
207  for(r=m_ind->nrl;r<=m_ind->nrh;r++) {
208  for(c=m_ind->ncl;c<=m_ind->nch;c++) {
209  m_ind->element[r][c]=novalue;
210  // printf("m_IND[%ld][%ld]: %ld\n",r,c,m_ind->element[r][c]);
211 
212  }
213  for(j=m_ind->ncl;j<=laddresses->index->element[r];j++){
214  m_ind->element[r][laddresses->element[r][j]]=(*index_pixel_from_a_bin)(r,j,s_index)+IBASE;
215 // printf("m_IND[%ld][%ld]: %ld\n",r,laddresses->element[r][j],m_ind->element[r][laddresses->element[r][j]]);
216  }
217 
218 
219  }
220 
221  free_longvector(s_index);
222  return m_ind;
223 
224 
225 
226 
227 }
228 
229 
235 LONGMATRIX *m_indices_from_mask(DOUBLEMATRIX *mask,long row_shift,long col_shift,long novalue,long IBASE,long (*index_pixel_from_a_bin)(long r, long c,LONGVECTOR *s_index),DOUBLEVECTOR *V, int (*check_novalues)(double x, DOUBLEVECTOR *V)){
252  LONGBIN *LB;
253  LONGMATRIX *m_ind;
254  LONGMATRIX *lmask;
255  long r,c;
256 
257  lmask=new_longmatrix(mask->nrh,mask->nch);
258  printf("1: lmask %ld %ld:\n",lmask->nrh,lmask->nch);
259 // stop_execution();
260  for (r=mask->nrl;r<=mask->nrh;r++) {
261  for (c=mask->ncl;c<=mask->nch;c++) {
262  // printf("value:%lf at %ld,%ld \n",mask->element[r][c],r,c);
263  if ((*check_novalues)(mask->element[r][c],V)==0) {
264  // printf("NONNULL value:%lf at %ld,%ld \n",mask->element[r][c],r,c);
265  lmask->element[r][c]=1;
266  // stop_execution();
267  }else if ((*check_novalues)(mask->element[r][c],V)==1) {
268  lmask->element[r][c]=0;
269  // printf("NULL value:%lf at %ld,%ld \n",mask->element[r][c],r,c);
270  }
271  }
272  }
273  printf("lmask %ld %ld :\n",lmask->nrh,lmask->nch);
274 
275 
276  //printf("novalue %ld :\n",novalue);
277  //stop_execution();
278  LB=addresses(lmask,row_shift,col_shift);
279  m_ind=m_indices_with_novalues(LB,lmask->nch,novalue,IBASE,index_pixel_from_a_bin);
280  printf("STOP!\n");
281  free_longbin(LB);
282  free_longmatrix(lmask);
283 
284  //stop_execution();
285  return m_ind;
286 
287 }
288 //LONGMATRIX *indices_conditioned(long nrh, long nch,long IBASE,long (*t_index)(long r, long c,long nrh, long nch)
289 
290 
291 
292 // POINT *new_point_from_raster(long r,long c, long nrh, long nch, double lx, double ly, double nsres, double ewres, double blc_x, double blc_y, long NBASE,long (*t_index)(long r, long c,long nrh, long nch)) {
293 
294 POINT *new_point_from_raster(long r,long c, long nrh, long nch, double lx, double ly, double nsres, double ewres, double blc_x, double blc_y, long index) {
316  POINT *P;
317  double x,y;
318 
319 // index=(*t_index)(r,c,nrh,nch)+NBASE;
320  x=((double)c+lx)*ewres+blc_x;
321  y=((double)(nrh-r)+ly)*nsres+blc_y;
322 
323  P=new_point(index,x,y,NO_ELEVATION);
324 
325  return P;
326 
327 }
328 
329 LINE *new_horizontal_line_from_raster(long r,long c, long nrh, long nch, double nsres, double ewres, double blc_x, double blc_y, long iline, long ivertex1, long ivertex2){
330 
331 
356  LINE *line;
357  POINT *P1, *P2;
358 
359 
360 // index=(*t_index)(r,c,nrh,nch)+NBASE;
361  P1=new_point_from_raster(r,c,nrh,nch,LEFT,TOP,nsres,ewres,blc_x,blc_y,ivertex1);
362  P2=new_point_from_raster(r,c,nrh,nch,RIGHT,TOP,nsres,ewres,blc_x,blc_y,ivertex2);
363 
364  line=new_line_from_points(iline,P1,P2);
365 
366  free_point(P1);
367  free_point(P2);
368 
369 
370  return line;
371 
372 }
373 
374 
375 LINE *new_vertical_line_from_raster(long r,long c, long nrh, long nch, double nsres, double ewres, double blc_x, double blc_y, long iline, long ivertex1, long ivertex2){
376 
377 
402  LINE *line;
403  POINT *P1, *P2;
404 
405 
406 // index=(*t_index)(r,c,nrh,nch)+NBASE;
407  P1=new_point_from_raster(r,c,nrh,nch,LEFT,BOTTOM,nsres,ewres,blc_x,blc_y,ivertex1);
408  P2=new_point_from_raster(r,c,nrh,nch,LEFT,TOP,nsres,ewres,blc_x,blc_y,ivertex2);
409 // printf("index=%ld r=%ld c=%ld \n",iline,r,c);
410 
411  line=new_line_from_points(iline,P1,P2);
412 
413  free_point(P1);
414  free_point(P2);
415 
416  return line;
417 
418 }
419 
420 
421 LINEVECTOR *get_linevector_from_raster_grid(LONGMATRIX *i_horizontal,LONGMATRIX *i_vertical,LONGMATRIX *i_vertex, double nsres, double ewres, double blc_x, double blc_y, long novalue){
422 
443 LINEVECTOR *lines;
444 long r,c,i,count;
446 count=0;
448 for(r=i_vertical->nrl;r<=i_vertical->nrh;r++){
449  for (c=i_vertical->ncl;c<=i_vertical->nch;c++){
452  if (i_vertical->element[r][c]!=novalue) count++;
453  }
454 }
455 //stop_execution();
458 for(r=i_horizontal->nrl;r<=i_horizontal->nrh;r++){
459  for (c=i_horizontal->ncl;c<=i_horizontal->nch;c++){
461  if (i_horizontal->element[r][c]!=novalue) count++ ;
462  }
463 }
466 lines=new_linevector(count);
467 count=0;
468 
470 for(r=i_vertical->nrl;r<=i_vertical->nrh;r++){
471  for (c=i_vertical->ncl;c<=i_vertical->nch;c++){
473 // printf("indexxx=%ld r=%ld c=%ld \n",i_vertical->element[r][c],r,c);
474  if (i_vertical->element[r][c]!=novalue) lines->element[i_vertical->element[r][c]]=new_vertical_line_from_raster(r,c,i_vertical->nrh,i_vertical->nch,nsres,ewres,blc_x,blc_y,i_vertical->element[r][c],i_vertex->element[r+1][c],i_vertex->element[r][c]);
475  }
476 }
477 //stop_execution();
480 for(r=i_horizontal->nrl;r<=i_horizontal->nrh;r++){
481  for (c=i_horizontal->ncl;c<=i_horizontal->nch;c++){
483  if (i_horizontal->element[r][c]!=novalue) lines->element[i_horizontal->element[r][c]]=new_horizontal_line_from_raster(r,c,i_horizontal->nrh,i_horizontal->nch,nsres,ewres,blc_x,blc_y,i_horizontal->element[r][c],i_vertex->element[r][c],i_vertex->element[r][c+1]);
484  }
485 }
488 for (i=lines->nl;i<=lines->nh;i++){
489  if (!lines->element[i]) printf("Warning:: element %ld corresponding of lines was not assigned !!\n",i);
490 }
491 
492 return lines;
493 }
494 
495 POLYGON *new_pixel_from_raster(long index,long r, long c ,LINEVECTOR *lines, LONGMATRIX *i_horizontal,LONGMATRIX *i_vertical,double nsres, double ewres, double blc_x, double blc_y,short print) {
518  LINEVECTOR *edges;
519  POLYGON *PO;
520  LONGVECTOR *ledges;
521  POINT *centroid;
522 
523  long nrh,nch;
524 
525  nrh=i_vertical->nrh;
526  nch=i_horizontal->nch;
527 
528 
529  centroid=new_point_from_raster(r,c,nrh,nch,LX_CENTER,LY_CENTER,nsres,ewres,blc_x,blc_y,index);
532  ledges=new_longvector(4);
533 
534 
535  if ((!i_vertical->element[r][c]) || (i_vertical->element[r][c]<=0)) {
536  printf ("Warning: left vertical line missing (polygon %ld) at r=%ld c=%ld \n ",index,r,c);
537  ledges->element[1]=i_vertical->element[r][c];
538  } else {
539  ledges->element[1]=i_vertical->element[r][c];
540 
541 
542  }
543  if ((!i_vertical->element[r][c+1]) || (i_vertical->element[r][c+1]<=0)) {
544  printf ("Warning: right vertical line missing (polygon %ld) at r=%ld c=%ld \n ",index,r,c);
545  ledges->element[2]=i_vertical->element[r][c+1];
546  } else {
547  ledges->element[2]=i_vertical->element[r][c+1];
548 
549  }
550 
551  if ((!i_horizontal->element[r][c]) || (i_horizontal->element[r][c]<=0)) {
552  printf ("Warning: top horizontal line missing (polygon %ld) at r=%ld c=%ld \n",index,r,c);
553  ledges->element[3]=i_horizontal->element[r][c];
554  } else {
555  ledges->element[3]=i_horizontal->element[r][c];
556 
557  }
558 
559 
560  if ((!i_horizontal->element[r+1][c]) || (i_horizontal->element[r+1][c]<=0)) {
561  printf ("Warning: bottom horizontal line missing (polygon %ld) at r=%ld c=%ld \n",index,r,c);
562  ledges->element[4]=i_horizontal->element[r+1][c];
563  } else {
564  ledges->element[4]=i_horizontal->element[r+1][c];
565 
566  }
567 
568 
569  if (print==1) printf(" Polygon : %ld [r= %ld , c = %ld ] lines (edges): %ld, %ld, %ld, %ld \n",index,r,c,ledges->element[1],ledges->element[2],ledges->element[3],ledges->element[4]);
570 
571  edges=extract_linvector_from_linevector(ledges,lines);
572  PO=new_polygon_from_a_linevector(edges,centroid);
573 
574 
575  if (print==1) printf(" Polygon : %ld [r= %ld , c = %ld ] lines (edges): %ld, %ld, %ld, %ld was allocated \n",index,r,c,ledges->element[1],ledges->element[2],ledges->element[3],ledges->element[4]);
576  free_point(centroid);
577  free_longvector(ledges);
578  free_linevector(edges);
579 
580  return PO;
581 
582 
583 }
584 
585 POLYGONVECTOR *get_polygonvector_from_raster(LINEVECTOR *lines,LONGMATRIX *i_pixels,LONGMATRIX *i_horizontal,LONGMATRIX *i_vertical,double nsres, double ewres, double blc_x, double blc_y,long novalue, short print) {
604  POLYGONVECTOR *pv;
605  long i,r,c,count;
607  count=0;
608 
609  for(r=i_pixels->nrl;r<=i_pixels->nrh;r++){
610  for(c=i_pixels->ncl;c<=i_pixels->nch;c++){
611  if (i_pixels->element[r][c]!=novalue) count++;
612 
613  }
614  }
615 
616  pv=new_polygonvector(count);
617  count=0;
618 
619  for(r=i_pixels->nrl;r<=i_pixels->nrh;r++){
620  for(c=i_pixels->ncl;c<=i_pixels->nch;c++){
621  if (i_pixels->element[r][c]!=novalue) {
622  i=i_pixels->element[r][c];
623  pv->element[i]=new_pixel_from_raster(i,r,c,lines,i_horizontal,i_vertical,nsres, ewres,blc_x,blc_y,print);
624  }
625  }
626  }
627 
628 
629 
630 
633  for (i=pv->nl;i<=pv->nh;i++){
634  if (!pv->element[i]) printf("Warning:: element %ld corresponding of polygons was not assigned !!\n",i);
635  }
636 
637  return pv;
638 
639 }
640 
643 DOUBLEVECTOR *get_doublevector_from_doublematrix(LONGMATRIX *indices,DOUBLEMATRIX *M, double novalue){
654  DOUBLEVECTOR *v;
655  long r,c,i,cnt;
656 
657  if ((indices->nrh!=M->nrh) || (indices->nch!=M->nch)) printf("Error:: in get_doublevector_from_doublematrix indices [%ld,%ld] and M [%ld,%ld] has different sizes! \n",indices->nrh,indices->nch,M->nrh,M->nch);
658  cnt=0;
659  for (r=M->nrl;r<=M->nrh;r++){
660  for (c=M->ncl;c<=M->nch;c++){
661  if (M->element[r][c]!=novalue) cnt++;
662 
663  }
664  }
665 
666  v=new_doublevector(cnt);
667 
668  for (i=v->nl;i<=v->nh;i++){
669  v->element[i]=INIT_VALUE;
670  }
671  if (cnt>M->nrh*M->nch) printf ("Error:: Error:: in get_doublevector_from_doublematrix vector elements [%ld] execeed number doublematrix elements [%ld]!!\n",v->nh,M->nrh*M->nch);
672  for (r=M->nrl;r<=M->nrh;r++){
673  for (c=M->ncl;c<=M->nch;c++){
674  i=indices->element[r][c];
675  if ((i<v->nl) || (i>v->nh)) {
676  if (M->element[r][c]!=novalue) printf ("Error:: in get_doublevector_from_doublevector index %ld exceeds size of matrix [%ld,%ld] at %ld,%ld \n",i,M->nrh,M->nch,r,c);
677  }else{
678  v->element[i]=M->element[r][c];
679  }
680  }
681  }
682 
683  for (i=v->nl;i<=v->nh;i++){
684  if (v->element[i]==INIT_VALUE) printf("Error:: in get_doublevector_from_doublematrix index %ld was not assigned (%lf) !!\n",i,v->element[i]);
685  }
686 
687  return v;
688 }
689 
690 
691 
692 DOUBLEMATRIX *get_doublematrix_from_doublevector(DOUBLEVECTOR *v, LONGMATRIX *indices, DOUBLEMATRIX *Mref, double novalue) {
708 DOUBLEMATRIX *M;
709 long i,r,c;
710 
711 M=new_doublematrix(Mref->nrh,Mref->nch);
712 
713 if ((indices->nrh!=M->nrh) || (indices->nch!=M->nch)) printf("Error:: in get_doublematrix_from_doublevector indices [%ld,%ld] and M [%ld,%ld] has different sizes! \n",indices->nrh,indices->nch,M->nrh,M->nch);
714 
715 for (r=M->nrl;r<=M->nrh;r++){
716  for (c=M->ncl;c<=M->nch;c++) {
717  M->element[r][c]=INIT_VALUE;
718  }
719 }
720 
721 for (r=M->nrl;r<=M->nrh;r++){
722  for (c=M->ncl;c<=M->nch;c++){
723  i=indices->element[r][c];
724 
725  if ((i>v->nh) || (i<v->nl)) {
726  if (Mref->element[r][c]!=novalue) printf ("Error:: in get_doublematrix_from_doublevector index %ld [%ld %ld] exceeds size of matrix [%ld,%ld] at %ld,%ld ",i,v->nl,v->nh,M->nrh,M->nch,r,c);
727  }else{
728  M->element[r][c]=v->element[i];
729  }
730  if ((M->element[r][c]==INIT_VALUE) && (Mref->element[r][c]==novalue)) M->element[r][c]=novalue;
731  }
732 }
733 
734 for (r=M->nrl;r<=M->nrh;r++){
735  for (c=M->ncl;c<=M->nch;c++){
736  if (M->element[r][c]==INIT_VALUE) printf("Error:: in get_doublematrix_from_doublevector matrix element %ld, %ld [%ld,%ld] was not assigned \n",r,c,M->nrh,M->nch);
737  }
738  }
739 
740 
741 return M;
742 
743 }
744 
745 DOUBLEVECTOR *read_doublevector_from_raster(short a, char *filename, DOUBLEMATRIX *Mref, T_INIT *UVref,LONGMATRIX *indices){
747 
762  DOUBLEMATRIX *M;
763  DOUBLEVECTOR *v;
764 
765  M=read_map(a,filename,Mref,UVref);
766  if (M==NULL) return NULL;
767 
768  v=get_doublevector_from_doublematrix(indices,M,UVref->V->element[2]);
769  free_doublematrix(M);
770 
771  return v;
772 
773 }
774 
775 int write_raster_from_doublevector(char *filename, DOUBLEVECTOR *v, T_INIT *UVref, LONGMATRIX *indices, DOUBLEMATRIX *Mref){
793  DOUBLEMATRIX *M;
794 
795  M=get_doublematrix_from_doublevector(v,indices,Mref,UVref->V->element[2]);
796  write_map(filename, FLOATING_POINT_TYPE,MAP_FORMAT,M,UVref);
797 
798  free_doublematrix(M);
799 
800  return 0;
801 
802 }
803 
804 DOUBLEMATRIX *get_doublematrix_from_mapseries(LONGMATRIX *indices,DOUBLETENSOR *mapseries, T_INIT *UVref){
811  long l,c;
812  DOUBLEMATRIX *map, *mv;
813  DOUBLEVECTOR *v;
814 
815  for (l=mapseries->ndl;l<=mapseries->ndh;l++){
816  map=extract_a_new_map(mapseries,l,UVref);
817  v=get_doublevector_from_doublematrix(indices,map,UVref->V->element[2]);
821  if (l==mapseries->ndl) mv=new_doublematrix(mapseries->ndh,v->nh);
822  for (c=mv->ncl;c<=mv->nch;c++){
823  mv->element[l][c]=v->element[c];
824  }
825  free_doublematrix(map);
826  free_doublevector(v);
827  }
828 
829  return mv;
830 
831 }
832 
833