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

葵花数据投影程序

admin3个月前 (01-20)遥感567

下载的这个链接 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://www.lovedm.club/?id=161

分享给朋友：

返回列表

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

没有最新的文章了...

个人学习笔记

葵花数据投影程序

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

Win10环境下6sV2.1模型编译

GDAL-3.0.1+PROJ-6.2.1+SQLite+VS2019编译

GDAL-3.0.1+PROJ-6.2.1+SQLite+HDF4+HDF5+VS2019编译

利用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.