在mac上使用gfortran读取可免费获取的UVR数据

Question

我想用fortran阅读日本宇航研究开发机构提供的紫外线辐射数据。该数据的日和月时间分辨率为2000-2010年，空间分辨率约为5公里。这一问题值得回答，因为这些数据可能对若干环境/卫生项目有用，并可免费获得，并适当确认资料来源，并分享今后任何出版物的预印本，网址如下：

05公里/月/uvb/

有一个自述文件可用，它提供关于如何使用fortran读取数据的说明，如下所示：

_le文件说明

标题

Read header (size= pixel size *2byte):
character head*14400

read(10,rec=1) head
read(head,'(2i6,2f8.2,f8.4,2e12.5,a1,a8,a1,a40)') 
     & npixel,nline,lon_min,lat_max,reso,slope,offset,',',
     & para,',',outfile

读取数据(例如，fortran77)

parameter(nl=7200, ml=3601)

... open file by "unformatted", "recl=nl*2(byte)" (,"bytereclen")

integer*2 i2buf(nl,ml)
do m=1,ml
 read(10,rec=1+m) (i2buf(n,m), n=1,nl)
 do n=1,nl
   par=i2buf(n,m)*slope+offset
   write(6,*) 'PAR[Ein/m^2/day]=',par
 enddo
enddo

斜率值

par__le :日PAR Ein/m^2/日= DN * 0.01

dpar_le :直接PAR = DN * 0.01

swr__le :日平均短波辐射W/m^2 = DN * 0.01

tip__le :中午瞬时PAR的透过率= DN * 0.0001

uva__le :日平均值UVA W/m^2 = DN * 0.001

uvb__le :每日平均UVB W/m^2 = DN * 0.0001

rpar_le :标准距离表面反射率(冠层/实心面顶部)= DN * 0.0001 (仅每月数据)

误差值

-1为有符号短整数(int16)

65535为无符号短整数(uint16)

进展到目前为止

我已经在mac上下载并成功安装了gfortran。我下载了一个测试文件(MOD02SSH_A20000224Av6_v601_7200_3601_uvb__le.gz)并对其进行了解压缩。我创建了一个程序文件：

PROGRAM readuvr
IMPLICIT NONE

!some code

END PROGRAM

然后，我将在命令行中输入以下内容来创建一个可执行文件，并运行它来提取数据。

gfortran -o executable

./executable

作为fortran的完全初学者，我的问题是:如何使用所提供的指令来构建一个程序来读取数据并将其输出到文本文件中？

Answer 1

这个文件扩展到51,868,800字节。注释意味着报头为14,400字节，这就留下了51,854,400字节的实际数据有效负载。

似乎有7200行数据，所以这意味着每行有7202字节。这里似乎有2个字节(16位样本)，所以如果我们假设为2个字节/样本，这意味着每行有3601个样本，这与ml=3601相匹配。

所以基本上，您需要读取14,400字节的标题，然后是7200行数据，每一行由3601个值组成，每个值都是2字节宽.

实际上，如果您对FORTRAN非常陌生，您可能希望使用Perl提取数据，因为Perl已经在OS X上安装并可用。我已经启动了一个非常简单的Perl程序，它读取dat并在每一行上打印前两个值：

#!/usr/bin/perl
use strict;
use warnings;

# Read 14,400 bytes of header
my $buffer;
my $nBytes = 14400; 
my $bytesRead = read (STDIN, $buffer, $nBytes) ;
my ($npixel,$nline,$lon_min,$lat_max,$reso,$slope,$offset,$junk)=split(' ',$buffer);
print "npixel:$npixel\n";
print "nline:$nline\n";
print "lon_min:$lon_min\n";
print "lat_max:$lat_max\n";
print "reso:$reso\n";
print "slope:$slope\n";
$offset =~ s/,.*//; # strip trailing comma and junk
print "offset:$offset\n";

# Read actual lines of data
my $line;
for(my $m=1;$m<=$nline;$m++){
   read(STDIN,$line,$npixel*2);
   my $x=$npixel*2;
   my @values=unpack("S$x",$line);
   printf "Line: %d",$m;
   for(my $j=0;$j<2;$j++){
      printf ",%f",$values[$j]*$slope+$offset;
   }
   printf "\n"; # newline
}

将其保存为go.pl，然后在终端中键入以下一次以使其可执行

chmod +x go.pl

然后像这样运行

./go.pl < MOD02SSH_A20000224Av6_v601_7200_3601_uvb__le

样本输出提取：

npixel:7200
nline:3601
lon_min:0.00
lat_max:90.00
reso:0.0500
slope:0.10000E-03
offset:0.00000E+00
...
...
Line: 3306,0.099800,0.099800
Line: 3307,0.099900,0.099900
Line: 3308,0.099400,0.074200
Line: 3309,0.098900,0.098900
Line: 3310,0.098400,0.098400
Line: 3311,0.074300,0.074200
Line: 3312,0.071300,0.071200

Answer 2

fortran (f2003左右)解决方案。(顺便说一句，链接说明很糟糕)

      implicit none
      character*80 para,outfile
      character(len=:),allocatable::header,infile      
      integer npixel,nline,blen,i
c note kind=2 is not standard. This needs to be a 2-byte integer.
      integer(kind=2),allocatable :: data(:,:)
      real lon_min,lat_max,reso,slope,off
c header is plain text, so first open formatted and
c directly read header data
      infile='MOD02SSH_A20000224Av6_v601_7200_3601_uvb__le'
      open(10,file=infile)
      read(10,*)npixel,nline,lon_min,lat_max,reso,slope,off,
     $     para,outfile
      close(10)
      write(*,*)npixel,nline,lon_min,lat_max,reso,slope,off,
     $     trim(para),' ',trim(outfile)
      blen=2*npixel
      allocate(character(len=blen)::header)
      allocate(data(npixel,nline))
      if( sizeof(data(1,1)).ne.2 )then
         write(*,*)'error kind=2 did not give a 2 byte integer'
         stop
      endif
c now close and reopen for binary read.
c direct access approach:
      open(20,file=infile,access='direct',recl=blen/4)      
c note the granularity of the recl= specifier is not standard.
c ifort uses 4 bytes.  (note this will break if npixel is not even )
      read(20,rec=1)header
      write(*,*)trim(header)
      do i=1,nline
         read(20,rec=i+1)data(:,i)
      enddo
c note streams if available is simpler: (we don't need to know rec len )
c      open(20,file=infile,access='stream')
c      read(20)header,data
      end

这实际上是没有验证的，因为我没有可与之比较的文件内容。