@@ -2443,7 +2443,7 @@ def createVirtualRaster(self, inputRasterList, output, bandNumberList = 'No', qu
24432443 offY = 0
24442444 except :
24452445 pass
2446- vRast .AddBand (cfg . gdalSCP . GDT_Float32 )
2446+ vRast .AddBand (gBand2 . DataType )
24472447 bandNumber = bandNumberList [x ]
24482448 band = vRast .GetRasterBand (x + 1 )
24492449 bsize = band .GetBlockSize ()
@@ -2537,8 +2537,8 @@ def createVirtualRaster(self, inputRasterList, output, bandNumberList = 'No', qu
25372537 offY = 0
25382538 except :
25392539 pass
2540+ vRast .AddBand (gBand2 .DataType )
25402541 gBand2 = None
2541- vRast .AddBand (cfg .gdalSCP .GDT_Float32 )
25422542 try :
25432543 errorCheck = 'Yes'
25442544 if bandNumberList == 'No' :
@@ -2982,7 +2982,9 @@ def rasterBlocks(self, gdalRaster, blockSizeX = 1, blockSizeY = 1, previewSize =
29822982
29832983
29842984 # read a block of band as array
2985- def readArrayBlock (self , gdalBand , pixelStartColumn , pixelStartRow , blockColumns , blockRow ):
2985+ def readArrayBlock (self , gdalBand , pixelStartColumn , pixelStartRow , blockColumns , blockRow , calcDataType = None ):
2986+ if calcDataType is None :
2987+ calcDataType = cfg .np .float32
29862988 try :
29872989 o = gdalBand .GetOffset ()
29882990 s = gdalBand .GetScale ()
@@ -2993,8 +2995,12 @@ def readArrayBlock(self, gdalBand, pixelStartColumn, pixelStartRow, blockColumns
29932995 except :
29942996 o = 0.0
29952997 s = 1.0
2998+ # logger
2999+ cfg .utls .logToFile (str (__name__ ) + '-' + str (cfg .inspectSCP .stack ()[0 ][3 ])+ ' ' + cfg .utls .lineOfCode (), 's ' + str (s )+ ' o ' + str (o ))
29963000 try :
2997- a = gdalBand .ReadAsArray (pixelStartColumn , pixelStartRow , blockColumns , blockRow )* s + o
3001+ a = gdalBand .ReadAsArray (pixelStartColumn , pixelStartRow , blockColumns , blockRow ) * cfg .np .asarray (s ).astype (calcDataType ) + cfg .np .asarray (o ).astype (calcDataType )
3002+ # logger
3003+ cfg .utls .logToFile (str (__name__ ) + '-' + str (cfg .inspectSCP .stack ()[0 ][3 ])+ ' ' + cfg .utls .lineOfCode (), 'a ' + str (a [0 ,0 ]))
29983004 except :
29993005 return None
30003006 return a
@@ -3789,15 +3795,26 @@ def calculateRasterWithStats(self, gdalBandList, rasterSCPArrayfunctionBand, col
37893795 return [outputRaster , [[min , max , mean , stdDev ], array , count ]]
37903796
37913797 # calculate raster unique values with sum
3792- def rasterUniqueValuesWithSum (self , gdalBandList , rasterSCPArrayfunctionBand , columnNumber , rowNumber , pixelStartColumn , pixelStartRow , outputArrayFile , functionBandArgumentNoData , functionVariableList , outputBandNumber ):
3793- o = cfg .np .array (cfg .np .unique (rasterSCPArrayfunctionBand [~ cfg .np .isnan (rasterSCPArrayfunctionBand )], return_counts = True ))
3798+ def rasterUniqueValuesWithSum (self , gdalBandList , rasterSCPArrayfunctionBand , nodataMask , rowNumber , pixelStartColumn , pixelStartRow , outputArrayFile , functionBandArgumentNoData , functionVariableList , outputBandNumber ):
3799+ try :
3800+ outputNoData = gdalBandList [2 ]
3801+ b = rasterSCPArrayfunctionBand [nodataMask != outputNoData ]
3802+ o = cfg .np .array (cfg .np .unique (b [~ cfg .np .isnan (b )], return_counts = True ))
3803+ except :
3804+ o = cfg .np .array (cfg .np .unique (rasterSCPArrayfunctionBand [~ cfg .np .isnan (rasterSCPArrayfunctionBand )], return_counts = True ))
37943805 return o
37953806
37963807 # calculate raster unique values
37973808 def rasterUniqueValues (self , gdalBandList , rasterSCPArrayfunctionBand , nodataMask , rowNumber , pixelStartColumn , pixelStartRow , outputArrayFile , functionBandArgumentNoData , functionVariableList , outputBandNumber ):
3798- a = rasterSCPArrayfunctionBand [:,:,0 ].ravel ()
3799- for i in range (1 , rasterSCPArrayfunctionBand .shape [2 ]):
3800- a = cfg .np .vstack ((a ,rasterSCPArrayfunctionBand [:,:,i ].ravel ()))
3809+ outputNoData = gdalBandList [2 ]
3810+ try :
3811+ a = rasterSCPArrayfunctionBand [:,:,0 ][nodataMask != outputNoData ].ravel ()
3812+ for i in range (1 , rasterSCPArrayfunctionBand .shape [2 ]):
3813+ a = cfg .np .vstack ((a ,rasterSCPArrayfunctionBand [:,:,i ][nodataMask != outputNoData ].ravel ()))
3814+ except :
3815+ a = rasterSCPArrayfunctionBand [:,:,0 ].ravel ()
3816+ for i in range (1 , rasterSCPArrayfunctionBand .shape [2 ]):
3817+ a = cfg .np .vstack ((a ,rasterSCPArrayfunctionBand [:,:,i ].ravel ()))
38013818 a = a .T
38023819 a = a [~ cfg .np .isnan (a ).any (axis = 1 )]
38033820 # adapted from Jaime answer at https://stackoverflow.com/questions/16970982/find-unique-rows-in-numpy-array
@@ -3871,7 +3888,7 @@ def crossRasters(self, gdalBandList, rasterSCPArrayfunctionBand, nodataMask, row
38713888 return [o , stats ]
38723889
38733890 # band calculation
3874- def bandCalculation (self , gdalBandList , rasterSCPArrayfunctionBand , columnNumber , rowNumber , pixelStartColumn , pixelStartRow , outputArrayFile , functionBandArgument , functionVariableList , outputBandNumber ):
3891+ def bandCalculation (self , gdalBandList , rasterSCPArrayfunctionBand , nodataMask , rowNumber , pixelStartColumn , pixelStartRow , outputArrayFile , functionBandArgument , functionVariableList , outputBandNumber ):
38753892 f = functionBandArgument
38763893 # logger
38773894 cfg .utls .logToFile (str (__name__ ) + '-' + str (cfg .inspectSCP .stack ()[0 ][3 ])+ ' ' + cfg .utls .lineOfCode (), 'f ' + str (f ) )
@@ -3900,9 +3917,8 @@ def bandCalculation(self, gdalBandList, rasterSCPArrayfunctionBand, columnNumber
39003917 if not isinstance (o , cfg .np .ndarray ):
39013918 return 'No'
39023919 # check nan
3903- #oo = cfg.np.nan_to_num(o) * 1.0
3904- #oo[cfg.np.isnan(o)] = cfg.np.nan
3905- #o = oo
3920+ if nodataMask is not None :
3921+ o [::, ::][nodataMask [::, ::] != 0 ] = nodataMask [::, ::][nodataMask [::, ::] != 0 ]
39063922 return o
39073923
39083924 # multiple where scatter raster calculation
@@ -4788,6 +4804,7 @@ def processRasterDev(self, raster = None, signatureList = None, functionBand = N
47884804 for b in range (0 , len (gdalBandList )):
47894805 sclB = 1
47904806 offsB = 0
4807+ ndvBand = None
47914808 '''
47924809 if nodataValue is None:
47934810 ndv = None
@@ -4801,9 +4818,9 @@ def processRasterDev(self, raster = None, signatureList = None, functionBand = N
48014818 pCount = 0
48024819 while True :
48034820 pCount = pCount + 1
4804- a = cfg .utls .readArrayBlock (gdalBandList [b ], x , y , bSX , bSY )
4821+ a = cfg .utls .readArrayBlock (gdalBandList [b ], x , y , bSX , bSY , calcDataType )
48054822 time .sleep (0.1 )
4806- a2 = cfg .utls .readArrayBlock (gdalBandList [b ], x , y , bSX , bSY )
4823+ a2 = cfg .utls .readArrayBlock (gdalBandList [b ], x , y , bSX , bSY , calcDataType )
48074824 checkA = cfg .np .allclose (a , a2 , equal_nan = True )
48084825 # logger
48094826 cfg .utls .logToFile (str (__name__ ) + '-' + str (cfg .inspectSCP .stack ()[0 ][3 ])+ ' ' + cfg .utls .lineOfCode (), 'checkA ' + str (checkA ) )
@@ -4814,10 +4831,13 @@ def processRasterDev(self, raster = None, signatureList = None, functionBand = N
48144831 break
48154832 a2 = None
48164833 else :
4817- a = cfg .utls .readArrayBlock (gdalBandList [b ], x , y , bSX , bSY )
4834+ a = cfg .utls .readArrayBlock (gdalBandList [b ], x , y , bSX , bSY , calcDataType )
4835+ # logger
4836+ cfg .utls .logToFile (str (__name__ ) + '-' + str (cfg .inspectSCP .stack ()[0 ][3 ])+ ' ' + cfg .utls .lineOfCode (), 'a[0, 0] ' + str (a [0 , 0 ] ) )
48184837 try :
48194838 # logger
48204839 cfg .utls .logToFile (str (__name__ ) + '-' + str (cfg .inspectSCP .stack ()[0 ][3 ])+ ' ' + cfg .utls .lineOfCode (), 'a ' + str (a [0 ,0 ]) )
4840+ cfg .utls .logToFile (str (__name__ ) + '-' + str (cfg .inspectSCP .stack ()[0 ][3 ])+ ' ' + cfg .utls .lineOfCode (), 'a.dtype ' + str (a .dtype ) )
48214841 except :
48224842 pass
48234843 # logger
@@ -4851,14 +4871,21 @@ def processRasterDev(self, raster = None, signatureList = None, functionBand = N
48514871 except :
48524872 ndvBand = None
48534873 # logger
4874+ cfg .utls .logToFile (str (__name__ ) + '-' + str (cfg .inspectSCP .stack ()[0 ][3 ])+ ' ' + cfg .utls .lineOfCode (), 'ndvBand ' + str (ndvBand ))
4875+ if ndvBand is not None :
4876+ # adapt NoData to dtype
4877+ ndvBand = cfg .np .asarray (ndvBand ).astype (a .dtype ).astype (calcDataType )
4878+ # logger
48544879 cfg .utls .logToFile (str (__name__ ) + '-' + str (cfg .inspectSCP .stack ()[0 ][3 ])+ ' ' + cfg .utls .lineOfCode (), 'ndvBand ' + str (ndvBand ) + ' sclB ' + str (sclB )+ ' offsB ' + str (offsB ))
48554880 if a is not None :
48564881 if functionBandArgument == cfg .multiAddFactorsVar :
48574882 multiAdd = functionVariable
48584883 a = cfg .utls .arrayMultiplicativeAdditiveFactors (a , multiAdd [0 ][b ], multiAdd [1 ][b ])
48594884 array [::, ::, b ] = a .astype (calcDataType )
48604885 # logger
4886+ cfg .utls .logToFile (str (__name__ ) + '-' + str (cfg .inspectSCP .stack ()[0 ][3 ])+ ' ' + cfg .utls .lineOfCode (), 'a[0, 0] ' + str (a [0 , 0 ] ) )
48614887 cfg .utls .logToFile (str (__name__ ) + '-' + str (cfg .inspectSCP .stack ()[0 ][3 ])+ ' ' + cfg .utls .lineOfCode (), 'array[0, 0, b] ' + str (array [0 , 0 , b ] ) )
4888+ cfg .utls .logToFile (str (__name__ ) + '-' + str (cfg .inspectSCP .stack ()[0 ][3 ])+ ' ' + cfg .utls .lineOfCode (), 'ndvBand ' + str (ndvBand .astype (calcDataType )) )
48624889 else :
48634890 procError = 'Error array none'
48644891 # set nodata value
@@ -4873,37 +4900,22 @@ def processRasterDev(self, raster = None, signatureList = None, functionBand = N
48734900 array [::, ::, b ][cfg .np .isnan (array [::, ::, b ])] = ndvBand
48744901 except :
48754902 pass
4876- '''
4877- if ndv is not None:
4878- try:
4879- array[::, ::, b][array[::, ::, b] == ndv] = cfg.np.nan
4880- except:
4881- pass
4882- '''
4883- else :
4903+ try :
4904+ nodataMask [0 , 0 ]
4905+ except :
4906+ nodataMask = cfg .np .zeros ((bSY , bSX ), dtype = calcDataType )
4907+ if skipReplaceNoData is not None :
4908+ pass
4909+ elif ndvBand is not None :
48844910 try :
4885- nodataMask [0 , 0 ]
4911+ nodataMask [::, ::][ array [::, ::, b ] == ndvBand ] = outputNoData
48864912 except :
4887- nodataMask = cfg .np .zeros ((bSY , bSX ), dtype = calcDataType )
4888- if skipReplaceNoData is not None :
48894913 pass
4890- elif ndvBand is not None :
4891- try :
4892- nodataMask [::, ::][a [::, ::] == ndvBand ] = outputNoData
4893- except :
4894- pass
4895- '''
4896- if ndv is not None:
4897- try:
4898- nodataMask[::, ::][a[::, ::, b] == ndv] = outputNoData
4899- except:
4900- pass
4901- '''
4902- # logger
4903- cfg .utls .logToFile (str (__name__ ) + '-' + str (cfg .inspectSCP .stack ()[0 ][3 ])+ ' ' + cfg .utls .lineOfCode (), 'ndvBand ' + str (ndvBand ))
4904- cfg .utls .logToFile (str (__name__ ) + '-' + str (cfg .inspectSCP .stack ()[0 ][3 ])+ ' ' + cfg .utls .lineOfCode (), 'outputNoData ' + str (outputNoData ))
4905- cfg .utls .logToFile (str (__name__ ) + '-' + str (cfg .inspectSCP .stack ()[0 ][3 ])+ ' ' + cfg .utls .lineOfCode (), 'a[0,0] ' + str (a [0 , 0 ]))
4906- cfg .utls .logToFile (str (__name__ ) + '-' + str (cfg .inspectSCP .stack ()[0 ][3 ])+ ' ' + cfg .utls .lineOfCode (), 'nodataMask[0, 0] ' + str (nodataMask [0 , 0 ]))
4914+ # logger
4915+ cfg .utls .logToFile (str (__name__ ) + '-' + str (cfg .inspectSCP .stack ()[0 ][3 ])+ ' ' + cfg .utls .lineOfCode (), 'outputNoData ' + str (outputNoData ))
4916+ cfg .utls .logToFile (str (__name__ ) + '-' + str (cfg .inspectSCP .stack ()[0 ][3 ])+ ' ' + cfg .utls .lineOfCode (), 'a[0,0] ' + str (a [0 , 0 ]))
4917+ cfg .utls .logToFile (str (__name__ ) + '-' + str (cfg .inspectSCP .stack ()[0 ][3 ])+ ' ' + cfg .utls .lineOfCode (), 'nodataMask[0, 0] ' + str (nodataMask ))
4918+ cfg .utls .logToFile (str (__name__ ) + '-' + str (cfg .inspectSCP .stack ()[0 ][3 ])+ ' ' + cfg .utls .lineOfCode (), 'array ' + str (array [::, ::, 0 ]))
49074919 # logger
49084920 cfg .utls .logToFile (str (__name__ ) + '-' + str (cfg .inspectSCP .stack ()[0 ][3 ])+ ' ' + cfg .utls .lineOfCode (), 'skipReplaceNoData ' + str (skipReplaceNoData ))
49094921 #if functionBand is not None and functionBand != 'No':
@@ -4984,7 +4996,7 @@ def processRasterDev(self, raster = None, signatureList = None, functionBand = N
49844996 writeOut = cfg .utls .writeRasterNew (wrtFile , x - roX , y - roY , bSX , bSY , oo , outputNoData , scl , offs , dataType )
49854997 oR = cfg .gdalSCP .Open (wrtFile , cfg .gdalSCP .GA_ReadOnly )
49864998 bO = oR .GetRasterBand (1 )
4987- oo2 = cfg .utls .readArrayBlock (bO , x - roX , y - roY , bSX , bSY )
4999+ oo2 = cfg .utls .readArrayBlock (bO , x - roX , y - roY , bSX , bSY , calcDataType )
49885000 bO = None
49895001 oR = None
49905002 checkOO = cfg .np .allclose (oo , oo2 , equal_nan = True )
0 commit comments