当前位置：首页 > 遥感 > 正文内容

葵花数据投影程序

admin6个月前 (01-20)遥感1021

下载的这个链接 https://download.csdn.net/download/weixin_42662293/86183693 的代码，原代码是用C++写的，平台无关，将原代码复刻为C#代码，代码如下，可以直接复制创建项目去用。

```csharp
using System;
using System.Collections.Generic;
using System.IO;
using System.Linq;
using System.Runtime.InteropServices;
using System.Text;
using System.Threading.Tasks;

namespace H9Proj
{
    internal class Program
    {
        const double omega_e = 1.00273790934; //solar day over sidereal day
        const double pi = 3.14159265358979323846;
        const double r2d = 180 / pi;
        const double d2r = pi / 180;
        const double Re = 6378.135;
        const double AU = 1.49597870691E8;  /* Astronomical unit - km (IAU 76) */

const double topdx = 0.008333333333333333, topdy = 0.008333333333333333;
        const double topulat = 70.0, topulon = 70.0;
        const int topnx = 8400, topny = 8400;
        const int topnn = topnx * topny;
        static int topflag = 0;
        static short[] topdata = new short[topnn];

struct hsdinfo
        {
            public byte[][] hdr;
            public byte[] headbuf;
            public int headlen;
            public byte[] data;
            public int datalen;
            public int firline;
            public int band;
            public int flag;
            public string file;
        }

[StructLayout(LayoutKind.Sequential, Pack = 1)]
        struct hs_head1
        { //offset=0
            /*0*/
            public byte bnum; //header block number =1
            /*1*/
            public ushort len; //block length =282 bytes
            /*3*/
            public ushort tnum; //total number of header blocks =11
            /*5*/
            public byte byteorder; //0: littel endian; 1: big endian
            /*6*/
            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 16)]
            public char[] satename; //himawari-8
            /*22*/
            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 16)]
            public char[] centername; //
            /*38*/
            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 4)]
            public char[] obarea; // observation area,e.g. R303
            /*42*/
            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 2)]
            public char[] obinfo; //
            /*44*/
            public ushort obtimeline; //hhmm(integer) hh(00-23),mm(00-50)
            /*46*/
            public double starttime; //modified julian date
            /*54*/
            public double endtime; //MJD
            /*62*/
            double filecreatetime; //MJD
            /*70*/
            public uint headlen; //total header length
            /*74*/
            public uint datalen; //total data length
            /*78*/
            public byte qflag1; //operation flag
            /*79*/
            public byte qflag2; //
            /*80*/
            public byte qflag3;
            /*81*/
            public byte qflag4;
            /*82*/
            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 32)]
            public char[] formatvers; //file format version
            /*114*/
            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 128)]
            public char[] filename;
            /*242*/
            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 40)]
            public byte[] spare;
        }

//data information block
        [StructLayout(LayoutKind.Sequential, Pack = 1)]
        struct hs_head2
        { //offset=282(0x11a)
            /*0*/
            public byte bnum; //=2
            /*1*/
            public ushort len; //=50
            /*3*/
            public ushort pixelbits; //=16
            /*5*/
            public ushort cols;
            /*7*/
            public ushort lines;
            /*9*/
            public byte compressflag; //0: no compression(default); 1: gzip; 2: bzip2
            /*10*/
            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 40)]
            public byte[] spare;
        }

//projection information block
        [StructLayout(LayoutKind.Sequential, Pack = 1)]
        struct hs_head3
        { //offset=332(0x14c)
            /*0*/
            public byte bnum; //=3
            /*1*/
            public ushort len; //=127
            /*3*/
            public double sublon;
            /*11*/
            public uint cfac; //column scaling factor
            /*15*/
            public uint lfac; //line scaling factor
            /*19*/
            public float coff; //column offset
            /*23*/
            public float loff; //line offset
            /*27*/
            public double rs; //distance from earth's center to virtual satellite=42164km
            /*35*/
            public double req; //earth's equatorial radius=6378.1370km
            /*43*/
            public double rpol; //earth's polar radius=6356.7523km
            /*51*/
            public double e2; //(req^2-rpol^2)/req^2 =0.00669438444
            /*59*/
            public double ew; //rpol^2/req^2 =0.993305616
            /*67*/
            public double el; //req^2/rpol^2 =1.006739501
            /*75*/
            double sd; //rs^2-req^2 = 1737122264
            /*83*/
            ushort resampletype; //
            /*85*/
            ushort resamplesize; //
            /*87*/
            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 40)]
            public byte[] spare;
        }

//navigation information block
        [StructLayout(LayoutKind.Sequential, Pack = 1)]
        struct hs_head4
        { //offset=459(0x1cb)
            /*0*/
            public byte bnum; //=4
            /*1*/
            public ushort len; //=139
            /*3*/
            public double navtime; //navigation information time [MJD]
            /*11*/
            public double ssplon;
            /*19*/
            public double ssplat;
            /*27*/
            public double rs; //distance from earth's center to satellite[km]
            /*35*/
            public double nadirlon;
            /*43*/
            public double nadirlat;
            /*51*/
            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 3)]
            public double[] sunpos; //[km] (x,y,z)(j2000 inertial coordinate)
            /*75*/
            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 3)]
            public double[] moonpos; //[km] (x,y,z)(j2000 inertial coordinate)
            /*99*/
            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 40)]
            public byte[] spare;
        }

[StructLayout(LayoutKind.Sequential, Pack = 1)]
        public struct ReflectCoef
        {
            /*35*/
            public double cc; //radiance to albedo; A=cc*I; cc[m^2 sr um/W];I[W/(m^2 sr um)]
            /*43*/
            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 104)]
            public byte[] spare; // Reserved space
        }

[StructLayout(LayoutKind.Sequential, Pack = 1)]
        public struct InfraCoef
        {
            /*35*/
            public double c0;
            /*43*/
            public double c1;
            /*51*/
            public double c2; //radiance to brightness T
            /*59*/
            public double cc0;
            /*67*/
            public double cc1;
            /*75*/
            public double cc2; //brightness T to radiance
            /*83*/
            public double c; //speed of light (m/s) 0x2a9
            /*91*/
            public double h; //planck constant (Js)
            /*99*/
            public double k; //boltzmann constant(k)
            /*107*/
            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 40)]
            public byte[] spare; // Reserved space
        }

[StructLayout(LayoutKind.Explicit, Pack = 1)]
        public struct Coefficients
        {
            [FieldOffset(0)]
            public InfraCoef infra;

[FieldOffset(0)]
            public ReflectCoef refect;
        }

//calibration information block
        [StructLayout(LayoutKind.Sequential, Pack = 1)]
        public struct hs_head5
        {//offset=598(0x256)
            /*0*/
            public byte bnum; //=5
            /*1*/
            public ushort len; //=147
            /*3*/
            public ushort bandnum; //1: vis(0.68um) 2: ir4(3.7um) 3: ir3(6.8um) 4: ir1(10.8um) 5: ir2(12.0um)
            /*5*/
            public double wavelen; // central wave length (um)
            /*13*/
            public ushort pixelbits; //valid number of bits per pixel
            /*15*/
            public ushort errvalue; //=65535
            /*17*/
            public ushort outvalue; //=65534
            /*19*/
            public double gain; //radiance= gain*count+constant [W/(m^2 sr um count)]
            /*27*/
            public double constant; //[W/(m^2 sr um)]

[MarshalAs(UnmanagedType.ByValArray, SizeConst = 112)]
            public byte[] other; // Reserved space
        }

//inter-calibration information block
        [StructLayout(LayoutKind.Sequential, Pack = 1)]
        struct hs_head6
        { //offset=745(0x2e9)
            /*0*/
            public byte bnum; //=6
            /*1*/
            public ushort len; //=259
            /*3*/
            public double intercept; //GSICS calibration coefficient(intercept)
            /*11*/
            public double errintercept; //standard error of GSICS calibration coefficient
            /*19*/
            public double slope; //GSICS calibration coefficient(slope)
            /*27*/
            public double errslope; //standard error of GSICS calibration coefficient(slope)
            /*35*/
            public double quadratic; //GSICS calibration coefficient(quadratic term)
            /*43*/
            public double errquadratic; //standard error of GSICS calibration coefficient
            /*51*/
            public double GSICSstarttime;
            /*59*/
            public double GSICSendtime;
            /*67*/
            public float upperlimit; //Radiance valid range of GSICS calibration coefficients
            /*71*/
            public float lowlimit; //
            /*75*/
            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 128)]
            public char[] GSICSfilename;
            /*203*/
            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 56)]
            public char[] spare;
        }

//segment information block
        [StructLayout(LayoutKind.Sequential, Pack = 1)]
        struct hs_head7
        { //offset=1004(0x3ec)
            /*0*/

public byte bnum; //=7
            /*1*/
            public ushort len; //=47
            /*3*/
            public byte segments; //total number of segments. 1: no devision
            /*4*/
            public byte segseqnum; //segment sequence number
            /*5*/
            public ushort firstline; //first line number of image segment
            /*7*/
            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 40)]
            public char[] spare;
        }

//navigation correction information block
        [StructLayout(LayoutKind.Sequential, Pack = 1)]
        struct hs_head8
        { //offset=1051(0x41b)
            /*0*/
            public byte bnum; //=8
            /*1*/
            public ushort len; //
            /*3*/
            public float centcol; //center column of rotation
            /*7*/
            public float centline; //center line of rotation
            /*11*/
            public double rotamount; //amount of rotational correction
            /*19*/
            public ushort corrnum; //number of correction information data for column and line direction
            /*21*/
            public ushort rotline; //line number after rotation
            /*23*/
            public float shiftcol; //shift amount for column direction
            /*27*/
            public float shiftline; //shift amount for line direction
            /*corrnum repeats of rotline,shiftcol,shiftline*/
            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 40)]
            public byte[] spare;
        }

//observation time information block
        [StructLayout(LayoutKind.Sequential, Pack = 1)]
        struct hs_head9
        {
            /*0*/
            public byte bnum; //=9
            /*1*/
            public ushort len;
            /*3*/
            public ushort obstimes; // number of observation times
            /*5*/
            public ushort line;
            /*7*/
            public double obstime; //obervation time
            /*  obstimes repeats of line and obstime */
            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 40)]
            public byte[] spare;
        }

//error information block
        [StructLayout(LayoutKind.Sequential, Pack = 1)]
        struct hs_head10
        {
            /*0*/
            public byte bnum; //=10
            /*1*/
            public ushort len;
            /*3*/
            public ushort errnum; //number of error information data
            /*5*/
            public ushort line; //line number
            /*7*/
            public ushort errpixels; //number of error pixels per line
            /* errnum repeats of  line errpixels*/
            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 40)]
            public byte[] spare;
        }

//spare block
        [StructLayout(LayoutKind.Sequential, Pack = 1)]
        struct hs_head11
        {
            /*0*/
            public byte bnum; //=11
            /*1*/
            public ushort len; //=259
            /*3*/
            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 256)]
            public byte[] spare;
        }

static void Main(string[] args)
        {
            string directory = AppDomain.CurrentDomain.BaseDirectory;
            Console.WriteLine(directory);

string[] RES = { "R10", "R10", "R05", "R10", "R20", "R20", "R20", "R20", "R20", "R20", "R20", "R20", "R20", "R20", "R20", "R20" };
            string[] bandNames = {"band1:0.46um", "band2:0.51um", "band3:0.64um", "band4:0.86um", "band5:1.6um", "band6:2.3um",
                                  "band7:3.9um", "band8:6.2um", "band9:7.0um", "band10:7.3um", "band11:8.6um", "band12:9.6um",
                                  "band13:10.4um","band14:11.2um","band15:12.3um","band16:13.3um"};

string file = null, dstfile = null, pstr, p;
            int firline, si, iio, segnum;
            double[] adfout = { 70, 0.01, 0, 60, 0, -0.01 };
            int nx = 8000, ny = 6000;
            float[] dstdata = null;
            int sateflag = 0, sunflag = 0, dnflag = 0;

string hdrfile = "", outfile = "z.img";

List<string> bands = new List<string>();	//选中的波段

if (args.Length == 0)
            {
                Console.WriteLine("<-dstwin lux luy dx dy> <-size nx ny> <-bands 1,2,3> <-addsat> <-addsun> <-dn> srcfile <dstfile>");
                return;
            }

//dealing with arguments
            for (int i = 0; i < args.Length; i++)
            {
                Console.WriteLine($"argv: {i}, {args[i]}");
                if (args[i] == "-dstwin" && i < args.Length - 3)
                {
                    adfout[0] = double.Parse(args[++i]); //LU corner x
                    adfout[3] = double.Parse(args[++i]); //LU corner y
                    adfout[1] = double.Parse(args[++i]); //dx
                    adfout[5] = double.Parse(args[++i]); //dy
                }
                else if (args[i] == "-size" && i < args.Length - 1)
                {
                    nx = int.Parse(args[++i]);
                    ny = int.Parse(args[++i]);
                }
                else if (args[i] == "-bands" && i < args.Length)
                {
                    string str_bands = args[++i];
                    bands = str_bands.Split(',').ToList();
                }
                else if (args[i] == "-addsat")
                {
                    sateflag = 1;
                }
                else if (args[i] == "-addsun")
                {
                    sunflag = 1;
                }
                else if (args[i] == "-dn")
                {
                    dnflag = 1;
                }
                else if (file == null)
                {
                    file = args[i];
                }
                else if (dstfile == null)
                {
                    dstfile = args[i];
                }
                else
                {
                    Console.WriteLine("<-dstwin lux luy dx dy> <-size nx ny> <-bands 1,2,3> srcfile <dstfile>\n");
                    return;
                }
            }

//dealing with outfile name
            if (dstfile == null)
            {
                dstfile = outfile;
            }
            int ii = 0;

hdrfile = Path.GetDirectoryName(dstfile) + "\\" + Path.GetFileNameWithoutExtension(dstfile) + ".hdr";
            Console.WriteLine($"outfile: {dstfile}\nhdrfile: {hdrfile}");

//没选择波段时：
            if (bands.Count < 1)
            {
                string fileName = Path.GetFileName(file);
                bands.Add(fileName.Substring(22, 2));
            }

FileStream fp = new FileStream(dstfile, FileMode.OpenOrCreate);
            if (fp == null)
            {
                Console.WriteLine($"error on open {dstfile}");
                return;
            }

string hdr_bandNames = "band names = {";
            int l = bands.Count;
            for (int i = 0; i < l; i++)
            {
                int bandId = int.Parse(bands[i]);
                string char_res = RES[bandId - 1];
                string char_bandName;
                char_bandName = string.Format("B{0:D2}", bandId);
                string bandfile = getbandfile(file, char_bandName, char_res);
                Console.WriteLine(bandfile);
                float[] bandData = projBandFile(bandfile, nx, ny, adfout, dnflag);
                byte[] bandDataByte = new byte[bandData.Length * 4];
                Buffer.BlockCopy(bandData, 0, bandDataByte, 0, bandData.Length * 4);
                Console.WriteLine("start");
                for (int j = 0; j < ny; j++)
                {

fp.Write(bandDataByte, j * nx * 4, nx * 4);
                    //fwrite(bandData + j * nx, sizeof(float), nx, fp);
                }
                //delete[] bandData;
                Console.WriteLine("check");
                hdr_bandNames += bandNames[bandId - 1];
                if (i < l - 1)
                    hdr_bandNames += ",";
            }

//投影卫星天顶角
            if (sateflag == 1)
            {
                float[] sate = null, sata = null;
                int ret = projSatZenith(file, nx, ny, adfout, ref sate, ref sata);
                //float *bandData=NULL;
                if (ret == 0)
                {
                    byte[] sateByte = new byte[sate.Length * 4];
                    Buffer.BlockCopy(sate, 0, sateByte, 0, sate.Length * 4);
                    for (int i = 0; i < ny; i++)
                    {
                        fp.Write(sateByte, i * nx * 4, nx * 4);
                    }
                    //GC.Collect();
                    //byte[] sataByte = new byte[sata.Length * 4];
                    //Buffer.BlockCopy(sata, 0, sataByte, 0, sata.Length * 4);
                    //for (int i = 0; i < ny; i++)
                    //{
                    //    fp.Write(sataByte, i * nx * 4, nx * 4);
                    //}
                    if (sate != null) sate = null;
                    //if (sata != null) sata = null;

//hdr_bandNames += ",satellite_zenith,satellite_azimuth";
                    hdr_bandNames += ",satellite_zenith";
                    l += 1;
                    GC.Collect();
                }
            }

//投影太阳高度角
            if (sunflag == 1)
            {
                float[] sune = null, suna = null;
                //float* sunData = SunElevation(file, nx, ny, adfout);
                int ret = SunElevation(file, nx, ny, adfout, ref sune, ref suna);
                //float *bandData=NULL;
                if (ret == 0)
                {
                    byte[] suneByte = new byte[sune.Length * 4];
                    Buffer.BlockCopy(sune, 0, suneByte, 0, sune.Length * 4);
                    for (int i = 0; i < ny; i++)
                    {
                        fp.Write(suneByte, i * nx * 4, nx * 4);
                    }

//byte[] sunaByte = new byte[suna.Length * 4];
                    //Buffer.BlockCopy(suna, 0, sunaByte, 0, suna.Length * 4);
                    //for (int i = 0; i < ny; i++)
                    //{
                    //    fp.Write(sunaByte, i * nx * 4, nx * 4);
                    //}
                    if (sune != null) sune = null;
                    //if (suna != null) suna = null;

//hdr_bandNames += ",sun_elevation,sun_azimuth";
                    hdr_bandNames += ",sun_elevation";
                    l += 1;
                }
            }

hdr_bandNames += "}\n";
            fp.Close();

fp = new FileStream(hdrfile, FileMode.OpenOrCreate);
            byte[] str = Encoding.Default.GetBytes("ENVI\n");
            fp.Write(str, 0, str.Length);
            str = Encoding.Default.GetBytes("description = {\nthis file using parameter as following:\n");
            fp.Write(str, 0, str.Length);
            str = Encoding.Default.GetBytes("}\n");
            fp.Write(str, 0, str.Length);
            str = Encoding.Default.GetBytes($"samples = {nx}\n");
            fp.Write(str, 0, str.Length);
            str = Encoding.Default.GetBytes($"lines   = {ny}\n");
            fp.Write(str, 0, str.Length);
            str = Encoding.Default.GetBytes($"bands   = {l}\nheader offset = 0\nfile type = ENVI Standard\n");
            fp.Write(str, 0, str.Length);
            str = Encoding.Default.GetBytes($"data type = 4\ninterleave = bsq\nbyte order = 0\n");
            fp.Write(str, 0, str.Length);
            str = Encoding.Default.GetBytes($"{hdr_bandNames}\n");
            fp.Write(str, 0, str.Length);
            str = Encoding.Default.GetBytes($"data ignore = 0\n");
            fp.Write(str, 0, str.Length);
            str = Encoding.Default.GetBytes("map info = {Geographic Lat/Lon, ");
            fp.Write(str, 0, str.Length);
            str = Encoding.Default.GetBytes($"1, 1, {adfout[0]}, {adfout[3]}, ");
            fp.Write(str, 0, str.Length);
            str = Encoding.Default.GetBytes($"{Math.Abs(adfout[1])}, {Math.Abs(adfout[5])}, WGS-84, units=Degrees}}\n");
            fp.Write(str, 0, str.Length);
            fp.Close();

Console.ReadKey();

}

static string getbandfile(string file, string band, string res)
        {

int len = file.Length;
            string srcfile = file;

string str;
            string path = Path.GetDirectoryName(srcfile);
            string name = Path.GetFileName(srcfile);
            int id = name.LastIndexOf("B");
            name = name.Remove(id, 3);
            name = name.Insert(id, band);
            id = name.LastIndexOf("R");
            name = name.Remove(id, 3);
            name = name.Insert(id, res);

str = path + "\\" + name;

return str;
        }

/*单波段投影*/
        static float[] projBandFile(string file, int nx, int ny, double[] adfout, int dnflag = 0)
        {

int line, col, lines, cols;
            string hdrfile = "", infile = "";
            int si, segnum, band, ii = 0;
            hsdinfo[] hsd;
            hsdinfo phsd;
            double lon, lat, cx = 0, ly = 0;
            //get segment order and total number from infile name
            infile = file;
            int pstr = infile.LastIndexOf('.');
            double h;

if (pstr != -1)
            {
                pstr -= 4; //locate pstr at first k of Skkll
                si = int.Parse(infile.Substring(pstr, 2));
                segnum = int.Parse(infile.Substring(pstr + 2, 2));
                //sscanf(pstr, "%2d%2d", &si, &segnum);
                //fprintf(stderr,"segnum=%d,si=%d\n",si,segnum);
            }
            else
            {
                segnum = 1;
                si = 1;
            }

if (segnum > 0 && segnum <= 10 && si > 0 && si <= segnum)
            {
                hsd = new hsdinfo[segnum];
                for (int i = 0; i < segnum; i++)
                {
                    hsd[i].hdr = new byte[11][];
                    hsd[i].headbuf = null;
                    hsd[i].data = null;
                    hsd[i].flag = 0;
                }
            }
            else
            {
                Console.WriteLine("wrong segnum or si.");
                return null;
            }

//read first hsd file and get lines and columns parameter         
            hsd[si - 1].file = infile;
            readhsdfile(ref hsd[si - 1]);
            phsd = hsd[si - 1];

//GCHandle handle = GCHandle.Alloc(hsd[0].hdr[0], GCHandleType.Pinned);
            //hs_head1 myStruct = (hs_head1)Marshal.PtrToStructure(handle.AddrOfPinnedObject(), typeof(hs_head1));
            //handle.Free();

GCHandle handle = GCHandle.Alloc(phsd.hdr[6], GCHandleType.Pinned);
            hs_head7 hs_head7 = (hs_head7)Marshal.PtrToStructure(handle.AddrOfPinnedObject(), typeof(hs_head7));
            handle.Free();
            handle = GCHandle.Alloc(phsd.hdr[4], GCHandleType.Pinned);
            hs_head5 hs_head5 = (hs_head5)Marshal.PtrToStructure(handle.AddrOfPinnedObject(), typeof(hs_head5));
            handle.Free();
            handle = GCHandle.Alloc(phsd.hdr[1], GCHandleType.Pinned);
            hs_head2 hs_head2 = (hs_head2)Marshal.PtrToStructure(handle.AddrOfPinnedObject(), typeof(hs_head2));
            handle.Free();
            handle = GCHandle.Alloc(phsd.hdr[2], GCHandleType.Pinned);
            hs_head3 hs_head3 = (hs_head3)Marshal.PtrToStructure(handle.AddrOfPinnedObject(), typeof(hs_head3));
            handle.Free();
            handle = GCHandle.Alloc(phsd.hdr[4], GCHandleType.Pinned);
            hs_head5 hs_head4 = (hs_head5)Marshal.PtrToStructure(handle.AddrOfPinnedObject(), typeof(hs_head5));
            handle.Free();
            handle = GCHandle.Alloc(hs_head5.other, GCHandleType.Pinned);
            ReflectCoef reflectCoef = (ReflectCoef)Marshal.PtrToStructure(handle.AddrOfPinnedObject(), typeof(ReflectCoef));
            InfraCoef infraCoef = (InfraCoef)Marshal.PtrToStructure(handle.AddrOfPinnedObject(), typeof(InfraCoef));
            handle.Free();

lines = hs_head2.lines; //*(unsigned short *)(phsd->hdr[1]+7);
            cols = hs_head2.cols; //*(unsigned short *)(phsd->hdr[1]+5);
            hsd[si - 1].firline = hs_head7.firstline - 1; //*(unsigned short *)(phsd->hdr[6]+5)-1;
            band = hs_head5.bandnum; //*(unsigned short*)(phsd->hdr[4]+3);

float[] dstdata = new float[nx * ny];
            if (dstdata == null)
            {
                Console.WriteLine("can't malloc enough space for dstdata.");
                return null;
            }

//do projection and calculate albedo or brightness temperature
            for (int i = 0; i < ny; i++)
            {
                for (int j = 0; j < nx; j++)
                {
                    lon = adfout[0] + adfout[1] * j;
                    lat = adfout[3] + adfout[5] * i;
                    h = getHeight(@"F:\CPlusPlusCode\Himawari8Proj\china.img", lon, lat);
                    ll2cl(lon, lat, ref cx, ref ly, hs_head3, h);
                    col = (int)(cx); line = (int)(ly);
                    ii = line / lines + 1;
                    ////fprintf(stdout,"segment %d.\n",ii);
                    if (ii > 0 && ii <= segnum)
                    {
                        if (ii != si && hsd[ii - 1].flag == 0)
                        {
                            infile = infile.Remove(pstr, 4);
                            infile = infile.Insert(pstr, string.Format("{0:D2}", ii) + string.Format("{0:D2}", segnum));
                            hsd[ii - 1].file = infile;
                            //strcpy(hsd[ii - 1].file, infile);
                            //fprintf(stderr,"%s\n",infile);
                            if (readhsdfile(ref hsd[ii - 1]) == 0)
                            {
                                handle = GCHandle.Alloc(hsd[ii - 1].hdr[6], GCHandleType.Pinned);
                                hs_head7 = (hs_head7)Marshal.PtrToStructure(handle.AddrOfPinnedObject(), typeof(hs_head7));
                                handle.Free();

hsd[ii - 1].firline = hs_head7.firstline - 1; //*(unsigned short *)(hsd[ii-1].hdr[6]+5)-1;
                                si = ii;
                            }
                        }
                    }

line -= hsd[ii - 1].firline;
                    if (col >= 0 && col < cols && line >= 0 && line < lines && hsd[ii - 1].data != null)
                    {
                        //double val = (double)*(short*)(hsd[ii - 1].data + (col + line * cols) * 2);
                        double val = BitConverter.ToInt16(hsd[ii - 1].data, (col + line * cols) * 2);
                        if (dnflag != 0) { dstdata[i * nx + j] = (float)val; }
                        else
                        {
                            if (band > 0 && band < 7)
                                dstdata[i * nx + j] = (float)(albedo(val, hs_head5, reflectCoef)/**1000*/);
                            else
                                dstdata[i * nx + j] = (float)(brightT(val, hs_head5, infraCoef)/**10*/);
                        }
                        if (i == ny / 2 && j == nx / 2)
                        {
                            Console.WriteLine($"band={band}");
                            albedo(val, hs_head5, reflectCoef, 1);
                        }
                    }
                }
            }

for (int i = 0; i < segnum; i++)
            {
                if (hsd[i].data != null) hsd[i].data = null;
                if (hsd[i].headbuf != null) hsd[i].headbuf = null;
            }
            //
            //if (hsd) delete[] hsd;
            return dstdata;

// return null;
        }

static int readhsdfile(ref hsdinfo hsd)
        {
            FileStream fp = new FileStream(hsd.file, FileMode.Open);
            short totalheads, i;
            byte[] p;

// read data from hsd file
            hsd.flag = 1;
            if (fp == null)
            {
                Console.WriteLine($"can't open %s {hsd.file}");
                return -1;
            }
            Console.WriteLine($"opened {hsd.file}.");
            byte[] buf = new byte[2];
            fp.Seek(3, SeekOrigin.Begin);
            fp.Read(buf, 0, sizeof(short));
            totalheads = BitConverter.ToInt16(buf, 0);
            //Console.WriteLine($"totalheads: {totalheads}");

buf = new byte[4];
            fp.Seek(70, SeekOrigin.Begin);
            fp.Read(buf, 0, sizeof(int));
            hsd.headlen = BitConverter.ToInt32(buf, 0);
            fp.Read(buf, 0, sizeof(int));
            hsd.datalen = BitConverter.ToInt32(buf, 0);

if (hsd.headbuf != null)
            {
                hsd.headbuf = null;
            }
            hsd.headbuf = new byte[hsd.headlen];

if (hsd.data != null)
            {
                hsd.data = null;
            }
            hsd.data = new byte[hsd.datalen];

if (hsd.data == null)
            {
                Console.WriteLine("can't malloc enough space.");
                return -1;
            }

fp.Seek(0, SeekOrigin.Begin);

if (hsd.headbuf != null && hsd.data != null)
            {
                fp.Read(hsd.headbuf, 0, hsd.headlen);
                fp.Read(hsd.data, 0, hsd.datalen);

p = hsd.headbuf;
                int pointIndex = 0;
                for (i = 0; i < totalheads; i++)
                {
                    short blocksize = BitConverter.ToInt16(hsd.headbuf, pointIndex + 1);
                    byte[] data = new byte[blocksize];
                    Array.Copy(hsd.headbuf, pointIndex, data, 0, blocksize);
                    hsd.hdr[i] = data;
                    // Console.WriteLine(((byte)BitConverter.ToChar(hsd.headbuf, pointIndex)));
                    pointIndex += blocksize;

//byte[] headerblock = hsd.headbuf.
                }
            }
            fp.Close();
            return 0;
        }

static double getHeight(string file, double lon, double lat)
        {

int i, j;
            double h = 0;

if (topflag == 0)
            {
                using (FileStream fp = new FileStream(file, FileMode.Open))
                {
                    if (fp == null)
                    {
                        Console.WriteLine($"can't open %s {file}");
                        topflag = -1;
                        return 0;
                    }
                    byte[] data = new byte[topnn * sizeof(short)];
                    fp.Read(data, 0, topnn * sizeof(short));
                    Buffer.BlockCopy(data, 0, topdata, 0, data.Length);
                    fp.Close();
                    topflag = 1;
                }
            }
            if (topflag == -1) return 0;
            j = (int)((lon - topulon) / topdx + 0.5);
            i = (int)((topulat - lat) / topdy + 0.5);
            if (i >= 0 && i < 8400 && j >= 0 && j < 8400)
            {
                if (topdata[i * 8400 + j] > 0 && topdata[i * 8400 + j] < 10000) h = topdata[i * 8400 + j] / 1000.0;
            }
            return h;

}

// geographical (longitude,latitude)->satellite scan's (X direction Angle,Z direction angle)
        // ->image's (column,line)
        static int ll2cl(double lon, double lat, ref double c, ref double l, hs_head3 header3, double h = 0)
        {

hs_head3 hd3 = header3;
            double d2r = 0.01745329252;
            double r2d = 57.295779513;
            double clat = Math.Atan(hd3.ew * Math.Tan(lat * d2r));
            //double sublon=*(double*)(hsd.hdr[2]+3);
            double dellon = d2r * (lon - hd3.sublon);
            double rp = hd3.rpol / (Math.Sqrt(1 - hd3.e2 * Math.Cos(clat) * Math.Cos(clat))) + h;
            double r1 = hd3.rs - rp * Math.Cos(clat) * Math.Cos(dellon);
            double r2 = -rp * Math.Cos(clat) * Math.Sin(dellon);
            double r3 = rp * Math.Sin(clat);
            double rn = Math.Sqrt(r1 * r1 + r2 * r2 + r3 * r3);
            double x = Math.Atan2(-r2, r1) * r2d;
            double y = Math.Asin(-r3 / rn) * r2d;

//double cfac=*(unsigned int*)(hsd.hdr[2]+11);
            //double lfac=*(unsigned int*)(hsd.hdr[2]+15);
            //double coff=*(float*)(hsd.hdr[2]+19);
            //double loff=*(float*)(hsd.hdr[2]+23);

c = hd3.coff + x * hd3.cfac * Math.Pow(2, -16); //round(x*pow(2,-16)*cfac);
            l = hd3.loff + y * hd3.lfac * Math.Pow(2, -16); //round(y*pow(2,-16)*lfac);

return 0;
        }

//count->radiance->aldebo
        //radiance = gain * count + constant
        //albedo = c0 * radiance
        static double albedo(double val, hs_head5 hd5, ReflectCoef reflectCoef, int ifag = 0)
        {

double rad, te, reflec;
            rad = hd5.gain * val + hd5.constant;

reflec = rad * reflectCoef.cc;
            if (ifag == 1)
            {
                Console.WriteLine($"count: {val}\tgain: {hd5.gain}\tconstant: {hd5.constant}\tc': {reflectCoef.cc}");
            }
            return reflec;
        }

//count->radiance->effective brightness->brightness temperature
        static double brightT(double val, hs_head5 hd5, InfraCoef infraCoef)
        {

double c0 = infraCoef.c0; //*(double*)(hsd.hdr[4]+35);
            double c1 = infraCoef.c1; //*(double*)(hsd.hdr[4]+43);
            double c2 = infraCoef.c2; //*(double*)(hsd.hdr[4]+51);
            double c = infraCoef.c; //*(double*)(hsd.hdr[4]+83);
            double h = infraCoef.h; //*(double*)(hsd.hdr[4]+91);
            double k = infraCoef.k; //*(double*)(hsd.hdr[4]+99);
            double rad, te;
            rad = hd5.gain * val + hd5.constant;
            //double te1=1.43883e4/wv/log(1+1.19107e8/(rad*pow(wv,5)));
            double wv = hd5.wavelen * 1e-6; //um->m
            rad = rad * 1e6;//um->m
                            //fprintf(stderr,"c=%lf,h=%lg,k=%lg,wv=%lg\n",c,h,k,wv);
                            //fprintf(stderr,"c0=%lf,c1=%lf,c2=%lf\n",c0,c1,c2);

te = h * c / k / wv / (Math.Log(2 * h * c * c / (Math.Pow(wv, 5) * rad) + 1));
            double bt = c0 + c1 * te + c2 * te * te;
            //fprintf(stderr,"gain=%lf,cont=%lf,val=%lf\nrad=%lf,te=%lf,bt=%lf\n",gain,cont,val,rad,te,bt);
            return bt;

}

/*卫星高度角投影；*/
        static int projSatZenith(string file, int nx, int ny, double[] adfout, ref float[] sate, ref float[] sata)
        {
            int i, j;

hsdinfo hsd;
            double lat, lon;
            double range = 0, azim = 0, elev = 0;

hsd = new hsdinfo();
            hsd.headbuf = null;
            hsd.data = null;
            hsd.flag = 0;
            hsd.file = file;
            hsd.hdr = new byte[11][];
            readhsdfile(ref hsd);

//float* dstdata=new float[nx*ny];
            if (sate == null) sate = new float[nx * ny];
            if (sata == null) sata = new float[nx * ny];
            if (sata == null || sate == null)
            {
                Console.WriteLine("error on allocating array satelev or satazim data.");
                return -1;
            }

var handle = GCHandle.Alloc(hsd.hdr[2], GCHandleType.Pinned);
            hs_head3 hd3 = (hs_head3)Marshal.PtrToStructure(handle.AddrOfPinnedObject(), typeof(hs_head3));
            handle.Free();

for (i = 0; i < ny; i++)
            {
                lat = adfout[3] + adfout[5] * i;
                for (j = 0; j < nx; j++)
                {
                    lon = adfout[0] + adfout[1] * j;
                    //dstdata[i*nx+j]=satZenith(lon,lat,*hsd);
                    obslook_sat(ref elev, ref azim, ref range, lon, lat, hd3, 0);
                    sate[i * nx + j] = (float)elev;
                    sata[i * nx + j] = (float)azim;
                }
            }

return 0;
        }

//calculating sat's elevation and azimuth angles
        static int obslook_sat(ref double el, ref double az, ref double rg, double lon, double lat, hs_head3 hd3, double h = 0)
        {
            double d2r = 0.01745329252;
            double r2d = 57.295779513;
            double x0, y0, z0;
            double rx, ry, rz;
            double xs, ys, zs = 0;
            double ls, le, lz;
            double clat = Math.Atan(0.993305616 * Math.Tan(lat * d2r));
            double clon = d2r * lon;

//double sublon=*(double*)(hsd.hdr[2]+3);
            double slon = d2r * hd3.sublon;
            double rp = 6356.7523 / (Math.Sqrt(1 - 0.00669438444 * Math.Cos(clat) * Math.Cos(clat))) + h;

x0 = rp * Math.Cos(clat) * Math.Cos(clon);
            y0 = rp * Math.Cos(clat) * Math.Sin(clon);
            z0 = rp * Math.Sin(clat);
            xs = 42164 * Math.Cos(slon);
            ys = 42164 * Math.Sin(slon);
            zs = 0;
            rx = xs - x0; ry = ys - y0; rz = zs - z0;

//return (acos(cos(clat)*cos(dellon))+acos(r1/rg))*r2d;

ls = Math.Sin(lat * d2r) * Math.Cos(clon) * rx + Math.Sin(lat * d2r) * Math.Sin(clon) * ry - Math.Cos(lat * d2r) * rz;
            le = -Math.Sin(clon) * rx + Math.Cos(clon) * ry;
            lz = Math.Cos(lat * d2r) * Math.Cos(clon) * rx + Math.Cos(lat * d2r) * Math.Sin(clon) * ry + Math.Sin(lat * d2r) * rz;
            az = Math.Atan(-le / ls);
            if (ls > 0) az += pi;
            if (az < 0) az += 2 * pi;
            rg = Math.Sqrt(rx * rx + ry * ry + rz * rz);
            el = Math.Asin(lz / rg);
            el = el * r2d;
            az = az * r2d;
            return 0;
        }

/*太阳仰角数据层*/
        static int SunElevation(string file, int nx, int ny, double[] adfout, ref float[] sune, ref float[] suna)
        {
            int i, j;

hsdinfo hsd;
            double lat, lon;
            double range = 0, azim = 0, elev = 0, jd, x, y, z;
            hs_head1 hs1;
            hs_head4 hs4;
            hsd = new hsdinfo();
            hsd.headbuf = null;
            hsd.data = null;
            hsd.flag = 0;
            hsd.file = file;
            hsd.hdr = new byte[11][];
            readhsdfile(ref hsd);
            var handle = GCHandle.Alloc(hsd.hdr[0], GCHandleType.Pinned);
            hs1 = (hs_head1)Marshal.PtrToStructure(handle.AddrOfPinnedObject(), typeof(hs_head1));
            handle.Free();
            handle = GCHandle.Alloc(hsd.hdr[3], GCHandleType.Pinned);
            hs4 = (hs_head4)Marshal.PtrToStructure(handle.AddrOfPinnedObject(), typeof(hs_head4));
            handle.Free();

//https://en.wikipedia.org/wiki/Julian_day(mjd为0h Nov 17, 1858（即2400000.5起算的jd数)
            jd = hs1.starttime + 2400000.5;
            x = hs4.sunpos[0]; y = hs4.sunpos[1]; z = hs4.sunpos[2];
            //obslook_sphere(jd,x,y,z,120,30,0,azim,elev,range);
            //fprintf(stderr,"(120,30,0)->(azimuth,elevation):(%f,%f)\n",azim,elev);
            if (sune == null) sune = new float[nx * ny];
            if (suna == null) suna = new float[nx * ny];
            if (suna == null || sune == null)
            {
                Console.WriteLine("error on allocating array sunelev or sunazim data.");
                return -1;
            }

for (i = 0; i < ny; i++)
            {
                lat = adfout[3] + adfout[5] * i;
                for (j = 0; j < nx; j++)
                {
                    lon = adfout[0] + adfout[1] * j;
                    obslook_sphere(jd, x, y, z, lon, lat, 0, ref azim, ref elev, ref range);
                    sune[i * nx + j] = (float)elev;
                    suna[i * nx + j] = (float)azim;
                }
            }
            //writeldf(dstfile,hdrfile,dstdata,adfout,nx,ny);

return 0;
        }

// see "orbital coordinate systems Part II" by Kelso
        /* If the satellite's position is the ECI coordinate system is defined as 
         [xs,ys,zs] and the observer is [x0,y0,z0],then the range vector is simply
         [rx,ry,rz]=[xs-x0,ys-y0,zs-z0]. This vector is the ECI system and to generate
         look angles, we need the topocentric-horizon system. That system has its z axis
         pointing toward the zenith, the x axis pointing south, and the y axis pointing
         east.
         it is same as the function Calculate_Obs() adopted from predict.c but calculating obs(x,y,z)
        */
        static int obslook_sphere(double jd, double xs, double ys, double zs,
                           double lon, double lat, double alt,
                           ref double az, ref double el, ref double rg)
        {
            double x0 = 0, y0 = 0, z0 = 0;
            double rx, ry, rz;
            double theta, ls, le, lz;
            ll2xyz_sphere(jd, lon, lat, alt, ref x0, ref y0, ref z0);
            theta = fmod(thetaG_JD(jd) + lon * d2r, 2 * pi);
            rx = xs - x0; ry = ys - y0; rz = zs - z0;
            ls = Math.Sin(lat * d2r) * Math.Cos(theta) * rx + Math.Sin(lat * d2r) * Math.Sin(theta) * ry - Math.Cos(lat * d2r) * rz;
            le = -Math.Sin(theta) * rx + Math.Cos(theta) * ry;
            lz = Math.Cos(lat * d2r) * Math.Cos(theta) * rx + Math.Cos(lat * d2r) * Math.Sin(theta) * ry + Math.Sin(lat * d2r) * rz;
            az = Math.Atan(-le / ls);
            if (ls > 0) az += pi;
            if (az < 0) az += 2 * pi;
            rg = Math.Sqrt(rx * rx + ry * ry + rz * rz);
            el = Math.Asin(lz / rg);
            el = el * r2d;
            az = az * r2d;
            return 0;
        }

// (longitude,latitude,altude,jd) to (x,y,z) of ECI coordinate on sphere earth
        // see "orbital coordinate systems Part II" by Kelso
        // Reference: the 1992 Astronomical Almanac,Page K11
        // x,y,z unit (km); lon,lat unit degree; alt unit m;
        static int ll2xyz_sphere(double jd, double lon, double lat, double alt, ref double x, ref double y, ref double z)
        {
            double theta = fmod(thetaG_JD(jd) + lon * d2r, 2 * pi);
            double r = (Re + alt / 1000) * Math.Cos(lat * d2r);
            x = r * Math.Cos(theta);
            y = r * Math.Sin(theta);
            z = (Re + alt / 1000) * Math.Sin(lat * d2r);
            return 0;
        }

//julian_date to mean sidereal time at Greenwich
        //see "orbital coordinate systems Part II" by Kelso
        //Reference: The 1992 Astronomical Almanac,page B6
        static double thetaG_JD(double jd)
        {
            double jd0 = 0;
            double ut = modf(jd + 0.5, ref jd0);
            jd0 = jd - ut;
            double T = (jd0 - 2451545.0) / 36525;
            double gmst = 24110.54841 + T * (8640184.812866 + T * (0.093104 - T * 6.2E-6));
            gmst = fmod(gmst + omega_e * 86400.0 * ut, 86400);
            if (gmst < 0) gmst += 86400;
            return gmst * 2 * pi / 86400.0;
        }

public static double fmod(double x, double y)
        {
            double result = Math.IEEERemainder(x, y);
            if (result != 0 && (result < 0) ^ (y < 0))
            {
                result += y;
            }
            return result;
        }

public static double modf(double x, ref double integerPart)
        {
            // 使用数学运算来获取整数部分和小数部分
            integerPart = Math.Truncate(x); // 或者使用 (int)x 强制转换，但这对于负数可能不是你所期望的行为
            return x - integerPart;
        }

}
}

```

扫描二维码推送至手机访问。

本文链接：https://lovedm.club/?id=161

分享给朋友：

返回列表

上一篇：Python GDAL In GDAL 4.0, exceptions will be enabled by default

没有最新的文章了...

“葵花数据投影程序” 的相关文章

个人学习笔记

葵花数据投影程序

“葵花数据投影程序” 的相关文章

AVIRIS数据下载与预处理相关

ENVI大气校正

利用GDAL进行波段合成

国产卫星绝对辐射定标系数（2008——2022）

Python GDAL In GDAL 4.0, exceptions will be enabled by default

鲁ICP备18000987号-1 鲁公网安备 37148202000241号

Powered By Z-BlogPHP. Theme by TOYEAN.