坐标转换

水草

想问一下我这个程序中输入相同纬度，不同经度时，为什么转换的结果 r 不相同呢？是程序写的有错误吗？

program cal_r

implicit none

character       :: grid

real                 ::    lat,  lon,  r,s

real,external    ::    ezlh_convert

write(*,*) ezlh_convert('Nl', 80., 10., r, s)

      write(*,*) r,s

end

     

function ezlh_convert(grid, lat, lon, r, s)

character*(*) grid

real lat, lon, r, s

!

! convert geographic coordinates (spherical earth) to

! azimuthal equal area or equal area cylindrical grid coordinates

!

! status = ezlh_convert (grid, lat, lon, r, s)

!

! input : grid - projection name '[NSM][lh]'

!               where l = "low"  = 25km resolution

!                     h = "high" = 12.5km resolution

!  lat, lon - geo. coords. (decimal degrees)

!

! output: r, s - column, row coordinates

!

! result: status = 0 indicates normal successful completion

!   -1 indicates error status (point not on grid)

!

!--------------------------------------------------------------------------

      integer cols, rows, scale

real Rg, phi, lam, rho

! radius of the earth (km), authalic sphere based on International datum

parameter (RE_km = 6371.228)

! nominal cell size in kilometers

parameter (CELL_km = 25.067525)

! 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

!

! r0,s0 are defined such that cells at all scales

! have coincident center points

!

        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

</font></div><div align=\"left\"><font style=\"font-size:12px\"><br /> </font></div><div align=\"left\">program cal_r</div><br /> <div align=\"left\">implicit none<br /> character       :: grid<br /> real                 ::    lat,  lon,  r,s<br /> real,external    ::    ezlh_convert</div><div align=\"left\"><br /> write(*,*) ezlh_convert('Nl', 80., 10., r, s)<br />       write(*,*) r,s<br /> end<br />      <br /> function ezlh_convert(grid, lat, lon, r, s)<br /> character*(*) grid<br /> real lat, lon, r, s<br /> !<br /> ! convert geographic coordinates (spherical earth) to <br /> ! azimuthal equal area or equal area cylindrical grid coordinates<br /> !<br /> ! status = ezlh_convert (grid, lat, lon, r, s)<br /> !<br /> ! input : grid - projection name '[NSM][lh]'<br /> !               where l = "low"  = 25km resolution<br /> !                     h = "high" = 12.5km resolution<br /> !  lat, lon - geo. coords. (decimal degrees)<br /> !<br /> ! output: r, s - column, row coordinates<br /> !<br /> ! result: status = 0 indicates normal successful completion<br /> !   -1 indicates error status (point not on grid)<br /> !<br /> !--------------------------------------------------------------------------<br />       integer cols, rows, scale<br /> real Rg, phi, lam, rho</div><div align=\"left\">! radius of the earth (km), authalic sphere based on International datum <br /> parameter (RE_km = 6371.228)<br /> ! nominal cell size in kilometers<br /> parameter (CELL_km = 25.067525)</div><div align=\"left\">! scale factor for standard paralles at +/-30.00 degrees<br /> parameter (COS_PHI1 = .866025403)</div><div align=\"left\">parameter (PI = 3.141592653589793)<br /> rad(t) = t*PI/180.<br /> deg(t) = t*180./PI</div><div align=\"left\">ezlh_convert = -1</div><div align=\"left\">if ((grid(1:1).eq.'N').or.(grid(1:1).eq.'S')) then<br />    cols = 721<br />    rows = 721<br /> else if (grid(1:1).eq.'M') then<br />    cols = 1383<br />    rows = 586<br /> else<br />    print *, 'ezlh_convert: unknown projection: ', grid<br />    return<br /> endif</div><div align=\"left\">if (grid(2:2).eq.'l') then<br />    scale = 1<br /> else if (grid(2:2).eq.'h') then<br />    scale = 2<br /> else<br />    print *, 'ezlh_convert: unknown projection: ', grid<br />    return<br /> endif</div><div align=\"left\">        Rg = scale * RE_km/CELL_km</div><div align=\"left\">!<br /> ! r0,s0 are defined such that cells at all scales<br /> ! have coincident center points<br /> !<br />         r0 = (cols-1)/2. * scale<br />         s0 = (rows-1)/2. * scale</div><div align=\"left\">phi = rad(lat)<br />         lam = rad(lon)</div><div align=\"left\">if (grid(1:1).eq.'N') then<br />    rho = 2 * Rg * sin(PI/4. - phi/2.)<br />    r = r0 + rho * sin(lam)<br />    s = s0 + rho * cos(lam)</div><div align=\"left\">else if (grid(1:1).eq.'S') then<br />    rho = 2 * Rg * cos(PI/4. - phi/2.)<br />    r = r0 + rho * sin(lam)<br />    s = s0 - rho * cos(lam)</div><div align=\"left\">        else if (grid(1:1).eq.'M') then<br />           r = r0 + Rg * lam * COS_PHI1<br />           s = s0 - Rg * sin(phi) / COS_PHI1</div><div align=\"left\">endif</div><div align=\"left\">ezlh_convert = 0<br /> return<br /> end</div><font style=\"font-size:12px\">

纬度40，经度10的时候结果是397.30444 ，571.5638

纬度80，经度10的时候结果是367.6932，403.6304

纬度90，经度10的时候结果是360.000,360.000

artmunich

提醒一下，论坛支持语法加亮等功能，贴代码的时候能用一下更便于别人阅读~

水草

artmunich 发表于 2012-4-29 20:12

                               
登录/注册后可看大图

提醒一下，论坛支持语法加亮等功能，贴代码的时候能用一下更便于别人阅读~

哦~谢谢！我暂时还不太了解论坛的功能，以后注意

水草

问题解决了啊~~

topmad

水草 发表于 2012-4-29 22:23

                               
登录/注册后可看大图

问题解决了啊~~

怎么解决的呢，发上来大家一起学习一下呗

godenflame135

水草 发表于 2012-4-29 20:30

                               
登录/注册后可看大图

哦~谢谢！我暂时还不太了解论坛的功能，以后注意

已高显！

godenflame135

水草 发表于 2012-4-29 22:23

                               
登录/注册后可看大图

问题解决了啊~~

怎么解决的，说说？

水草

topmad 发表于 2012-4-29 22:27

                               
登录/注册后可看大图

怎么解决的呢，发上来大家一起学习一下呗

我上面那个程序是想把地理坐标转化成格点坐标，这样转化的结果是格点会出现小数，取舍很麻烦，无法准确的找出需要的格点，所以我那个思路就是错的！正确的思路应该是把格点坐标先转化成地理坐标，然后把在需要的地理区域的格点挑出来，我说的是不是很混乱？

水草

godenflame135 发表于 2012-4-29 23:58

                               
登录/注册后可看大图

怎么解决的，说说？

这个程序本身是没有错的，但我这个思路是错的。应该是先把格点坐标转化成地理坐标，再选择需要的格点数据

水草

这是北半球EASE-grid的格点坐标转化成地理坐标的程序
[code] implicit none
      character grid
        integer a
        real r,s,lat,lon
        real,external::ezlh_inverse
      open(2,file='E:\geo.txt')
      open(3,file='E:\selectedgrid.txt')
        do r=1,721
          do s=1,721
          a=ezlh_inverse ('Nl', r, s, lat, lon)
          write(2,*) lat,lon
          if(lat>=40.0)then
            if(lon>=10.0.and.lon<=140.0)then
          write(3,*) r,s
            end if
          end if
          end do
        end do
        close(2)
        close(3)
        end
        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[code]