2

I have three datasets now. One is SMOS level3 soil moisture dataset which is in EASE grid 2.0 ~25km cylindrical projection. (EPSG: 6933). TRMM precipitation dataset which is in 0.25 degree. SMAP soil moisture dataset which is in EASE grid ~3km cylindrical projection. Now I want to do some data mining work so I want to downscale the resolution of the SMOS and TRMM dataset to be the same as the SMAP dataset. All the datasets are only for Oklahoma state. Say, the SMOS dataset is like a 22 * 43 matrix, now I want to expand it to a 145 * 338 matrix with interpolation with the latitude and longitude being considered.

I am currently trying to use cKDTree in python to find k nearest neighbors for one specific grid. I use the code from this link (the blog only describes the methods but it has an associate notebook with code there). The cKDTree can return four points. Then I can use the latitude and longitude to compute the distance in order to do the bilinear interpolation. Here is the code:

def distance(lat_center, lon_center, lat0, lon0):
    rad_factor = pi/180.0
    lat_center *= rad_factor
    lon_center *= rad_factor
    lat0_rad = lat0 * rad_factor
    lon0_rad = lon0 * rad_factor

    clat, clon = cos(lat_center), cos(lon_center)
    slat, slon = sin(lat_center), sin(lon_center)
    delX = cos(lat0_rad) * cos(lon0_rad) - clat * clon
    delY = cos(lat0_rad) * sin(lon0_rad) - clat * slon
    delZ = sin(lat0_rad) - slat
    dist_sq = delX ** 2 + delY ** 2 + delZ ** 2
    return dist_sq

def bilinear(lat_center, lon_center, iy, ix, latvar, lonvar, target_var):
    area = []
    for x, y, i in zip(ix, iy, range(len(iy))):
        area.append(distance(lat_center, lon_center, latvar[y], lon_center) * distance(lat_center, lon_center, lat_center, lonvar[x]))
    print area

I do process the data a little bit before construct the KDTree like this:

class Kdtree_fast(object):
def __init__(self, ncfile, latvarname, lonvarname, target_varname):
    self.ncfile = ncfile
    self.latvar = self.ncfile.variables[latvarname]
    self.lonvar = self.ncfile.variables[lonvarname]
    self.target_var = self.ncfile.variables[target_varname]

    rad_factor = pi / 180.0
    self.latvals = self.latvar[:] * rad_factor
    self.lonvals = self.lonvar[:] * rad_factor
    ny = len(self.latvals)
    nx = len(self.lonvals)
    self.latvals = np.repeat(self.latvals, nx, axis=0).reshape(ny, nx)
    self.lonvals = np.tile(self.lonvals, ny).reshape(ny, nx)
    self.shape = self.latvals.shape

    clat, clon = cos(self.latvals), cos(self.lonvals)
    slat, slon = sin(self.latvals), sin(self.lonvals)
    clat_clon = clat * clon
    clat_slon = clat * slon
    triples = list(zip(np.ravel(clat_clon), np.ravel(clat_slon), np.ravel(slat)))
    self.kdt = cKDTree(triples)

def query(self, lat0, lon0):
    rad_factor = pi / 180.0
    lat0_rad = lat0 * rad_factor
    lon0_rad = lon0 * rad_factor
    clat0, clon0 = cos(lat0_rad), cos(lon0_rad)
    slat0, slon0 = sin(lat0_rad), sin(lon0_rad)
    dist_sq_min, minindex_1d = self.kdt.query(
        [clat0 * clon0, clat0 * slon0, slat0], k=4)
    iy_min, ix_min = np.unravel_index(minindex_1d, self.shape)
    return dist_sq_min, iy_min, ix_min

The dataset offers the coordinate information like this:

latitude: 33.04897 33.28237 ... 37.844597 38.092167

longitude: -103.87608 -103.616714 ... -93.24207 -92.98271

I want to implement the bilinear by myself because later I may use some machine learning method to do the interpolation which will still needs to compute the distance between grids correctly.

Now my trouble is the dist_sq_min returned by the cKDTree doesn't match the distance returned by the distance function if I call them manually.

The reason I need to write the distance function is the cKDTree only returns the distance between the target point and the neighbors while the bilinear needs the distance vertically and horizontally to compute the area as the weight.

  • @PolyGeo I am new to this field. I googled this problem and found ArcGIS seems to be able to deal with my issue. But it cannot be installed on my mac (neither will I install windows system). Do you have any suggestions? I am currently trying to use cKDTree in python to find k nearest neighbors for one specific grid. Then get a weighted value for it. Hopefully it will work. – HannaMao Jul 16 '16 at 1:30
  • 1
    If you are trying to do it using cKDTree then I think you should edit your question to focus it on that. For focussed Q&A the idea is to tell us precisely what you want to do and what you have tried, then ask about where you are stuck. – PolyGeo Jul 16 '16 at 1:41
1

You can use the KD-tree based method that you propose. I've split the answer into two parts, addressing your questions about projection and how to interpolate with the KDTree.

Projection

Provided that the SMOS and SMAP datasets are in a cylindrical projected coordinate system, you will need to first convert them to longitude and latitude to use the linked method. pyproj is a good way to do this in Python, or read https://en.wikipedia.org/wiki/Cylindrical_equal-area_projection. However, if your data are gridded to the EASE grid but have longitude and latitude data already associated with with grid cell, then it's done for you.

Interpolation

I will assume that cKDTree refers to the scipy implementation of a KD-tree.

You seem to be converting the data to 3d cartesian space during your preprocessing step, which is fine. When you query scipy.spatial.cKDTree, it will return a tuple of distances and indices, so there is no need to recompute distance with your distance function.

I can't see how you're performing the queries, but here's where I think your trouble might lie. In general you can't mix scipy's cKDTree with geographical coordinates, which is why the Unidata Notebook you link to converts to 3d cartesian coordinates first. Your KD-tree is now build with points in 3d cartesian space (assuming a spherical earth). The points that you query it with will have to also be in 3d cartesian space, not geographical lon/lat.

Since you're only doing bilinear interpolation, an alternative solution to consider (which I'm speculating will be much higher performance) is to convert all of your data to the target cylindrical coodinate system rather than lon/lat, and to use scipy.interpolate.griddata to perform the interpolation.

I noticed that you tile your coordinate data. One pitfall to using regularly-spaced geographical coordinates is that they are not regularly spaced in cartesian space. For example, the lon/lat points (0, 0) and (1, 0) are about 111 km apart, but the points (0, 60) and (1, 60) are much closer. As a result, where you expect for the position * in the diagram below ot be nearest to points 1, 2, 4, and 5

1       2       3
      * 

4       5       6

it might be the case that 1, 2, 3 and 5 (or some other set) are actually closest because the horizontal and vertical spacing are not the same. This would break your bilinear interpolation function. If you're going to rely on k nearest points with the uneven grid you get from lon/lat values, your interpolation may need to handle this. There are interpolation techniques that do, including IDW and kriging.

Edit

Here is a mocked version of the code in question that seems to work.

class Kdtree_fast(object):

    def __init__(self, lat, lon):
        self.latvals = np.radians(lat)
        self.lonvals = np.radians(lon)
        ny = len(self.latvals)
        nx = len(self.lonvals)
        self.latvals = np.repeat(self.latvals, nx, axis=0).reshape(ny, nx)
        self.lonvals = np.tile(self.lonvals, ny).reshape(ny, nx)
        self.shape = self.latvals.shape

        clat, clon = np.cos(self.latvals), np.cos(self.lonvals)
        slat, slon = np.sin(self.latvals), np.sin(self.lonvals)
        triples = list(zip(np.ravel(clat*clon),
                           np.ravel(clat*slon),
                           np.ravel(slat)))
        self.kdt = cKDTree(triples)

    def query(self, lat0, lon0):
        lat0_rad = np.radians(lat0)
        lon0_rad = np.radians(lon0)
        clat0, clon0 = np.cos(lat0_rad), np.cos(lon0_rad)
        slat0, slon0 = np.sin(lat0_rad), np.sin(lon0_rad)
        dist_sq_min, minindex_1d = self.kdt.query(
            [clat0 * clon0, clat0 * slon0, slat0], k=4)
        iy_min, ix_min = np.unravel_index(minindex_1d, self.shape)
        return dist_sq_min, iy_min, ix_min

def distance(lat_center, lon_center, lat0, lon0):
    lat_center = np.radians(lat_center)
    lon_center = np.radians(lon_center)
    lat0_rad = np.radians(lat0)
    lon0_rad = np.radians(lon0)

    clat, clon = np.cos(lat_center), np.cos(lon_center)
    slat, slon = np.sin(lat_center), np.sin(lon_center)
    delX = np.cos(lat0_rad) * np.cos(lon0_rad) - clat * clon
    delY = np.cos(lat0_rad) * np.sin(lon0_rad) - clat * slon
    delZ = np.sin(lat0_rad) - slat
    dist_sq = delX ** 2 + delY ** 2 + delZ ** 2
    return dist_sq

lats = np.arange(33, 38, 0.2)
lons = np.arange(-103, -92, 0.3)
kd = Kdtree_fast(lats, lons)

# should print four zeros
d, i, j = kd.query(35.1, -100.2)
for _d, _i, _j in zip(d, i, j):
    print(_d - np.sqrt(distance(35.1, -100.2, lats[_i], lons[_j])))
  • Thanks for the answer! I think I will try cKDTree first to implement this from scratch as it gives me more flexibility. I am totally new to this area and quite cautious on every step I take. It did take me long time to find cKDTree. Would you please tell me whether I am in the right direction or not? – HannaMao Jul 16 '16 at 15:30
  • @HannaMao A KD-tree is a data structure that may be useful for implementing various interpolation schemes, but it's an implementation detail. The key is ensuring that you are interpolating points at the right positions on the relevant coordinate system, which will require reprojection between the EASE grid and geographical longitude and latitude. – Nat Wilson Jul 16 '16 at 16:17
  • The smos dataset is in EASE grid 2.0 with WGS84. Do I still need to reproject it? I thought it is already in geographical latitude and longitude. – HannaMao Jul 16 '16 at 17:26
  • I believe you will need to reproject. Two things are important. (1) Cylindrical coordinates (i.e. the EASE grid) are not geographical coordinates; they are projected to a cylindrical coordinate system. (2) You probably don't want to use geographical coordinates anyway. cKDTree does not know how to represent positions on an ellipsoid surface as needed for geographical coordinates. This is why map projections are useful - they project geographical coordinates to a more manageable space. I would read up on map projections and coordinate systems if you are intent on coding this yourself. – Nat Wilson Jul 16 '16 at 17:46
  • I've edited my post with detailed code. I am sorry I don't know the difference between cylindrical coordinate and geographical coordinate and what they mean to my program. Do you think I should do what was discussed here: link? I ignored this before is what was discussed there is for EASE original. Now I am using EASE 2.0 with WGS84. I thought I didn't need to process the coordinates. – HannaMao Jul 16 '16 at 18:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.