天天看點

[GDAL]讀取HDF格式的calipso資料

探測地球雲層分布的CloudSat和CALIPSO衛星

http://www.nasa.gov/mission_pages/calipso/main/index.html

http://www.nasa.gov/topics/earth/features/calipso-1billion.html

Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO)

The occasion was on CALIPSO's 9,491st orbit of Earth. Since its launch, on April 23, 2006, the satellite has traveled 260,198,685.3 miles. It has generated data that would fill almost 3,500 DVDs or more than 24,000 CDs, delighting scientists and researchers from 35 countries who are using the information to prove any number of theories about the atmosphere, about weather, about the past and how it will relate to the future.

資料格式

[GDAL]讀取HDF格式的calipso資料
[GDAL]讀取HDF格式的calipso資料
[GDAL]讀取HDF格式的calipso資料
讀取代碼:
[GDAL]讀取HDF格式的calipso資料
[GDAL]讀取HDF格式的calipso資料

1             double lonFrom = 0;
  2             double latFrom = 0;
  3             double lonTo = 0;
  4             double latTo = 0;
  5             OSGeo.GDAL.Gdal.AllRegister();
  6             Gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "Yes");
  7             string path = string.Format(@"C:\Users\Desktop\lidar_track_image\CAL_LID_L1-ValStage1-V3-02.2013-01-01T18-30-59ZN_Subset.hdf", WorldSettings.StartupDirectory);
  8             //path = string.Format(@"{0}\data\CAL_LID_L1-ValStage1-V3-02.2013-01-01T18-30-59ZN_Subset.hdf", WorldSettings.StartupDirectory);
  9             byte[] pathUtf16 = Encoding.Unicode.GetBytes(path);
 10 
 11             byte[] byteUtf8 = Encoding.Convert(Encoding.Unicode, Encoding.UTF8, pathUtf16);
 12 
 13             string pathUtf8 = Encoding.UTF8.GetString(byteUtf8);
 14 
 15             Encoding encoud = Encoding.GetEncoding(936);
 16             byte[] byteGb2312 = Encoding.Convert(Encoding.Unicode, encoud, pathUtf16);
 17             string pathGb2312 = encoud.GetString(byteGb2312);
 18 
 19             if (!File.Exists(path))
 20             {
 21                 MessageBox.Show("Calpso資料不存在", "提示", MessageBoxButtons.OK, MessageBoxIcon.Information);
 22                 return;
 23             }
 24             //讀取檔案,擷取子資料集資訊
 25         Dataset ds = Gdal.Open(path, Access.GA_ReadOnly);
 26             string[] stMetadata = ds.GetMetadata("");
 27             string[] strMetadata = ds.GetMetadata("SUBDATASETS");
 28             string info = string.Empty;
 29             for (int i = 0; i < strMetadata.Length; i++)
 30             {
 31                 info += strMetadata[i] + "\n";
 32             }
 33             info += "讀取波段資料:----------------\n";
 34             //緯度資料
 35             string tmpstr = strMetadata[0];
 36             tmpstr = tmpstr.Substring(tmpstr.IndexOf("=") + 1);
 37             Dataset dsLat = Gdal.Open(tmpstr, Access.GA_ReadOnly);
 38             string strMetadata2 = dsLat.GetDescription();
 39 
 40             info += string.Format("像素數目:{0} X {1}", dsLat.RasterXSize, dsLat.RasterYSize);
 41             int count = dsLat.RasterCount;
 42             Band bandLat = dsLat.GetRasterBand(1);
 43             int Cols = dsLat.RasterXSize;
 44             int Rows = dsLat.RasterYSize;
 45             int width = dsLat.RasterXSize;
 46             int height = dsLat.RasterYSize;
 47             double min, max, no;
 48             int hasvalue;
 49             info += string.Format("類型Lat:{0}\n", bandLat.DataType);
 50             if (bandLat.DataType == DataType.GDT_Float32)
 51             {
 52                 bandLat.GetMaximum(out max, out hasvalue);
 53                 bandLat.GetMinimum(out min, out hasvalue);
 54                 bandLat.GetNoDataValue(out no, out hasvalue);
 55                 float[] data = new float[Cols * Rows];
 56                 bandLat.ReadRaster(0, 0, width, height, data, Cols, Rows, 0, 0);
 57 
 58                 string ddd = string.Empty;
 59                 double min1, max1, mean, hasvalue1;
 60                 OSGeo.GDAL.Gdal.GDALProgressFuncDelegate dele = new Gdal.GDALProgressFuncDelegate(GDALProgressFunc);
 61                 bandLat.ComputeStatistics(false, out min1, out max1, out hasvalue1, out mean, dele, ddd);
 62                 lonFrom = data[0];
 63                 lonTo = data[data.Length - 1];
 64             }
 65             //經度資料
 66             string tmpstrLon = strMetadata[2];
 67             tmpstrLon = tmpstrLon.Substring(tmpstrLon.IndexOf("=") + 1);
 68             Dataset dsLon = Gdal.Open(tmpstrLon, Access.GA_ReadOnly);
 69             info += string.Format("像素數目:{0} X {1}", dsLon.RasterXSize, dsLon.RasterYSize);
 70             Band bandLon = dsLon.GetRasterBand(1);
 71             int Cols2 = dsLon.RasterXSize;
 72             int Rows2 = dsLon.RasterYSize;
 73             int width2 = dsLon.RasterXSize;
 74             int height2 = dsLon.RasterYSize;
 75             double min2, max2, no2;
 76             int hasvalue2;
 77             info += string.Format("類型Lon:{0}\n", bandLon.DataType);
 78             if (bandLon.DataType == DataType.GDT_Float32)
 79             {
 80                 bandLon.GetMaximum(out max2, out hasvalue2);
 81                 bandLon.GetMinimum(out min2, out hasvalue2);
 82                 bandLon.GetNoDataValue(out no2, out hasvalue2);
 83                 float[] data2 = new float[Cols2 * Rows2];
 84                 bandLon.ReadRaster(0, 0, width2, height2, data2, Cols2, Rows2, 0, 0);
 85 
 86                 string ddd = string.Empty;
 87                 double min1, max1, mean, hasvalue1;
 88                 OSGeo.GDAL.Gdal.GDALProgressFuncDelegate dele = new Gdal.GDALProgressFuncDelegate(GDALProgressFunc);
 89                 bandLon.ComputeStatistics(false, out min1, out max1, out hasvalue1, out mean, dele, ddd);
 90                 latFrom = data2[0];
 91                 latTo = data2[data2.Length - 1];
 92             }
 93             //讀取雷達衰減系數
 94             string tmpStrBackScatter = strMetadata[70];
 95             tmpStrBackScatter = tmpStrBackScatter.Substring(tmpStrBackScatter.IndexOf("=") + 1);
 96             Dataset dsBackScatter = Gdal.Open(tmpStrBackScatter, Access.GA_ReadOnly);
 97             info += string.Format("像素數目:{0} X {1}", dsBackScatter.RasterXSize, dsBackScatter.RasterYSize);
 98             Band bandBackScatter = dsBackScatter.GetRasterBand(1);
 99             int ColsBackScatter = dsBackScatter.RasterXSize;
100             int RowsBackScatter = 3000;// ds270.RasterYSize;讀取3000個像素
101             int widthBackScatter = dsBackScatter.RasterXSize;
102             int heightBackScatter = 3000;// ds270.RasterYSize;
103             double min270, max270, no270;
104             int hasvalue270;
105             info += string.Format("類型衰減系數:{0}\n", bandBackScatter.DataType);
106             if (bandBackScatter.DataType == DataType.GDT_Float32)
107             {
108                 bandBackScatter.GetMaximum(out max270, out hasvalue270);
109                 bandBackScatter.GetMinimum(out min270, out hasvalue270);
110                 bandBackScatter.GetNoDataValue(out no270, out hasvalue270);
111                 string ddd = string.Empty;
112                 double mean, has;
113 
114                 OSGeo.GDAL.Gdal.GDALProgressFuncDelegate dele = new Gdal.GDALProgressFuncDelegate(GDALProgressFunc);
115                 bandBackScatter.ComputeStatistics(false, out min270, out max270, out has, out mean, dele, ddd);
116                 float[] data270 = new float[ColsBackScatter * RowsBackScatter];
117                 bandBackScatter.ReadRaster(0, 0, widthBackScatter, heightBackScatter, data270, ColsBackScatter, RowsBackScatter, 0, 0);
118 
119                 int x1width = heightBackScatter;
120                 int y1height = widthBackScatter;//583
121                 Bitmap bitmap = new Bitmap(x1width, y1height, System.Drawing.Imaging.PixelFormat.Format32bppArgb);
122                 int iPixelSize = 4;
123 
124                 CreateColorRamp createColorRamp = CreateColorRamps();
125                 BitmapData bitmapdata = bitmap.LockBits(new Rectangle(0, 0, x1width, y1height), ImageLockMode.ReadWrite, bitmap.PixelFormat);
126                 try
127                 {
128                     unsafe
129                     {
130                         for (int y = 0; y < y1height; y++)
131                         {
132                             byte* row = (byte*)bitmapdata.Scan0 + (y * bitmapdata.Stride);
133 
134                             for (int x = 0; x < x1width; x++)
135                             {
136                                 //int posSrc = x + y * x1width;
137                                 int posSrc = y + x * y1height;//圖像順時針旋轉90度
138                                 double item = data270[posSrc];
139                                 if (item == no270)
140                                 {
141                                     Color? getColor = getColor = createColorRamp.GetColor(0); ;
142                                     row[x * iPixelSize] = (byte)getColor.Value.B;
143                                     row[x * iPixelSize + 1] = (byte)getColor.Value.G;
144                                     row[x * iPixelSize + 2] = (byte)getColor.Value.R;
145                                     row[x * iPixelSize + 3] = getColor.Value.A;
146                                 }
147                                 else
148                                 {
149                                     Color? getColor = null;
150                                     if (item < 0.0001 || item > 0.1)
151                                     {
152                                         getColor = createColorRamp.GetColor(0);
153                                     }
154                                     else
155                                     {
156                                         getColor = createColorRamp.GetColor(item);
157                                     }
158                                     if (getColor != null)
159                                     {
160                                         row[x * iPixelSize] = getColor.Value.B;
161                                         row[x * iPixelSize + 1] = getColor.Value.G;
162                                         row[x * iPixelSize + 2] = getColor.Value.R;
163                                         row[x * iPixelSize + 3] = getColor.Value.A;
164                                     }
165                                 }
166                             }
167                         }
168                     }
169                 }
170                 catch
171                 {
172                 }
173                 finally
174                 {
175                     bitmap.UnlockBits(bitmapdata);
176                 }
177                 bitmap.Save("D:\\b.png");
178                 Vector3d v1 = new Vector3d(lonFrom, latFrom, 0);
179                 Vector3d v2 = new Vector3d(lonTo, latTo, 0);
180                 Vector3d dir = v2 - v1;
181 
182                 Vector3d v3 = v1 + (double)heightBackScatter / dsBackScatter.RasterYSize * dir.Normalize();
183                 CalpsoRenderObject obj = new CalpsoRenderObject("", v1, v3);
184             }      

View Code

[GDAL]讀取HDF格式的calipso資料
[GDAL]讀取HDF格式的calipso資料
1  public CreateColorRamp CreateColorRamps()
 2         {
 3             CreateColorRamp createColorRamp = new CreateColorRamp();
 4             List<ColorRamp> mColorRamp = new List<ColorRamp>();
 5             double[] valueGray = new double[35] {0.0,   0.0001,0.0002,0.0003,0.0004,0.0005,
 6                                                  0.0006,0.0007,0.0008,0.0009,0.001,
 7                                                  0.0015,0.002, 0.0025,0.003, 0.0035,
 8                                                  0.004, 0.0045,0.005, 0.0055, 0.006,
 9                                                  0.0065,0.007, 0.0075,0.008, 0.01,
10                                                  0.02,  0.03,  0.04,  0.05,  0.06,
11                                                  0.07,  0.08,  0.09,  0.1};
12             int[] Rv = new int[35]             {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
13                                                 255,255,255,255,255,255,255,255,255,255,
14                                                 70,100,130,155,180,220,225,235,240,242,
15                                                 245,249,253,255};
16             int[] Gv = new int[35]             {43,43,43,128,170,213,255,255,255,128,
17                                                 170,255,255,213,170,128,85,0,43,85,128,
18                                                 70,100,130,155,180,220,225,235,240,242,
19                                                 245,249,253,255};
20             int[] Bv = new int[35]             {128,170,170,255,255,255,255,213,170,128,
21                                                 85,0,0,0,0,0,0,0,85,128,170,
22                                                 70,100,130,155,180, 220,225,235,240,242,
23                                                 245,249,253,255};
24 
25             for (int i = 0; i < 35; i++)
26             {
27                 ColorRamp clr = new ColorRamp();
28                 clr.start = valueGray[i];
29                 if (i < 34)
30                 {
31                     clr.end = valueGray[i + 1];
32                 }
33                 clr.ColorS = Color.FromArgb(255, Rv[i], Gv[i], Bv[i]);
34                 mColorRamp.Add(clr);
35             }
36             mColorRamp[mColorRamp.Count - 1].end = 100;
37             createColorRamp.ColorRampCollection = mColorRamp;
38             return createColorRamp;
39         }
40         int GDALProgressFunc(double Complete, IntPtr Message, IntPtr Data)
41         {
42             return 1;
43         }      

CreateColorRamp

[GDAL]讀取HDF格式的calipso資料
[GDAL]讀取HDF格式的calipso資料
1  public class ColorRamp
 2     {
 3         public double start;
 4         public double Interval;
 5         public double end;
 6 
 7         public Color ColorS;
 8     }
 9     public class CreateColorRamp
10     {
11         public List<ColorRamp> ColorRampCollection { get; set; }
12         public Color? GetColor(double val)
13         {
14             if (ColorRampCollection == null)
15             {
16                 throw new ArgumentNullException("顔色分級設定為空");
17             }
18             for (int i = 0; i < ColorRampCollection.Count; i++)
19             {
20                 ColorRamp colorRam = ColorRampCollection[i];
21                 if (val >= colorRam.start && val < colorRam.end)
22                 {
23                     return colorRam.ColorS;
24                 }
25             }
26             return null;
27 
28         }
29     }
30 }      

GDAL C#封裝還是沒有解決中文字元的問題。C# 采用Unicode,GDAL的Dll采用多位元組字元。是以隻能在C++的代碼中做些工作,将Unicode轉換為多位元組字元或者是Utf8。準備測試一下。

NASA上的圖檔:

[GDAL]讀取HDF格式的calipso資料