雷达基数据中儒略日换算成公历

Lauretta

自己根据SA雷达基数据说明摸索了下，参考朋友的C语言编写的转换程序，写了用matlab读取雷达基数据中儒略日，然后转成公历的方法。希望有用~
clear;clc;         
fid=fopen('D:\data\Z_RADR_I_Z9250_20090605110000_O_DOR_SA_CAP.bin','rb','b');     %2009年6月5日11时雷达基数据   
[rda_data,COUNT]=fread(fid,'ubit8');  
fclose(fid);
rda_data=reshape(rda_data,2432,COUNT/2432);  %   数据格式是一列2432个字节，一列为一条径向数据
A=rda_data(1:2432,1:100);
juliandata=(rda_data(34,1))*256+rda_data(33,1);   %读取儒略日数据
hms=(rda_data(32,:))*256*256*256+(rda_data(31,:))*256*256+(rda_data(30,:))*256+rda_data(29,:);       %读取时间数据
jd=juliandata;
jd=jd+2440587;jd=double(int32(jd));
L=jd+68569;L=double(int32(L));
N=(4*L)/146097;N=double(int32(N));
L=L-(146097*N+3)/4;
I=4000*(L+1)/1461001;I=double(int32(I));
L=L-(1461*I)/4+31;L=double(int32(L));
J=80*L/2447;J=double(int32(J));
day=L-fix(2447*J/80);
L=J/11;L=double(int32(L));
month=J+2-12*L;
year=100*(N-49)+I+L;
%day=double(int32(day));
month=double(int32(month));
year=double(int32(year));
seconds=hms/1000;seconds=double(int32(seconds));
h=seconds/3600;h=double(int32(h));
m=floor((seconds-h*3600)/60);
s=seconds-(60*h+m)*60;

mofangbao

感谢分享，雷达数据转换要熟悉还得下点时间才行

topmad

这个算法好复杂....
我就喜欢DAY OF YEAR

Lauretta

topmad 发表于 2012-10-25 14:27

                               
登录/注册后可看大图

这个算法好复杂....
我就喜欢DAY OF YEAR

一年的第几天啊。。。不过SA雷达数据说明里的那个儒略日是从1970年1月1日开始，显示出来得写成公历才行。有个fortran写的计算日期是一年第几天的程序，matlab没写过。。
program main
        implicit none
        integer iyear,iDDD,imonth,iday       
        character*2 cyear*4,cDDD*3,cmonth,cday
        integer mday(12),leap
        integer remn
        iyear = 2007
        write(cyear,"(I4)") iyear
        print*,"INPUT THE year and DDD: year  DDD"
        read(*,*)iyear,iDDD
      print*,"Processing ",iyear,iDDD
        write(cDDD,"(I3)") iDDD
      if ( iDDD .LT. 10 )  cDDD(1:2)="00"
      if ( iDDD .LT. 100 )  cDDD(1:1)="0"

      iday =  iDDD
        data mday/31,28,31,30,31,30,31,31,30,31,30,31/
        imonth = 01
        leap = 0
        remn = 0
        remn = mod(iyear,100)

        if (remn == 0 ) then
                          remn =mod(iyear,400)

                if(remn == 0)then
                leap = 1
                endif
        else
                 remn =mod(iyear,4)

                        if(remn == 0)then
                leap = 1
                else
                leap = 0
                endif

        endif
        if ( leap == 1)then
         mday(2)=mday(2)+1
        endif

      do while (iday > mday(imonth))
         iday = iday - mday(imonth)
       
         imonth = imonth+1
      enddo
        print*,iyear,imonth, iday
        write(cmonth,"(I2)") imonth
        write(cday,"(I2)") iday

        if ( imonth .LT. 10 )  cmonth(1:1)="0"
        if ( iday   .LT. 10 )  cday  (1:1)="0"
        write(*,"(A6,A5,2A4)")cyear,cDDD,cmonth,cday
end

topmad

Lauretta 发表于 2012-10-25 14:37

                               
登录/注册后可看大图

一年的第几天啊。。。不过SA雷达数据说明里的那个儒略日是从1970年1月1日开始，显示出来得写成公历才行。 ...

你这个是计算某一天是一年中的第几天的程序吗？
论坛有个很简单的算法计算day of year的

Aires

不大懂这数据，感谢楼主分享自己的程序~

rosewei

LZ，LZ，请教请教，如何将基数据的坐标系转成直接坐标系？

Lauretta

rosewei 发表于 2012-10-26 17:12

                               
登录/注册后可看大图

LZ，LZ，请教请教，如何将基数据的坐标系转成直接坐标系？

转笛卡尔坐标？我也在头疼这事呢，matlab里有自带函数可以转，但是我那样写的就不是转到固定格点上。[px1,py1,pz1]=sph2cart(X,Y,Z);读基数据，然后这样转的。

rosewei

Lauretta 发表于 2012-11-7 14:17

                               
登录/注册后可看大图

转笛卡尔坐标？我也在头疼这事呢，matlab里有自带函数可以转，但是我那样写的就不是转到固定格点上。=sph ...

就是转成笛卡尔坐标，雷达数据编程的第一步。我也是个没人带的孩子，要多走很多弯路地。。。。

康理想

请问楼主解决好雷达基数据读取和转换到固定格点的问题没，急需！！