Demo entry 6338987

computeCN

   

Submitted by anonymous on Dec 18, 2016 at 19:12
Language: Python. Code size: 3.1 kB.

from netCDF4 import Dataset
import datetime
import numpy
from PIL import Image
Image.MAX_IMAGE_PIXELS = None
c2c_dict = {11:0,14:1,20:7,30:8,40:45,50:45,60:44,70:44,90:43,100:43,110:42,120:40,130:38,140:36,150:30,160:0,170:0,180:0,190:46,200:48,210:0,220:0}
lc_value = c2c_dict.keys()
CNbase = [77,86,91,94,76,85,90,93,74,83,88,90,72,81,88,91,67,78,85,89,71,80,87,90,64,75,82,85,70,79,84,88,65,75,82,86,69,78,83,87,64,74,81,85,66,74,80,82,62,71,78,81,65,73,79,81,61,70,77,80,65,76,84,88,63,75,83,87,64,75,83,86,60,72,80,84,63,74,82,85,61,73,81,84,62,73,81,84,60,72,80,83,61,72,79,82,59,70,78,81,60,71,78,81,58,69,77,80,66,77,85,89,58,72,81,85,64,75,83,85,55,69,78,83,63,73,80,83,51,67,76,80,68,79,86,89,49,69,79,84,39,61,74,80,30,58,71,78,48,67,77,83,35,56,70,77,30,48,65,73,57,73,82,86,43,65,76,82,32,58,72,79,45,66,77,83,36,60,73,79,30,55,70,77,59,74,82,86,72,82,87,89,76,85,89,91,-9999,-9999,-9999,-9999,-9999,-9999,-9999,-9999,-9999,-9999,-9999]

#read nc: K_s and CN lookup table
#read tiff: landcover
file1 = 'HSG.nc'
file2 = 'LookupCN.nc'
file3 = 'GLOBCOVER_L4_200901_200912_V2.3_China.tif'

nc_K = Dataset(file1,'r+',diskless=False)
print 'nc read OK'
#nc_lookup = Dataset(file2,'r')
tiff_lc = Image.open(file3)
a_lc = numpy.array(tiff_lc)
print 'TIFF read OK'

dim_lat = len(nc_K.dimensions['lat'])
dim_lon = len(nc_K.dimensions['lon'])
HSG = nc_K.variables['HSG']

#create new nc variables to save CN results
cn = nc_K.createVariable('CN_7k','i',('lat','lon'))

#save mismatch[HSGLat,HSGLon], valid in HSG but not in lc
#due to lc position can be computed by HSG position
#only HSG position is saved
mm = []

#traverse K_s to determine CN
def determincn(argi,argj):
    soilCode = HSG[argi,argj]

    lcPosLat = 4319-argi
    lcPosLon = argj
    lcCode = a_lc[lcPosLat,lcPosLon]
    #print argi,',',argj,':',lcCode
    if lcCode in lc_value:
        c2cCode = c2c_dict[lcCode]
        rv = CNbase[c2cCode*4+soilCode]
    else:
        rv = -9999
        mm.append((argi,argj))
        # print 'HSG:',argi,',',argj,' LC: ',lcPosLat,',',lcPosLon
    
    return rv

##def deleteNCvariable(inFile, outFile):
##    inNC = Dataset(infile, 'r')
##    outNC = Dataset(outfile, 'w', format = 'NETCDF4_CLASSIC')
##
##    dim1 = inNC.dimensions['lat']
##    dim2 = inNC.dimensions['lon']
##    var1 = inNC.variables['HSG']
##    var2 = inNC.variables['lat']
##    var3 = inNC.variables['lon']
##
##    outNC.createDimension('lat')
##    outNC.createDimension('lon')
##
##    latitude = outNC.createVariable('lat','d',('lat','lon'))
##    longitude = outNC.createVariable('lon','d',('lat','lon'))
##    HSG = outNC.createVariable('lon','h',('lat','lon'))
##    latitude = var2
##    longitude = var3
##    HSG = var1
##    outNC.close()
##
##    return 1
    

for i in range(0,dim_lat):
    print i
    for j in range(0,dim_lon):
        if HSG[i,j]==-9999:
            cn[i,j] = -9999
        else:
            cn[i,j] = determincn(i,j)
        #print 'CN[',i,',',j,']=',cn[i,j]

#close CN
nc_K.close()
f1 = file('pyOutMismatch.txt','w')
f1.write(str(mm))

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).