numpy 이동 창 퍼센트 커버

Question

나는 분류가 작은 배열로 읽고있는 분류 래스터를 가지고 있습니다. (n 클래스)

창이있는 각 클래스의 % 커버를 저장하는 n 차원 벡터를 만들기 위해 2 차원 이동 윈도우 (예 : 3 x 3)를 사용하고 싶습니다. 래스터가 크기 때문에이 정보를 매번 다시 계산하지 않도록 저장하는 것이 유용 할 것입니다. 따라서 최상의 솔루션은 벡터로 작동하는 3D 배열을 만드는 것이라고 생각합니다. 이러한 %/개수 값을 기준으로 새 래스터가 만들어집니다.

1) 차원 어레이 N + 1 '밴드'

2) 대역 1 = 원래 장비 래스터 작성할 :

내 아이디어이다. 서로 '밴드'값 = 예 .... 창 (클래스마다 즉, 하나 개의 밴드) 내의 값의 셀 수 :

[[2 0 1 2 1] 
[2 0 2 0 0] 
[0 1 1 2 1] 
[0 2 2 1 1] 
[0 1 2 1 1]] 
[[2 2 3 2 2] 
[3 3 3 2 2] 
[3 3 2 2 2] 
[3 3 0 0 0] 
[2 2 0 0 0]] 
[[0 1 1 2 1] 
[1 3 3 4 2] 
[1 2 3 4 3] 
[2 3 5 6 5] 
[1 1 3 4 4]] 
[[2 3 2 2 1] 
[2 3 3 3 2] 
[2 4 4 3 1] 
[1 3 5 3 1] 
[1 3 3 2 0]] 

4)을 VRT로 그렇게 만 요구 이들 밴드 판독 생성 할 것을 한 번 ... 추가 모듈을 위해 읽을 수 있습니다.

질문 : 창 안에서 '계산하는'가장 효율적인 '움직이는 창'방법은 무엇입니까?

현재 - 나는 노력하고, 다음과 같은 코드로 실패하고있다 :

def lcc_binary_vrt(raster, dim, bands): 
    footprint = np.zeros(shape = (dim,dim), dtype = int)+1 
    g = gdal.Open(raster) 
    data = gdal_array.DatasetReadAsArray(g) 

    #loop through the band values 
    for i in bands: 
     print i 
     # create a duplicate '0' array of the raster 
     a_band = data*0 
     # we create the binary dataset for the band   
     a_band = np.where(data == i, 1, a_band) 
     count_a_band_fname = raster[:-4] + '_' + str(i) + '.tif'   
     # run the moving window (footprint) accross the band to create a 'count' 
     count_a_band = ndimage.generic_filter(a_band, np.count_nonzero(x), footprint=footprint, mode = 'constant') 
     geoTiff.create(count_a_band_fname, g, data, count_a_band, gdal.GDT_Byte, np.nan) 

어떤 제안이 아주 많이 감사합니다.

베키

Answer 1

나는 공간 과학 물건에 대해 아무것도 몰라, 그래서 그냥 '창을 이동'가장 효율적인 방법 '이 무엇인지 주요 질문 :)

에 초점을 맞출 것 창문 안에서 '카운트'하시겠습니까?

NumPy와 함께 창 통계를 이동 할 수있는 일반적인 방법은, numpy.lib.stride_tricks.as_strided를 사용하는 예를 this answer를 참조하는 것입니다. 당신이 반복해서 같은 번호를 합산하고 있기 때문에,

from numpy.lib.stride_tricks import as_strided 

... 

m, n = a_band.shape 
newshape = (m-dim+1, n-dim+1, dim, dim) 
newstrides = a_band.strides * 2 # strides is a tuple 
count_a_band = as_strided(ar, newshape, newstrides).sum(axis=(2,3)) 

그러나 귀하의 사용 경우에이 방법은 비효율적이다 : 기본적으로, 아이디어는 메모리 사용량의 증가없이, 모든 윈도우가 포함되고있는 배열을 만드는 것입니다 특히 창 크기가 커질 경우 다시 한 번 강조합니다. 더 좋은 방법은 this answer에서처럼 cumsum 트릭을 사용하는 것입니다 : 위 두 코드의 가장자리 경우를 처리하는 지루한 될 것이라고

가

def windowed_sum_1d(ar, ws, axis=None): 

    if axis is None: 
     ar = ar.ravel() 
    else: 
     ar = np.swapaxes(ar, axis, 0) 

    ans = np.cumsum(ar, axis=0) 
    ans[ws:] = ans[ws:] - ans[:-ws] 

    ans = ans[ws-1:] 

    if axis: 
     ans = np.swapaxes(ans, 0, axis) 

    return ans 


def windowed_sum(ar, ws): 
    for axis in range(ar.ndim): 
     ar = windowed_sum_1d(ar, ws, axis) 
    return ar 

... 

count_a_band = windowed_sum(a_band, dim) 

참고. 다행히,이 포함하고 두 번째 코드와 동일한 효율을 얻을 수있는 쉬운 방법이있다 :

count_a_band = ndimage.uniform_filter(a_band, size=dim, mode='constant') * dim**2 

이미 한 것과 매우 유사한 비록를이 훨씬 빨라집니다! 단점은 부동 소수점 반올림 오류를 제거하기 위해 정수로 반올림해야 할 수도 있다는 것입니다. 최종 참고로

, 코드는

# create a duplicate '0' array of the raster 
a_band = data*0 
# we create the binary dataset for the band   
a_band = np.where(data == i, 1, a_band) 

비트 중복 : 당신은 a_band = (data == i)를 사용할 수 있습니다.