在 mac OSX 上使用 gfortran 读取免费提供的 UVR 数据
Reading freely available UVR data using gfortran on mac OSX
我想用fortran 读取日本宇宙航空研究开发机构制作的紫外线辐射数据。该数据采用 2000-2010 年的每日和每月时间分辨率,空间分辨率约为 5 公里。这个问题值得回答,因为数据可能对许多 environment/health 项目有用,并且可以免费获得,并适当确认来源并共享任何后续出版物的预印本,来自:
ftp://suzaku.eorc.jaxa.jp/pub/GLI/glical/Global_05km/monthly/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:PAR 范围表面反射率(canopy/solid 表面的顶部)= DN * 0.0001(仅月度数据)
错误值
-1 作为有符号短整数 (int16)
65535 作为无符号短整数 (uint16)
目前进度
我已经在 mac OSX 上成功下载并安装了 gfortran。我已经下载了一个测试文件(MOD02SSH_A20000224Av6_v601_7200_3601_uvb__le.gz)并解压了。我已经创建了一个程序文件:
PROGRAM readuvr
IMPLICIT NONE
!some code
END PROGRAM
然后我将在命令行中键入以下内容以创建一个可执行文件并运行它来提取数据。
gfortran -o executable
./executable
作为 Fortran 的完全初学者,我的问题是:如何使用提供的说明构建一个可以读取数据并将其输出到文本文件的程序?
好吧,该文件扩展到 51,868,800 字节。评论暗示 header 是 14,400 字节,剩下 51,854,400 字节的实际数据负载。
好像有7200行数据,也就是说每行有7202字节。似乎有 2 个字节(16 位样本)所以如果我们假设 2 bytes/sample,这意味着每行有 3601 个样本,这与 ml=3601
.
匹配
所以基本上,您需要读取 14,400 个字节的 header,然后是 7200 行数据,每行包含 3601 个值,每个值都是 2 个字节宽...
实际上,如果您对 FORTRAN 不太熟悉,您可能会喜欢使用 Perl 来提取数据,Perl 已经安装并且可以在 OS X 上使用。我已经启动了一个非常简单的 Perl 程序,它读取数据并在每行上打印前 2 个值:
#!/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
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
这实际上没有经过验证,因为我没有已知的文件内容可以与之比较。
我想用fortran 读取日本宇宙航空研究开发机构制作的紫外线辐射数据。该数据采用 2000-2010 年的每日和每月时间分辨率,空间分辨率约为 5 公里。这个问题值得回答,因为数据可能对许多 environment/health 项目有用,并且可以免费获得,并适当确认来源并共享任何后续出版物的预印本,来自:
ftp://suzaku.eorc.jaxa.jp/pub/GLI/glical/Global_05km/monthly/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:PAR 范围表面反射率(canopy/solid 表面的顶部)= DN * 0.0001(仅月度数据)
错误值
-1 作为有符号短整数 (int16)
65535 作为无符号短整数 (uint16)
目前进度
我已经在 mac OSX 上成功下载并安装了 gfortran。我已经下载了一个测试文件(MOD02SSH_A20000224Av6_v601_7200_3601_uvb__le.gz)并解压了。我已经创建了一个程序文件:
PROGRAM readuvr
IMPLICIT NONE
!some code
END PROGRAM
然后我将在命令行中键入以下内容以创建一个可执行文件并运行它来提取数据。
gfortran -o executable
./executable
作为 Fortran 的完全初学者,我的问题是:如何使用提供的说明构建一个可以读取数据并将其输出到文本文件的程序?
好吧,该文件扩展到 51,868,800 字节。评论暗示 header 是 14,400 字节,剩下 51,854,400 字节的实际数据负载。
好像有7200行数据,也就是说每行有7202字节。似乎有 2 个字节(16 位样本)所以如果我们假设 2 bytes/sample,这意味着每行有 3601 个样本,这与 ml=3601
.
所以基本上,您需要读取 14,400 个字节的 header,然后是 7200 行数据,每行包含 3601 个值,每个值都是 2 个字节宽...
实际上,如果您对 FORTRAN 不太熟悉,您可能会喜欢使用 Perl 来提取数据,Perl 已经安装并且可以在 OS X 上使用。我已经启动了一个非常简单的 Perl 程序,它读取数据并在每行上打印前 2 个值:
#!/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
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
这实际上没有经过验证,因为我没有已知的文件内容可以与之比较。