多分辨率风场数据可视化

/************中心环形矢量场*马鞍矢量场*****卷积白噪声纹理***********/#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <cv.h>
#include <highgui.h>
#include <iostream>
#include "netcdfcpp.h"
using namespace std;#define  SQUARE_FLOW_FIELD_SZ   400
#define  DISCRETE_FILTER_SIZE   2048        //离散的滤波尺寸
#define  LOWPASS_FILTR_LENGTH   10.00000f   //低通滤波长度
#define  LINE_SQUARE_CLIP_MAX   100000.0f   //线性平方夹
#define  VECTOR_COMPONENT_MIN   0.050000f //矢量分量最小值void     SyntheszSaddle(int  n_xres,  int     n_yres,  float*   pVectr);
void     NormalizVectrs(int  n_xres,  int     n_yres,  float*   pVectr,float* vecSize);
void     GenBoxFiltrLUT(int  LUTsiz,  float*  p_LUT0,  float*   p_LUT1);
void     MakeWhiteNoise(int  n_xres,  int     n_yres,  unsigned char*  pNoise);
void     FlowImagingLIC(int  n_xres,  int     n_yres,  float*   pVectr,   unsigned char*  pNoise,  unsigned char*  pImage,  float*  p_LUT0,  float*  p_LUT1,  float  krnlen);
void     WriteImage2PPM(int  n_xres,  int     n_yres,  unsigned char*  pImage,char*  f_name);
void  color(int n_xres, int n_yres,unsigned char *pImage,float* vecSize);double maxvecmag;void  main()
{
//  int             n_xres = 721;
//  int             n_yres = 281;int               n_xres = 1441;int              n_yres = 561;//    int             n_xres = 2881;
//  int             n_yres = 1121;
//
//  int             n_xres = 5761;
//  int             n_yres = 2241;//   int             n_xres = 11521;
//      int             n_yres = 2241;float*           pVectr = (float*         ) malloc( sizeof(float        ) * n_xres * n_yres * 2 );float*            p_LUT0 = (float*        ) malloc( sizeof(float        ) * DISCRETE_FILTER_SIZE);float*         p_LUT1 = (float*        ) malloc( sizeof(float        ) * DISCRETE_FILTER_SIZE);unsigned char* pNoise = (unsigned char* ) malloc( sizeof(unsigned char) * n_xres * n_yres     );unsigned char*    pImage = (unsigned char* ) malloc( sizeof(unsigned char) * n_xres * n_yres     );float*    vecSize = (float* ) malloc( sizeof(float) * n_xres*n_yres );SyntheszSaddle(  n_xres,   n_yres,   pVectr);//CenterFiled(n_xres, n_yres, pVectr);//包含矢量归一化NormalizVectrs(n_xres, n_yres, pVectr,vecSize);//所以这就不用矢量归一化了，因为之前的马鞍矢量场生成函数里没有归一化，才有此步的MakeWhiteNoise(n_xres, n_yres, pNoise);GenBoxFiltrLUT(DISCRETE_FILTER_SIZE, p_LUT0, p_LUT1);FlowImagingLIC(n_xres, n_yres, pVectr, pNoise, pImage, p_LUT0, p_LUT1, LOWPASS_FILTR_LENGTH);
//  color(n_xres, n_yres,pImage,vecSize);//WriteImage2PPM(n_xres, n_yres, pImage, "LIC.ppm");
//  WriteImage2PPM(n_xres, n_yres, pImage, "LIC_721_281.ppm");WriteImage2PPM(n_xres, n_yres, pImage, "LIC1441_561.ppm");
//  WriteImage2PPM(n_xres, n_yres, pImage, "LIC2281_1121.ppm");//WriteImage2PPM(n_xres, n_yres, pImage, "LIC5761_1121.ppm");//WriteImage2PPM(n_xres, n_yres, pImage, "LIC_11521_2241.ppm");//system("pause");free(pVectr);  pVectr = NULL;free(p_LUT0);    p_LUT0 = NULL;free(p_LUT1);    p_LUT1 = NULL;free(pNoise);    pNoise = NULL;free(pImage);    pImage = NULL;
}///    中心环形矢量场图形   synthesize a saddle-shaped vector field     ///
void    SyntheszSaddle(int  n_xres,  int  n_yres,  float*  pVectr)
{       static const int LatNum = n_yres;static const int LonNum = n_xres;static const int Time = 1;static const int TMP = Time*LonNum*LatNum;//NcFile dataReadFile("global_721_281.nc",NcFile::ReadOnly);NcFile dataReadFile("global_1441_561.nc",NcFile::ReadOnly);
//  NcFile dataReadFile("global_2881_1121.nc",NcFile::ReadOnly);//    NcFile dataReadFile("global_5761_2241.nc",NcFile::ReadOnly);//NcFile dataReadFile("global_11521_2241.nc",NcFile::ReadOnly);if (!dataReadFile.is_valid()){std::cout<<"couldn't open file!"<<std::endl;}double *Tmp_UX = new double[TMP];double *Tmp_VY = new double[TMP];double *Tmp_LAT = new double[TMP];double *Tmp_LON = new double[TMP];NcVar *dataTmp_LAT = dataReadFile.get_var("LAT"); NcVar *dataTmp_LON = dataReadFile.get_var("LON");    NcVar *dataTmp_UX = dataReadFile.get_var("UX");  NcVar *dataTmp_VY = dataReadFile.get_var("VY");  dataTmp_LAT->get(Tmp_LAT,LatNum,LatNum);dataTmp_LON->get(Tmp_LON,LonNum,LonNum);dataTmp_UX->get(Tmp_UX,Time,LatNum,LonNum);dataTmp_VY->get(Tmp_VY,Time,LatNum,LonNum);for(int  j = 0;  j < n_yres;  j ++)for(int  i = 0;  i < n_xres;  i ++){   int  index = (  (n_yres - 1 - j) * n_xres + i  )  <<  1;//int index = j*n_yres+i;pVectr[index    ] = Tmp_UX[j*LonNum+i];pVectr[index +1   ]= Tmp_VY[j*LonNum+i];}    delete []Tmp_UX;delete []Tmp_VY;delete []Tmp_LAT;}///       normalize the vector field     ///
// void    NormalizVectrs(int  n_xres,  int  n_yres,  float*  pVectr)
// {
//
//
//  for(int  j = 0;  j < n_yres;  j ++)
//      for(int  i = 0;  i < n_xres;  i ++)
//      {
//          int     index = (j * n_xres + i) << 1;
//          float   vcMag = float(  sqrt( double(pVectr[index] * pVectr[index] + pVectr[index + 1] * pVectr[index + 1]) )  );
//
//          float   scale = (vcMag == 0.0f) ? 0.0f : 1.0f / vcMag;//矢量大小归一化后的矢量值
//          //pVectr[index    ] *= scale;
//          pVectr[index    ]=pVectr[index    ]*scale;
//          //  cout<<"pVectr[index    ]="<<pVectr[index    ];
//          pVectr[index + 1] *= scale;
//          //cout<<"pVectr[index    ]="<<pVectr[index  +1  ];
//
//      }
// }
void    NormalizVectrs(int n_xres,  int  n_yres,  float*  pVectr,float* vecSize)
{   for(int  j = 0;  j < n_yres;  j ++)for(int     i = 0;  i < n_xres;  i ++){  int     index = (j * n_xres + i) << 1;float vcMag = float(  sqrt( double(pVectr[index] * pVectr[index] + pVectr[index + 1] * pVectr[index + 1]) )  );vecSize[j * n_xres + i]=vcMag;if (vcMag<100&&vcMag>maxvecmag){maxvecmag=vcMag;}float  scale = (vcMag == 0.0f) ? 0.0f : 1.0f / vcMag;pVectr[index    ] *= scale*5.5;//????????????????????????????????????????????????????????????原来问题出在这pVectr[index + 1] *= scale*5.5;}
}///        make white noise as the LIC input texture     ///
// void MakeWhiteNoise(int  n_xres,  int  n_yres,  unsigned char*  pNoise)
// {
//  IplImage * NoiseImg=cvCreateImage(cvSize(n_xres,n_yres),IPL_DEPTH_8U,1);
//  CvScalar s;
//
//  for(int  j = 0;   j < n_yres;  j ++)
//  {
//      for(int  i = 0;   i < n_xres;  i ++)
//
// //   for(int  j = 0;   j < n_yres;  j=j +10)//产生稀疏白噪声
// //   {
// //       for(int  i = 0;   i < n_xres;  i=i+10)
//
//      {
//          int  r = rand();
//          r = (  (r & 0xff) + ( (r & 0xff00) >> 8 )  ) & 0xff;
//          pNoise[j * n_xres + i] = (unsigned char) r;
//          s = cvGet2D(NoiseImg,i,j);
//          s.val[0]=r;
//          s.val[1]=r;
//          s.val[2]=r;
//          cvSet2D(NoiseImg,i,j,s);
//      }
//  }
//
// }void    MakeWhiteNoise(int  n_xres,  int  n_yres,  unsigned char*  pNoise)
{       for(int  j = 0;   j < n_yres;  j ++)for(int  i = 0;   i < n_xres;  i ++){   int  r = rand();r = (  (r & 0xff) + ( (r & 0xff00) >> 8 )  ) & 0xff;pNoise[j * n_xres + i] = (unsigned char) r;}
}
///     generate box filter LUTs     ///
void    GenBoxFiltrLUT(int  LUTsiz,  float*  p_LUT0,  float*  p_LUT1)
{       for(int  i = 0;  i < LUTsiz;  i ++)  p_LUT0[i] = p_LUT1[i] = i;
}void color(int n_xres,int n_yres, unsigned char* pImage,float* vecSize)
{IplImage * licImage = cvCreateImage(cvSize(n_xres,n_yres),IPL_DEPTH_8U,3);IplImage* img = cvLoadImage("11.jpg",1);int k = 0;double magind;double mag;double newMag;double x = 0.1;//x为非线性映射因子，且x！=1CvScalar colorTable[500];CvScalar s,s1;for (int i = 0;i < img->width;i++){s = cvGet2D(img,1,i);colorTable[i] =s;}for (int j=0;j<n_yres;j++){for (int i= 0;i<n_xres;i++){if (k>=img->width){k=0;}double  scale=  pImage[j * n_xres + i];mag = vecSize[j * n_xres + i];//********矢量大小归一化******if (mag<100){magind = (mag/maxvecmag);}//非线性颜色增强LICnewMag =(pow(x,magind)-1)/(x-1);s = cvGet2D(licImage,j,i);//渐变颜色映射表int k = int(newMag*446); s1.val[0]=colorTable[k].val[0]*(k+1-newMag*446)+colorTable[k+1].val[0]*(newMag*446-k);s1.val[1]=colorTable[k].val[1]*(k+1-newMag*446)+colorTable[k+1].val[1]*(newMag*446-k);s1.val[2]=colorTable[k].val[2]*(k+1-newMag*446)+colorTable[k+1].val[2]*(newMag*446-k);s1.val[0]*=scale;s1.val[1]*=scale;s1.val[2]*=scale;/*cout<<"s1.val[3]="<<s1.val[3]<<endl;*/cvSet2D(licImage,j,i,s1);}}//Mat AlphaImage= imread("s.jpg");//cv::Mat AlphaImage = imread("licImage",1);cvNamedWindow("lic_three channles",0);cvShowImage("lic_three channles",licImage);cvWaitKey(0);system("pause");cvDestroyWindow("lic_three channles");cvReleaseImage(&licImage);
}
///     write the LIC image to a PPM file     ///
void    WriteImage2PPM(int  n_xres,  int  n_yres,  unsigned char*  pImage,  char*  f_name)
{   FILE*   o_file;if(   ( o_file = fopen(f_name, "w") )  ==  NULL   )  {  printf("Can't open output file\n");  return;  }fprintf(o_file, "P6\n%d %d\n255\n", n_xres, n_yres);for(int  j = 0;  j < n_yres;  j ++)for(int  i = 0;  i < n_xres;  i ++){unsigned  char   unchar = pImage[j * n_xres + i];//某点像素的灰度纹理值fprintf(o_file, "%c%c%c", unchar, unchar, unchar);//}fclose (o_file);   o_file = NULL;
}///        flow imaging (visualization) through Line Integral Convolution     ///
void    FlowImagingLIC(int     n_xres,  int     n_yres,  float*  pVectr,  unsigned char*  pNoise,  unsigned char*  pImage,  float*  p_LUT0,  float*  p_LUT1,  float   krnlen)
{   int     vec_id;                     ///ID in the VECtor buffer (for the input flow field)int        advDir;                     ///ADVection DIRection (0: positive;  1: negative)int       advcts;                     ///number of ADVeCTion stepS per direction (a step counter)int      ADVCTS = int(krnlen * 3);  ///MAXIMUM number of advection steps per direction to break dead loops  //跳出死循环的条件float vctr_x;                     ///x-component  of the VeCToR at the forefront pointfloat   vctr_y;                     ///y-component  of the VeCToR at the forefront pointfloat   clp0_x;                     ///x-coordinate of CLiP point 0 (current)float  clp0_y;                     ///y-coordinate of CLiP point 0 (current)float  clp1_x;                     ///x-coordinate of CLiP point 1 (next   )float  clp1_y;                     ///y-coordinate of CLiP point 1 (next   )float  samp_x;                     ///x-coordinate of the SAMPle in the current pixelfloat samp_y;                     ///y-coordinate of the SAMPle in the current pixelfloat tmpLen;                     ///TeMPorary LENgth of a trial clipped-segmentfloat segLen;                     ///SEGment   LENgthfloat    curLen;                     ///CURrent   LENgth of the streamlinefloat  prvLen;                     ///PReVious  LENgth of the streamline       float   W_ACUM;                     ///ACcuMulated Weight from the seed to the current streamline forefrontfloat    texVal;                     ///TEXture VALuefloat   smpWgt;                     ///WeiGhT of the current SaMPlefloat    t_acum[2];                  ///two ACcUMulated composite Textures for the two directions, perspectively 两个方向的卷积和float   w_acum[2];                  ///two ACcUMulated Weighting values   for the two directions, perspectively 两个方向的权重和float*  wgtLUT = NULL;             ///WeiGhT Look Up Table pointing to the target filter LUT权重查找表float len2ID = (DISCRETE_FILTER_SIZE - 1) / krnlen;  ///map a curve LENgth TO an ID in the LUT///for each pixel in the 2D output LIC image///for(int  j = 0;     j < n_yres;  j ++)for(int  i = 0;     i < n_xres;  i ++){   ///init the composite texture accumulators and the weight accumulators///每一个像素点为起始点，初始化一次权重和卷积和t_acum[0] = t_acum[1] = w_acum[0] = w_acum[1] = 0.0f;//初始化正反方向卷积和及权重和///for either advection direction///分别计算正反方向的卷积和及权重和for(advDir = 0;  advDir < 2;  advDir ++){  ///init the step counter, curve-length measurer, and streamline seed/////初始化当前方向上前进的步数和当前流线的总长advcts = 0;//前进的步数curLen = 0.0f;clp0_x = i + 0.5f;clp0_y = j + 0.5f;///access the target filter LUT///LUT显示查找表wgtLUT = (advDir == 0) ? p_LUT0 : p_LUT1;///until the streamline is advected long enough or a tightly  spiralling center / focus is encountered///while( curLen < krnlen && advcts < ADVCTS ) //??????{///access the vector at the sample///vec_id = ( int(clp0_y) * n_xres + int(clp0_x) ) << 1;vctr_x = pVectr[vec_id    ];vctr_y = pVectr[vec_id + 1];///in case of a critical point///遇到零向量，结束循环if( vctr_x == 0.0f && vctr_y == 0.0f ){    t_acum[advDir] = (advcts == 0) ? 0.0f : t_acum[advDir];         ///this line is indeed unnecessaryw_acum[advDir] = (advcts == 0) ? 1.0f : w_acum[advDir];break;}///negate the vector for the backward-advection case///相反的方向取相反的方向vctr_x = (advDir == 0) ? vctr_x : -vctr_x;vctr_y = (advDir == 0) ? vctr_y : -vctr_y;///clip the segment against the pixel boundaries --- find the shorter from the two clipped segments//////replace  all  if-statements  whenever  possible  as  they  might  affect the computational speed///segLen = LINE_SQUARE_CLIP_MAX;//cout<<"segLen="<<segLen<<endl;//cout<<"VECTOR_COMPONENT_MIN="<<LINE_SQUARE_CLIP_MAX<<endl;segLen = (vctr_x < -VECTOR_COMPONENT_MIN) ? ( int(     clp0_x         ) - clp0_x ) / vctr_x : segLen;//int（0.5）=0segLen = (vctr_x >  VECTOR_COMPONENT_MIN) ? ( int( int(clp0_x) + 1.5f ) - clp0_x ) / vctr_x : segLen;segLen = (vctr_y < -VECTOR_COMPONENT_MIN) ?    (      (    (  tmpLen = ( int(     clp0_y)          - clp0_y ) / vctr_y  )  <  segLen    ) ? tmpLen : segLen      ) : segLen;segLen = (vctr_y >  VECTOR_COMPONENT_MIN) ?(      (    (  tmpLen = ( int( int(clp0_y) + 1.5f ) - clp0_y ) / vctr_y  )  <  segLen    ) ? tmpLen : segLen      ) : segLen;///update the curve-length measurers///prvLen = curLen;curLen+= segLen;segLen+= 0.0004f;///check if the filter has reached either end///segLen = (curLen > krnlen) ? ( (curLen = krnlen) - prvLen ) : segLen;///obtain the next clip point///clp1_x = clp0_x + vctr_x * segLen;clp1_y = clp0_y + vctr_y * segLen;///obtain the middle point of the segment as the texture-contributing sample///samp_x = (clp0_x + clp1_x) * 0.5f;samp_y = (clp0_y + clp1_y) * 0.5f;///obtain the texture value of the sample///texVal = pNoise[ int(samp_y) * n_xres + int(samp_x) ];///update the accumulated weight and the accumulated composite texture (texture x weight)///W_ACUM = wgtLUT[ int(curLen * len2ID) ];smpWgt = W_ACUM - w_acum[advDir];         w_acum[advDir]  = W_ACUM;                              t_acum[advDir] += texVal * smpWgt;///update the step counter and the "current" clip point///advcts ++;clp0_x = clp1_x;clp0_y = clp1_y;///check if the streamline has gone beyond the flow field///if( clp0_x < 0.0f || clp0_x >= n_xres || clp0_y < 0.0f || clp0_y >= n_yres)  break;}    }///normalize the accumulated composite texture///texVal = (t_acum[0] + t_acum[1]) / (w_acum[0] + w_acum[1]);///clamp the texture value against the displayable intensity range [0, 255]texVal = (texVal <   0.0f) ?   0.0f : texVal;texVal = (texVal > 255.0f) ? 255.0f : texVal; pImage[j * n_xres + i] = (unsigned char) texVal;}
}

多分辨率风场数据可视化相关推荐

python--转换wrf输出的风场数据为网页可视化的json格式
前言: 一般网页可视化风场中的数据都是json格式,而如果我们希望将wrf模式模拟输出的风场数据在网页中进行展示,这就需要先将wrfoutput数据转换为网页可以识别的json格式. 这里主要需要用到 ...
风场可视化：风场数据
引子了解 WebGL 基础之后,接着去看获取解析风场数据的逻辑,又遇到问题. 系统:MacOS 版本:11.6 Origin My GitHub 安装 ecCodes 在文章示例源库的说明中,首先要 ...
百度大数据可视化产品矩阵
百度数据可视化实验室的产品矩阵如下图所示,内容涵盖基础库.各种可视化产品以及应用产品. 在官网中,百度数据可视化实验室分享了其发展和规划: 基于基础的可视化规范,依托 ZRender.ClayGL 基 ...
DataV数据可视化功能特性
使用DataV制作实时销售数据可视化大屏 (本课程可以帮助数据分析师学习数据可视化大屏的制作,包括制作的方法.设计原则等基础知识,并提供一个微项目,使用数加的DataV基于ABC公司的经营数据,快速构 ...
数据可视化、信息可视化与知识可视化
数据可视化.信息可视化与知识可视化 (2011-07-23 12:28:17) 标签: 校园分类: 工作篇数据可视化简介数据可视化是关于数据之视觉表现形式的研究:其中,这种数据的视觉表现形 ...
31个有点意思数据可视化作品！
在一个信息大爆炸的时代,每天都有很多的新消息.新发现.新趋势向我们狂轰乱炸而来.在这个过程中,我们既是数据的生产者,也是数据的使用者,然而初次获取和存储的原始数据总是杂乱无章的. 要想数据达到生动有趣 ...
数据可视化作品有哪些
在一个信息大爆炸的时代,每天都有很多的新消息.新发现.新趋势向我们狂轰乱炸而来.在这个过程中,我们既是数据的生产者,也是数据的使用者,然而初次获取和存储的原始数据总是杂乱无章的. 要想数据达到生动有趣 ...
【案例教程】Python气象海洋数据可视化到常见数据分析方法（折线图、柱状图、errorbar图、流场矢量、散点图、风玫瑰图、流场矢量、填色及等值线+地图）
[查看原文]Python在气象与海洋中的实践技术应用 Python是功能强大.免费.开源,实现面向对象的编程语言,能够在不同操作系统和平台使用,简洁的语法和解释性语言使其成为理想的脚本语言.除了标准库 ...
大数据可视化之气象数据可视化（雷达、云图、落区、等值面）
气象数据包含多类数据,包括地面气象资料地面逐小时资料高空气象资料高空定时资料卫星探测资料大海区云图(IR1) | 扫描辐射计L1数据天气雷达探测资料基本反射率数值预报模式产品T639模式产品 ...
python气象数据可视化学习记录1——基于ERA5数据画风场和海平面气压填色叠加图
python气象数据可视化学习记录1--基于ERA5数据画风场和海平面气压填色叠加图 1. 写在前面 2. 图片效果 3. 逐步代码解析 3.1导入库 3.2 读取NC格式数据 3.3 对数据进行加工 ...

多分辨率风场数据可视化

多分辨率风场数据可视化相关推荐

最新文章

热门文章