风云四号行列号与经纬度转换

一阵妖风

长恨晴天夕照晚 发表于 2019-3-18 10:19
您好，能问一下raw格式的文件怎么读取吗？

您好，请问raw文件的读取问题解决了吗？有几个问题想问一下？

一阵妖风

laoping 发表于 2019-8-12 10:46
我也想请问老师raw格式的文件怎么打开？

您好，请问raw文件的读取问题解决了吗？有几个问题想问一下

李山山0103

你好，楼主这个问题解决了吗？

laoping

open(10,file="FullMask_Grid_4000.raw",form="unformatted",access="direct",recl=Nc*Nl*16)
        read(10,rec=1) ((lat(i,j),lon(i,j),i=1,Nc),j=1,Nl)
        close(10)
Nl=Nc=2748

hpzheng

请问楼主这个问题解决了吗？

laoping

hpzheng 发表于 2020-3-15 19:52
请问楼主这个问题解决了吗？

可以直接用fortran读取，就当作无格式二进制文件来读，不过需要你知道他的数组大小及字符长度。
CC Nc: Number of colums
CC Nl: Number of lines
        integer*4 Nc,Nl
        parameter(Nc=2748,Nl=2748)
        real*8 lon(Nc,Nl),lat(Nc,Nl)

        open(10,file="./FullMask_Grid_4000.raw",form="unformatted",
     &                  access="direct",recl=Nc*Nl*16)
        read(10,rec=1) ((lat(i,j),lon(i,j),i=1,Nc),j=1,Nl)
        close(10)

五蕴

fileName = "E:\科研\卫星\FullMask_Grid_4000\FullMask_Grid_4000.raw";
fileID=fopen(fileName);
BB=fread(fileID,[2748*2 2748],'double');
fclose(fileID);
LAT=BB(1:2:2*2748,:);
LON=BB(2:2:2*2748,:);

五蕴

#将FullMask_Grid_4000.raw读取为NC格式的经纬度文件，经纬度维度为2748*2748，以下为MATLAB语言
file ID=fopen('FullMask_Grid_4000.raw');
BB=fread(fileID,[2748*2 2748],'double');
fclose(fileID);
LAT=BB(1:2:2*2748,:);
LON=BB(2:2:2*2748,:);
LAT(find(abs(LAT)>1000))=NaN;
LON(find(abs(LON)>1000))=NaN;
LAT=LAT';
LON=LON';

五蕴

用公式将行列号转换成经纬度，行列号（默认是int16）需要变成double类型数据，这样就不会出现复数了。ColumnNumber = double(ColumnNumber);
LineNumber = double(LineNumber);
LineNumber(find(LineNumber==-1))=NaN;
ColumnNumber(find(ColumnNumber==-1))=NaN;
记得把无效值-1改成NaN

李嗣源

我用的是Python,尝试了一下，好像是成功了，只需要把两个数组先网格化之后再去按公式计算（但是好像算出来是反的，然后我把矩阵转置了一下），把计算结果和用raw文件对照表得出来的结果对比画出来的图是一样的。我贴出来

l=np.arange(0,2748)
c=np.arange(0,2748)
L,C=np.meshgrid(l,c)#网格化
coff=1373.5
cfac=10233137
loff=1373.5
lfac=10233137
h=42164
ea=6378.137
eb=6356.7523
ld=104.7

x=(np.pi*(C-coff))/(180*2**-16*cfac)
y=(np.pi*(L-loff))/(180*2**-16*lfac)
sd=((h*np.cos(x)*np.cos(y))**2-(np.cos(y)**2+(ea**2/eb**2)*np.sin(y)**2)*(h**2-ea**2))**0.5
sn=(h*np.cos(x)*np.cos(y)-sd)/(np.cos(y)**2+(ea**2/eb**2)*np.sin(y)**2)
s1=h-sn*np.cos(x)*np.cos(y)
s2=sn*np.sin(x)*np.cos(y)
s3=-1*sn*np.sin(y)
sxy=(s1**2+s2**2)**0.5
lon0=(180/np.pi)*np.arctan(s2/s1)+ld
lat0=(180/np.pi)*np.arctan((ea**2/eb**2)*(s3/sxy))
lat0=np.transpose(lat0)
lon0=np.transpose(lon0)