NOAA的卫星反演欧亚大陆雪盖资料如何处理？

水草

我做毕业论文用的是美国冰雪资料中心提供的1966年10月-2007年5月北半球25km等面积NSIDC卫星反演雪盖资料，在每个网格上都包括了逐周的雪盖信息：1表示有陆地雪覆盖，0表示无陆地雪覆盖等，下面是我读其中一周资料的程序
implicit none
integer i,j,n,k,m
parameter(n=721)
character*3 swe(n,n)
open(1,file='E:\data\NL19661017-19661023.v03.SI',form='binary')
do i=1,n
do j=1,n
read(1,end=10)swe(i,j)
print *,(ichar(swe(i,j)))
enddo
enddo
10 close(1)
open(2,file='E:\NL19661017-19661023.txt')
do i=1,n
write(2,*)(ichar(swe(i,j)),j=1,n)
enddo
close(2)
print*,i-1,j-1
end
读出的数据中只有255，254,0,1,5这几个数字，我不知道255,254，和5代表什么？或者是我的程序出问题了？求用过这个资料的人帮帮忙，否则我的毕业论文无法进行下去啊！

mofangbao

没用过这个资料，不过你应该仔细看数据说明，至少明白数据是以什么格式存放的

水草

The Northern Hemisphere EASE-Grid Weekly Snow Cover and Sea Ice Extent Version 3 product and Version 3.1 update combine snow cover and sea ice extent at weekly intervals from 23 October 1978 through 24 June 2007, and snow cover alone from 3 October 1966 through 23 October 1978. The data set is the first representation of combined snow and sea ice measurements derived from satellite observations for the period of record. Snow cover extent is based on the digital NOAA-NESDIS Weekly Northern Hemisphere Snow Charts, revised by D. Robinson (Rutgers University) and regridded to the EASE-Grid. The original NOAA-NESDIS weekly snow charts are derived from the manual interpretation of AVHRR, GOES, and other visible-band satellite data. Sea ice extent is regridded to EASE-Grid from Sea Ice Concentrations from Nimbus-7 SMMR and DMSP SSM/I Passive Microwave Data. Designed to facilitate study of Northern Hemisphere seasonal fluctuations of snow cover and sea ice extent, the data set also includes monthly climatologies describing average extent, probability of occurrence, and variance. Data are provided in the Northern Hemisphere 25 km equal-area grid (NSIDC Nl EASE-Grid). Data are provided in flat, unsigned binary files and are available via FTP. 这是数据描述 我没在描述中发现每个格点到底代表什么

map999

一般用 grib或nc格式存放文件，请告诉下载数据的地址，我下载几个数据进行解析后估计可以解决你关心的问题

水草

http://nsidc.org/data/nsidc-0046.html

urbanwang

是啊，我也没弄明白这个数据是怎么处理的，伤心ing

urbanwang1986

请问楼主，该资料如何进行投影转换？谢谢

水草

urbanwang1986 发表于 2012-5-11 16:51

                               
登录/注册后可看大图

请问楼主，该资料如何进行投影转换？谢谢

这是我从网上下的坐标互相转换的函数子程序，我也不太懂，你看一下吧
C==========================================================================
C ezlhconv.for - FORTRAN routines for conversion of azimuthal
C  equal area and equal area cylindrical grid coordinates
C
C 30-Jan.-1992 H.Maybee
C 20-Mar-1992 Ken Knowles  303-492-0644  knowles@kryos.colorado.edu
C       16-Dec-1993 MJ Brodzik   303-492-8263  brodzik@jokull.colorado.edu
C                   Copied from nsmconv.f, changed resolutions from
C                   40-20-10 km to 25-12.5 km
C       21-Dec-1993 MJ Brodzik   303-492-8263  brodzik@jokull.colorado.edu
C                   Fixed sign of Southern latitudes in ease_inverse.
C 12-Sep-1994 David Hoogstrate 303-492-4116 hoogstra@jokull.colorado.edu
C      Changed grid cell size. Changed "c","f" to "l","h"
C 25-Oct-1994 David Hoogstrate 303-492-4116 hoogstra@jokull.colorado.edu
C      Changed row size from 587 to 586 for Mercator projection
C      Changed function names to "ezlh-.."
C$Log: ezlhconv.f,v $
CRevision 1.3  1994/11/01 23:40:43  brodzik
CReplaced all references to 'ease' with 'ezlh'
C
C==========================================================================
C--------------------------------------------------------------------------
function ezlh_convert (grid, lat, lon, r, s)
character*(*) grid
real lat, lon, r, s
C
C convert geographic coordinates (spherical earth) to
C azimuthal equal area or equal area cylindrical grid coordinates
C
C status = ezlh_convert (grid, lat, lon, r, s)
C
C input : grid - projection name '[NSM][lh]'
C               where l = "low"  = 25km resolution
C                     h = "high" = 12.5km resolution
C  lat, lon - geo. coords. (decimal degrees)
C
C output: r, s - column, row coordinates
C
C result: status = 0 indicates normal successful completion
C   -1 indicates error status (point not on grid)
C
C--------------------------------------------------------------------------
        integer cols, rows, scale
real Rg, phi, lam, rho
C radius of the earth (km), authalic sphere based on International datum
parameter (RE_km = 6371.228)
C nominal cell size in kilometers
parameter (CELL_km = 25.067525)
C scale factor for standard paralles at +/-30.00 degrees
parameter (COS_PHI1 = .866025403)
parameter (PI = 3.141592653589793)
rad(t) = t*PI/180.
deg(t) = t*180./PI
ezlh_convert = -1
if ((grid(1:1).eq.'N').or.(grid(1:1).eq.'S')) then
   cols = 721
   rows = 721
else if (grid(1:1).eq.'M') then
   cols = 1383
   rows = 586
else
   print *, 'ezlh_convert: unknown projection: ', grid
   return
endif
if (grid(2:2).eq.'l') then
   scale = 1
else if (grid(2:2).eq.'h') then
   scale = 2
else
   print *, 'ezlh_convert: unknown projection: ', grid
   return
endif
        Rg = scale * RE_km/CELL_km
C
C r0,s0 are defined such that cells at all scales
C have coincident center points
C
        r0 = (cols-1)/2. * scale
        s0 = (rows-1)/2. * scale
phi = rad(lat)
        lam = rad(lon)
if (grid(1:1).eq.'N') then
   rho = 2 * Rg * sin(PI/4. - phi/2.)
   r = r0 + rho * sin(lam)
   s = s0 + rho * cos(lam)
else if (grid(1:1).eq.'S') then
   rho = 2 * Rg * cos(PI/4. - phi/2.)
   r = r0 + rho * sin(lam)
   s = s0 - rho * cos(lam)
        else if (grid(1:1).eq.'M') then
          r = r0 + Rg * lam * COS_PHI1
          s = s0 - Rg * sin(phi) / COS_PHI1
endif
ezlh_convert = 0
return
end
C--------------------------------------------------------------------------
function ezlh_inverse (grid, r, s, lat, lon)
character*(*) grid
real r, s, lat, lon
C
C convert azimuthal equal area or equal area cylindrical
C grid coordinates to geographic coordinates (spherical earth)
C
C status = ezlh_inverse (grid, r, s, lat, lon)
C
C input : grid - projection name '[NSM][lh]'
C               where l = "low"  = 25km resolution
C                     h = "high" = 12.5km resolution
C  r, s - grid column and row coordinates
C
C output: lat, lon - geo. coords. (decimal degrees)
C
C result: status = 0 indicates normal successful completion
C   -1 indicates error status (point not on grid)
C
C--------------------------------------------------------------------------
        integer cols, rows, scale
real Rg, phi, lam, rho
real gamma, beta, epsilon, x, y
real sinphi1, cosphi1
C radius of the earth (km), authalic sphere based on International datum
parameter (RE_km = 6371.228)
C nominal cell size in kilometers
parameter (CELL_km = 25.067525)
C       scale factor for standard paralles at +/-30.00 degrees
        parameter (COS_PHI1 = .866025403)
        parameter (PI = 3.141592653589793)
        rad(t) = t*PI/180.
        deg(t) = t*180./PI
ezlh_inverse = -1
if ((grid(1:1).eq.'N').or.(grid(1:1).eq.'S')) then
   cols = 721
   rows = 721
else if (grid(1:1).eq.'M') then
   cols = 1383
   rows = 586
else
   print *, 'ezlh_inverse: unknown projection: ', grid
   return
endif
if (grid(2:2).eq.'l') then
   scale = 1
else if (grid(2:2).eq.'h') then
   scale = 2
else
   print *, 'ezlh_inverse: unknown projection: ', grid
   return
endif
        Rg = scale * RE_km/CELL_km
        r0 = (cols-1)/2. * scale
        s0 = (rows-1)/2. * scale
        x = r - r0
        y = -(s - s0)
        if ((grid(1:1).eq.'N').or.(grid(1:1).eq.'S')) then
          rho = sqrt(x*x + y*y)
          if (rho.eq.0.0) then
            if (grid(1:1).eq.'N') lat = 90.0
            if (grid(1:1).eq.'S') lat = -90.0
            lon = 0.0
   else
            if (grid(1:1).eq.'N') then
              sinphi1 = sin(PI/2.)
              cosphi1 = cos(PI/2.)
              if (y.eq.0.) then
                if (r.le.r0) lam = -PI/2.
                if (r.gt.r0) lam = PI/2.
              else
                lam = atan2(x,-y)
       endif
            else if (grid(1:1).eq.'S') then
              sinphi1 = sin(-PI/2.)
              cosphi1 = cos(-PI/2.)
              if (y.eq.0.) then
                if (r.le.r0) lam = -PI/2.
                if (r.gt.r0) lam = PI/2.
              else
                lam = atan2(x,y)
       endif
          endif
            gamma = rho/(2 * Rg)
     if (abs(gamma).gt.1.) return
            c = 2 * asin(gamma)
            beta = cos(c) * sinphi1 + y * sin(c) * (cosphi1/rho)
     if (abs(beta).gt.1.) return
            phi = asin(beta)
            lat = deg(phi)
            lon = deg(lam)
          endif
else if (grid(1:1).eq.'M') then
C
C   allow .5 cell tolerance in arcsin function
C   so that grid coordinates which are less than .5 cells
C   above 90.00N or below 90.00S are given a lat of 90.00
C
   epsilon = 1 + 0.5/Rg
          beta = y*COS_PHI1/Rg
          if (abs(beta).gt.epsilon) return
   if (beta.le.-1.) then
     phi = -PI/2.
   else if (beta.ge.1.) then
     phi = PI/2.
          else
     phi = asin(beta)
   endif
          lam = x/COS_PHI1/Rg
          lat = deg(phi)
          lon = deg(lam)
endif
ezlh_inverse = 0
return
end

urbanwang1986

水草 发表于 2012-5-11 20:32

                               
登录/注册后可看大图

这是我从网上下的坐标互相转换的函数子程序，我也不太懂，你看一下吧
C=============================== ...

好的。赞。我研究研究，如搞明白，及时给您回复