我如何转换纬度-经度以缓解美国宇航局的网格？

Question

我正在处理从美国宇航局eosdis下载的he5文件。

我使用r中的rhdf5包成功地读取了文件。

文件具有子数据集，由一个矩阵，dim 721x721组成。

据我所知，该文件是作为easegrid构建的，它没有任何关于坐标的信息，因此我无法找到在哪个网格(矩阵中的元素)中我的工作场所。

是否有任何方法来转换经度值以简化网格？

非常感谢

Answer 1

我发现了一个python软件，它能够在杨百翰大学FTP服务器上完成不同地图系统之间的几次转换。

您可能希望在此文件中使用ease2grid函数。

我还修复了文件上的几个衬里问题。这是文件；

"""
Created on Sep 21, 2011
Revised on Apr 7, 2017 + include EASE2 support

@author: Bradley, DGL
"""
from __future__ import division
from numpy import cos, sin, tan, mod, sqrt, all
import numpy as np
from .ease2helper import ease2_map_info, easeconv_normalize_degrees


def latlon2pix(alon, alat, head):
    """Latitude/longitude to pixels
    (x, y) = latlon2pix(lon,lat,head)

    Convert a lat,lon coordinate (lon,lat) to an image pixel location
    (x,y) (in floating point, matlab convention).
    To compute integer pixel indices (ix,iy): check to insure
    1 <= x < nsx+1 and 1 <= x < nsx+1 then ix=floor(x) iy=floor(y)

    INPUTS:
        lon,lat - longitude, latitude
        head - header array from load sir

    OUTPUTS:
        x,y - pixel location (matlab coordinates y_matlab=nxy-y_sir+1)
    """
    nsx = head[0]  # noqa: F841
    nsy = head[1]
    iopt = head[16]
    xdeg = head[2]
    ydeg = head[3]
    ascale = head[5]
    bscale = head[6]
    a0 = head[7]
    b0 = head[8]

    if iopt == -1:  # image only (can't transform!)
        x = ascale * (alon - a0)
        y = bscale * (alat - b0)
    elif iopt == 0:  # rectalinear lat/lon
        thelon = alon
        thelat = alat
        x = ascale * (thelon - a0)
        y = bscale * (thelat - b0)
    elif (iopt == 1) or (iopt == 2):  # lambert
        thelon, thelat = lambert1(alat, alon, ydeg, xdeg, iopt)
        x = ascale * (thelon - a0)
        y = bscale * (thelat - b0)
    elif iopt == 5:  # polar stereographic
        thelon, thelat = polster(alon, alat, xdeg, ydeg)
        x = (thelon - a0) / ascale
        y = (thelat - b0) / bscale
    elif (iopt == 8) or (iopt == 9) or (iopt == 10):  # EASE2
        thelon, thelat = ease2grid(iopt, alat, alon, ascale, bscale)
        x = thelon + 1.0 - a0
        y = thelat + 1.0 + b0
    elif (iopt == 11) or (iopt == 12) or (iopt == 13):  # EASE
        thelon, thelat = easegrid(iopt, alat, alon, ascale)
        thelon = thelon + xdeg
        thelat = thelat + ydeg
        x = thelon - (xdeg + a0)
        y = thelat - (ydeg + b0)
    else:
        print("*** Unknown SIR transformation: %d" % iopt)

    y = nsy - y - 1.0  # convert from matlab coordinates to SIR coordinates
    return x, y


def lambert1(lat, lon, orglat, orglon, iopt):
    """Lambert azimuthal equal-area projection
    (x,y)=lambert1(lat,lon,orglat,orglon,iopt)

    Computes the transformation from lat/lon to x/y for the
    lambert azimuthal equal-area projection

    inputs:
        lat    (r): latitude +90 to -90 deg with north positive
        lon    (r): longitude 0 to +360 deg with east positive
            or -180 to +180 with east more positive
        orglat    (r): origin parallel +90 to -90 deg with north positive
        orglon    (r): central meridian (longitude) 0 to +360 deg
            or -180 to +180 with east more positive
        iopt    (i): earth radius option
            for iopt=1 a fixed, nominal earth radius is used.
            for iopt=2 the local radius of the earth is used.
    outputs:
        x,y    (r): rectangular coordinates in km

    see "map projections used by the u.s. geological survey"
    geological survey bulletin 1532, pgs 157-173
    for this routine, a spherical earth is assumed for the projection
    the 1972 wgs ellipsoid model (bulletin pg 15).
    the error will be small for small-scale maps.
    """

    radearth = 6378.135  # equitorial earth radius
    f = 298.26  # 1/f wgs 72 model values
    dtr = 3.141592654 / 180.0

    lon1 = mod(lon + 720.0, 360.0)
    orglon1 = mod(orglon + 720.0, 360.0)
    #
    # compute local radius of the earth at center of image
    #
    eradearth = 6378.0  # use fixed nominal value
    if iopt == 2:  # local radius
        era = 1.0 - 1.0 / f
        eradearth = (
            radearth
            * era
            / sqrt(era * era * cos(orglat * dtr) ** 2 + sin(orglat * dtr) ** 2)
        )

    denom = (
        1.0
        + sin(orglat * dtr) * sin(lat * dtr)
        + cos(orglat * dtr) * cos(lat * dtr) * cos(dtr * (lon1 - orglon1))
    )
    if all(denom > 0.0):
        ak = sqrt(2.0 / denom)
    else:
        print("*** division error in lambert1 routine ***")
        ak = 1.0

    x = ak * cos(lat * dtr) * sin(dtr * (lon1 - orglon1))
    y = ak * (
        cos(dtr * orglat) * sin(dtr * lat)
        - sin(dtr * orglat) * cos(dtr * lat) * cos(dtr * (lon1 - orglon1))
    )
    x = x * eradearth
    y = y * eradearth
    return x, y


def polster(alon, alat, xlam, slat):
    """Polar stereographic trasnformation
    (x,y)=polster(lon,lat,xlam,slat)

    computes the polar sterographic transformation for a lon,lat
    input of (alon,alat) with reference origin  lon,lat=(xlam,slat).
    output is (x,y) in km


    algorithm is the same as used for processing ers-1 sar images
    as received from m. drinkwater (1994)
    """
    # ported from polster.m by JPB 21 Sept 2011
    e2 = 0.006693883
    re = 6378.273
    dtr = 3.141592654 / 180.0
    e = sqrt(e2)
    if slat < 0:
        sn = -1.0
        rlat = -alat
    else:
        sn = 1.0
        rlat = alat

    t = ((1.0 - e * sin(rlat * dtr)) / (1.0 + e * sin(rlat * dtr))) ** (e * 0.5)
    ty = tan(dtr * (45.0 - 0.5 * rlat)) / t
    if slat < 0:
        rlat = -slat
    else:
        rlat = slat

    t = ((1.0 - e * sin(dtr * rlat)) / (1.0 + e * sin(dtr * rlat))) ** (e * 0.5)
    tx = tan(dtr * (45.0 - 0.5 * rlat)) / t
    cm = cos(dtr * rlat) / sqrt(1.0 - e2 * sin(dtr * rlat) ** 2)
    rho = re * cm * ty / tx
    x = (sn * sin(dtr * (sn * alon - xlam))) * rho
    y = -(sn * cos(dtr * (sn * alon - xlam))) * rho
    return x, y


def easegrid(iopt, alat, alon, ascale):
    """EASE grid transformation
    (thelon thelat)=easegrid(iopt,lat,lon,ascale)

    computes the forward "ease" grid transform

    given a lat,lon (alat,alon) and the scale (ascale) the image
    transformation coordinates (thelon,thelat) are comuted
    using the "ease grid" (version 1.0) transformation given in fortran
    source code supplied by nsidc.

    the radius of the earth used in this projection is imbedded into
    ascale while the pixel dimension in km is imbedded in bscale
    the base values are: radius earth= 6371.228 km
                 pixel dimen =25.067525 km
    then, bscale = base_pixel_dimen
          ascale = radius_earth/base_pixel_dimen

    iopt is ease type: iopt=11=north, iopt=12=south, iopt=13=cylindrical
    """
    # ported from easegrid.m by JPB 21 Sept 2011
    pi2 = np.pi / 2.0
    dtr = pi2 / 90.0

    if iopt == 11:  # ease grid north
        thelon = ascale * sin(alon * dtr) * sin(dtr * (45.0 - 0.5 * alat))
        thelat = ascale * cos(alon * dtr) * sin(dtr * (45.0 - 0.5 * alat))
    elif iopt == 12:  # ease grid south
        thelon = ascale * sin(alon * dtr) * cos(dtr * (45.0 - 0.5 * alat))
        thelat = ascale * cos(alon * dtr) * cos(dtr * (45.0 - 0.5 * alat))
    elif iopt == 13:  # ease cylindrical
        thelon = ascale * pi2 * alon * cos(30.0 * dtr) / 90.0
        thelat = ascale * sin(alat * dtr) / cos(30.0 * dtr)

    return thelon, thelat


def ease2grid(iopt, alat, alon, ascale, bscale):
    """EASE2 grid transformation

    (thelon thelat)=ease2grid(iopt,lat,lon,ascale,bscale)

    given a lat,lon (alat,alon) and the scale (ascale) the image
        transformation coordinates (thelon,thelat) are comuted
    using the "ease2 grid" (version 2.0) transformation given in IDL
    source code supplied by MJ Brodzik

    RADIUS EARTH=6378.137 KM (WGS 84)
    MAP ECCENTRICITY=0.081819190843 (WGS84)

    inputs:
      iopt: projection type 8=EASE2 N, 9-EASE2 S, 10=EASE2 T/M
      alon, alat: lon, lat (deg) to convert (can be outside of image)
          ascale and bscale should be integer valued)
      ascale: grid scale factor (0..5)  pixel size is (bscale/2^ascale)
      bscale: base grid scale index (ind=int(bscale))

          see ease2helper.py for definitions of isc and ind

    outputs:
      thelon: X coordinate in pixels (can be outside of image)
      thelat: Y coordinate in pixels (can be outside of image)
    """

    DTR = 0.01745329241994
    ind = round(bscale)
    isc = round(ascale)
    dlon = alon
    phi = DTR * alat
    lam = dlon

    # get base EASE2 map projection parameters
    (
        map_equatorial_radius_m,
        map_eccentricity,
        e2,
        map_reference_latitude,
        map_reference_longitude,
        map_second_reference_latitude,
        sin_phi1,
        cos_phi1,
        kz,
        map_scale,
        bcols,
        brows,
        r0,
        s0,
        epsilon,
    ) = ease2_map_info(iopt, isc, ind)

    dlon = dlon - map_reference_longitude
    dlon = easeconv_normalize_degrees(dlon)
    lam = DTR * dlon

    sin_phi = np.sin(phi)

    q = (1.0 - e2) * (
        (sin_phi / (1.0 - e2 * sin_phi * sin_phi))
        - (1.0 / (2.0 * map_eccentricity))
        * np.log(
            (1.0 - map_eccentricity * sin_phi) / (1.0 + map_eccentricity * sin_phi)
        )
    )

    if iopt == 8:  # EASE2 grid north
        qp = 1.0 - (
            (1.0 - e2)
            / (2.0 * map_eccentricity)
            * np.log((1.0 - map_eccentricity) / (1.0 + map_eccentricity))
        )
        rho = map_equatorial_radius_m * np.sqrt(qp - q)
        if np.size(rho) > 1:
            rho[np.abs(qp - q) < epsilon] = 0.0
        else:
            if np.abs(qp - q) < epsilon:
                rho = 0

        x = rho * np.sin(lam)
        y = -rho * np.cos(lam)

    elif iopt == 9:  # EASE2 grid south
        qp = 1.0 - (
            (1.0 - e2)
            / (2.0 * map_eccentricity)
            * np.log((1.0 - map_eccentricity) / (1.0 + map_eccentricity))
        )
        rho = map_equatorial_radius_m * np.sqrt(qp + q)
        if np.size(rho) > 1:
            rho[np.abs(qp - q) < epsilon] = 0.0
        else:
            if np.abs(qp - q) < epsilon:
                rho = 0

        x = rho * np.sin(lam)
        y = rho * np.cos(lam)

    elif iopt == 10:  # EASE2 cylindrical
        x = map_equatorial_radius_m * kz * lam
        y = (map_equatorial_radius_m * q) / (2.0 * kz)

    else:
        print("*** invalid EASE2 projection specificaion in ease2grid")

    thelon = r0 + (x / map_scale) + 0.5
    thelat = s0 + (y / map_scale) + 0.5

    return (thelon, thelat)

这段Python代码依赖于另一个名为ease2helper的python代码。我在同一FTP服务器上发现的。

我也修复了这个文件上的几个衬里问题。这是ease2helper代码；

#!/usr/bin/env python
""" EASE2 grid helper utility functions for sirpy"""

######
# Imports
######
from __future__ import division
import numpy as np

##############################
# The EASE2 helper functions
##############################


def easeconv_normalize_degrees(dlon):
    #
    # Return dlon to within the range -180 <= dlon <= 180
    #  can handle array inputs
    #
    out = dlon
    if np.size(out) > 1:
        while (out < -180.0).sum() > 0:
            out[out < -180.0] = out[out > 180.0] + 360.0
        while (out > 180.0).sum() > 0:
            out[out > 180.0] = out[out > 180.0] - 360.0
    else:
        while out < -180.0:
            out = out + 360.0
        while out > 180.0:
            out = out - 360.0
    return out


def ease2_map_info(iopt, isc, ind):
    """ internally used routine

 (map_equatorial_radius_m, map_eccentricity, \
  e2, map_reference_latitude, map_reference_longitude, \
  map_second_reference_latitude, sin_phi1, cos_phi1, kz, \
  map_scale, bcols, brows, r0, s0, epsilon) = ease2_map_info(iopt, isc, nd)

       defines EASE2 grid information

    inputs:
         iopt: projection type 8=EASE2 N, 9=EASE2 S, 10=EASE2 T/M
         isc:  scale factor 0..5 grid size is (basesize(ind))/2^isc
         ind:  base grid size index   (map units per cell in m

        NSIDC .grd file for isc=0
           project type    ind=0     ind=1         ind=2       ind=3
          N         EASE2_N25km EASE2_N30km EASE2_N36km  EASE2_N24km
              S         EASE2_S25km EASE2_S30km EASE2_S36km  EASE2_S24km
              T/M       EASE2_T25km EASE2_M25km EASE2_M36km  EASE2_M24km

          cell size (m) for isc=0 (scale is reduced by 2^isc)
           project type  ind=0     ind=1            ind=2        ind=3
          N        25000.0     30000.0         36000.0      24000.0
              S        25000.0     30000.0         36000.0      24000.0
              T/M   T25025.26 M25025.26000 M36032.220840584 M24021.480560389347

      for a given base cell size (e.g., ind=0) isc is related to
          NSIDC CETB .grd file names according to
         isc        N .grd name   S .grd name   T .grd name
          0       EASE2_N25km     EASE2_S25km     EASE2_T25km
          1       EASE2_N12.5km   EASE2_S12.5km   EASE2_T12.5km
          2       EASE2_N6.25km   EASE2_S6.25km   EASE2_T6.25km
          3       EASE2_N3.125km  EASE2_S3.125km  EASE2_T3.125km
          4       EASE2_N1.5625km EASE2_S1.5625km EASE2_T1.5625km

  outputs
    map_equatorial_radius_m  EASE2 Earth equitorial radius (km) [WGS84]
    map_eccentricity         EASE2 Earth eccentricity [WGS84]
    map_reference_latitude   Reference latitude (deg)
    map_reference_longitude  Reference longitude (deg)
    map_second_reference_latitude Secondary reference longitude* (deg)
    sin_phi1, cos_phi1 kz    EASE2 Cylin parameters*
    map_scale                EASE2 map projection pixel size (km)
    bcols, brows,            EASE2 grid size in pixels
    r0, s0                   EASE2 base projection size in pixels
    epsilon                  EASE2 near-polar test factor
    """

    DTR = 0.01745329241994

    m = 2 ** np.floor(isc)  # compute power-law scale factor

    map_equatorial_radius_m = 6378137.0  # WGS84
    map_eccentricity = 0.081819190843  # WGS84
    e2 = map_eccentricity * map_eccentricity
    map_reference_longitude = 0.0
    epsilon = 1.0e-6

    #  map-specific parameters
    if iopt == 8:  # EASE2 grid north
        map_reference_latitude = 90.0
        if ind == 1:  # EASE2_N30km.gpd
            base = 30000.0
            nx = 600
            ny = 600
        elif ind == 2:  # EASE2_N36km.gpd
            base = 36000.0
            nx = 500
            ny = 500
        elif ind == 3:  # EASE2_N24km.gpd
            base = 24000.0
            nx = 750
            ny = 750
        else:  # EASE2_N25km.gpd
            base = 25000.0
            nx = 720
            ny = 720
        map_second_reference_latitude = 0.0
        sin_phi1 = 0.0
        cos_phi1 = 1.0
        kz = cos_phi1

    elif iopt == 9:  # EASE2 grid south
        map_reference_latitude = -90.0
        if ind == 1:  # EASE2_S30km.gpd
            base = 30000.0
            nx = 600
            ny = 600
        elif ind == 2:  # EASE2_S36km.gpd
            base = 36000.0
            nx = 500
            ny = 500
        elif ind == 3:  # EASE2_S24km.gpd
            base = 24000.0
            nx = 750
            ny = 750
        else:  # EASE2_S25km.gpd
            base = 25000.0
            nx = 720
            ny = 720
        map_second_reference_latitude = 0.0
        sin_phi1 = 0.0
        cos_phi1 = 1.0
        kz = cos_phi1

    elif iopt == 10:  # EASE2 cylindrical
        map_reference_latitude = 0.0
        map_second_reference_latitude = 30.0
        sin_phi1 = np.sin(DTR * map_second_reference_latitude)
        cos_phi1 = np.cos(DTR * map_second_reference_latitude)
        kz = cos_phi1 / np.sqrt(1.0 - e2 * sin_phi1 * sin_phi1)
        if ind == 1:  # EASE2_M25km.gpd
            base = 25025.26000
            nx = 1388
            ny = 584
        elif ind == 2:  # EASE2_M36km.gpd
            base = 36032.220840584
            nx = 964
            ny = 406
        elif ind == 3:  # EASE2_M24km.gpd
            base = 24021.480560389347
            nx = 1446
            ny = 609
        else:  # EASE2_T25km.gpd
            base = 25025.26000
            nx = 1388
            ny = 540
    else:
        print("*** invalid EASE2 projection code ***")

    # grid info
    if isc >= 0:
        map_scale = base / m
        bcols = np.ceil(nx * m)
        brows = np.ceil(ny * m)
        r0 = (nx * m - 1) / 2
        s0 = (ny * m - 1) / 2
    else:
        map_scale = base * m
        bcols = np.ceil(nx / np.float(m))
        brows = np.ceil(ny / np.float(m))
        r0 = (nx / np.float(m) - 1) / 2
        s0 = (ny / np.float(m) - 1) / 2

    return (
        map_equatorial_radius_m,
        map_eccentricity,
        e2,
        map_reference_latitude,
        map_reference_longitude,
        map_second_reference_latitude,
        sin_phi1,
        cos_phi1,
        kz,
        map_scale,
        bcols,
        brows,
        r0,
        s0,
        epsilon,
    )

Answer 2

来自CATDS的SMOS数据也有同样的问题。这就是为什么我为python创建了龙拉特包，它将地理坐标(经度、纬度)转换为(2)-grid坐标(col、row)，反之亦然。它使用吡普罗季 (PROJ)库将经度和纬度转换为EASE(2)投影的x和y坐标。然后使用NSIDC中的参数计算特定网格的行和列索引。

我已经用EASE2 2测试了它，分辨率为9公里(北、南和全球预测)和36公里(北和全球)的SMAP数据。

看起来，r也支持PROJ，所以您可以尝试同样的方法。