Python ve gdal kullanılarak nokta verilerinin IDW enterpolasyonu

Question

Lat, Long ve Rainfall Bilgilerine sahip bir CSV dosyasına sahibim. Bu noktaya enterpolasyon yapmak ve tiff dosyası oluşturmak istiyorum. Herhangi biri bana bunu yapmanın en kolay yolunu önerebilir.

gdal_grid'i kullanmaya çalışıyorum. Python'da gdal kullanmanın çok yeni.

Answer 1

Bu aslında birkaç soru. Lateksler ve uzunlar için bazı dağınık verilere sahip olduğunuzu varsayarak tüm konumunuzu tahmin edersiniz (Tiff imajınızın pikselleri için tüm lat ve longs).

Eğer (başka question bir son örneği kullanarak) veriler üzerinde IWD yapmak etrafında çözümlerin herhangi kullanabileceği olduğunda:

class Estimation(): 
    # IWD. Check: https://stackoverflow.com/questions/36031338/interpolate-z-values-in-a-3d-surface-starting-from-an-irregular-set-of-points/36037288#36037288 
    def __init__(self,lon,lat,values): 
     self.x = lat 
     self.y = lon 
     self.v = values 

    def estimate(self,x,y,using='ISD'): 
     """ 
     Estimate point at coordinate x,y based on the input data for this 
     class. 
     """ 
     if using == 'ISD': 
      return self._isd(x,y) 

    def _isd(self,x,y): 
     #d = np.sqrt((x-self.x)**2+(y-self.y)**2) 
     d = x.copy() 
     for i in range(d.shape[0]): 
      d[i] = haversine(self.x[i],self.y[i],x,y) 
     if d.min() > 0: 
      v = np.sum(self.v*(1/d**2)/np.sum(1/d**2)) 
      return v 
     else: 
      return self.v[d.argmin()] 

Yukarıdaki kod aslında ile mesafe hesaplamak için uyarlanmıştır

Haversine formula (bu, bir kürenin uzunlamasına ve enlemlerinden iki nokta arasında büyük daire mesafeleri verir). Bildirim tekrar this one gibi haversinüs mesafe için çözümlerin her türlü bulabilirsiniz:

def haversine(lon1, lat1, lon2, lat2): 
    """ 
    Check: https://stackoverflow.com/questions/15736995/how-can-i-quickly-estimate-the-distance-between-two-latitude-longitude-points 
    Calculate the great circle distance between two points 
    on the earth (specified in decimal degrees) 
    """ 
    # convert decimal degrees to radians 
    lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2]) 
    # haversine formula 
    dlon = lon2 - lon1 
    dlat = lat2 - lat1 
    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2 
    c = 2 * asin(sqrt(a)) 
    km = 6367 * c 
    return km 

Son olarak size sadece Tiff kullanarak GDAL inşa etmeliyiz diziniz hazır bir kez. Bunun için aşağıdaki çözümün bir parçası olan aşağıdaki question kontrolünü yapın:

driver = gdal.GetDriverByName('GTiff') 
ds = driver.Create('output.tif',xsize, ysize, 1, gdal.GDT_Float32,) 
# this assumes the projection is Geographic lat/lon WGS 84 
srs = osr.SpatialReference() 
srs.ImportFromEPSG(4326) 
ds.SetProjection(srs.ExportToWkt()) 
gt = [ulx, xres, 0, uly, 0, yres ] 
ds.SetGeoTransform(gt) 
outband=ds.GetRasterBand(1) 
outband.SetStatistics(np.min(mag_grid), np.max(mag_grid), np.average(mag_grid), np.std(mag_grid)) 
outband.WriteArray(mag_grid)