Как ускорить создание регулярной сетки в чанках и медленная запись в файл - PullRequest
Новая
       9

Как ускорить создание регулярной сетки в чанках и медленная запись в файл

0 голосов
/ 08 июня 2018

Я пытаюсь интерполировать разбросанные точки в регулярную сетку.Для небольших доменов, т. Е. Небольшое количество строк и столбцов np.meshgrid работает нормально.Если строки и столбцы большие, выдается MemoryError.Поэтому я попытался обработать весь домен небольшими порциями и применить функцию интерполяции и записать ее в файл geotiff с помощью gdal.
Ниже приведен код, и я дал пояснения в комментариях. 

import numpy as np
from osgeo import gdal
import csv
import scipy.spatial as spatial

## loading lat, lon and values from csv file
lat = []
lon = []
ele = []

## The csv file contains the lat/lon/ele 
with open('data.csv', 'r') as data:
    for row in data:
        row = row.strip().split()
        lat.append(float(row[2]))
        lon.append(float(row[1]))
        ele.append(float(row[3]))
## creating a numpy array to feed into KDTree for spatial indexing      
xycoord = np.c_[lon,lat]
ele_arr = np.array(ele)

## Generating KDTree for nearest neighour search
point_tree = spatial.cKDTree(xycoord, leafsize=15)

## Getting domain extents
## 23.204 , 146.447, -61.509, 25.073
xmin, xmax, ymin, ymax = min(lon),max(lon), min(lat), max(lat)  

res = 0.003 ## Grid spacing ~~330 meters 

x = np.arange(xmin, xmax, res, dtype=np.float16)
y = np.arange(ymin, ymax, res, dtype=np.float16)

nx = x.shape[0]
ny = y.shape[0]
print (nx, ny) # ~ (41081 28861)

## Creating of geotiff file using gdal
outFile = "test.tif"
format_ = "GTiff"
driver = gdal.GetDriverByName(format_)
outRas = driver.Create(outFile, nx, ny, 1, gdal.GDT_Float32, options=['COMPRESS=DEFLATE'])
outRas.SetGeoTransform((xmin, res, 0, ymax, 0, -res))

## No of rows and columns in each chunk
step = 2000 

## starting and ending indices for row and column for each chunk
xstart = []
xend = []
ystart = []
yend = []

for i in range(0,nx,step):
    for j in range(0,ny,step):
        xstart.append(i)
        xend.append(i+step)
        ystart.append(j)
        yend.append(j+step)

## Actual loop      
for i in range(len(xstart)):
    t = np.meshgrid(x[xstart[i]:xend[i]],y[ystart[i]:yend[i]]) ## Creating a meshgrid 
    ## Actual intended flow
    ## xy = np.stack(np.meshgrid(x[xstart[i]:xend[i]],y[ystart[i]:yend[i]]), axis = -1)
    ## distances, points_idx = point_tree.query(xy, k=3, eps=0)
    ## z = interpFn(distances, ele_arr[points_idx])

    ##To test the speed, not using above code and using a simple fn which 
    ##takes our input matrix and return matrix with same dimensions. even np.ones() will do

    z = fn(t) ## this could be any function
    outRas.GetRasterBand(1).WriteArray(z,xstart[i],ystart[i]) ## Writing to geotiff file 

outRas = None

Теперь *Разрешение 1007 * разрешено, но для матриц большого размера это очень медленно и занимает много времени, чтобы написать простую матрицу даже до применения какой-либо функции. Пожалуйста, предоставьте предложения по ускорению процесса

1 Ответ

0 голосов
/ 09 июня 2018

Относительно вашей первоначальной проблемы (в зависимости от содержимого ваших данных) установка параметра sparce на True потенциально может сэкономить вам много памяти: 

np.meshgrid(data, data, sparse=True)
...