import os

xllcorner_latlon = 2
yllcorner_latlon = 42
cellsize__latlon = 1/120.0
print cellsize__latlon

xll_coord_latlon = xllcorner_latlon
yll_coord_latlon = xllcorner_latlon

scol_latlon = 1 			# starting column
srow_latlon = 1                         # starting row
ncol_latlon = 1440                      # last column
nrow_latlon = 1440                      # last row
#Note: The origin (lower left) is irow=1,icol=1 (yb2t, like in PCRaster).

#print "Starting a new calculation will delete SUM_CELLAREA.ASC?"
#print "Do you want to start a new calculation?"
#os.system('pause')
#os.system('copy SUM_CELLAREA.ASC backup\*.ASC')
os.system('del SUM_CELLAREA.ASC')
area_summary = open('SUM_CELLAREA.ASC','w')
area_summary.write("cell geometry properties \n")
area_summary.write("14 \n")
area_summary.write("id \n")
area_summary.write("lon \n")
area_summary.write("lat \n")
area_summary.write("Ax(m) \n")
area_summary.write("Ay(m) \n")
area_summary.write("Bx(m) \n")
area_summary.write("By(m) \n")
area_summary.write("Cx(m) \n")
area_summary.write("Cy(m) \n")
area_summary.write("Dx(m) \n")
area_summary.write("Dy(m) \n")
area_summary.write("width(m) \n")
area_summary.write("length(m) \n")
area_summary.write("area(m2) \n")

i = 0
for irow_latlon in range(srow_latlon,nrow_latlon+1):
	for icol_latlon in range(scol_latlon,ncol_latlon+1):
		i = i + 1		# id for the coordinate
		Axll_latlon  = xllcorner_latlon + (icol_latlon-1)*cellsize__latlon
		Ayll_latlon  = yllcorner_latlon + (irow_latlon-1)*cellsize__latlon
		Bxll_latlon  = Axll_latlon + cellsize__latlon
		Byll_latlon  = Ayll_latlon
		Cxll_latlon  = Axll_latlon + cellsize__latlon
		Cyll_latlon  = Ayll_latlon + cellsize__latlon
		Dxll_latlon  = Axll_latlon
		Dyll_latlon  = Ayll_latlon + cellsize__latlon
                os.system('del inpfilecs2cs.txt')
                #os.system('del inpfilecs2cs.txt')
                os.system('del outfilecs2cs.txt')
                #os.system('del outfilecs2cs.txt')
                inpfilecs2cs  = open("inpfilecs2cs.txt",'w')
                A_latlon_text = "%20.15f %20.15f %8f \n" % (Axll_latlon, Ayll_latlon, i)
                B_latlon_text = "%20.15f %20.15f %8f \n" % (Bxll_latlon, Byll_latlon, i)
                C_latlon_text = "%20.15f %20.15f %8f \n" % (Cxll_latlon, Cyll_latlon, i)
                D_latlon_text = "%20.15f %20.15f %8f \n" % (Dxll_latlon, Dyll_latlon, i)
                inpfilecs2cs.write(str(A_latlon_text))
                inpfilecs2cs.write(str(B_latlon_text))
                inpfilecs2cs.write(str(C_latlon_text))
                inpfilecs2cs.write(str(D_latlon_text))
                inpfilecs2cs.close()
                os.system('cs2cs +proj=latlong +a=6370997 +f=0 +to +proj=laea +a=6370997 +f=0 +lat_0=55 +lon_0=20 inpfilecs2cs.txt >> outfilecs2cs.txt')
                outfilecs2cs  = open("outfilecs2cs.txt",'r')
                A_LAEAEU_text = outfilecs2cs.readline()
                B_LAEAEU_text = outfilecs2cs.readline()
                C_LAEAEU_text = outfilecs2cs.readline()
                D_LAEAEU_text = outfilecs2cs.readline()
                outfilecs2cs.close()
                A_LAEAEU_coor = A_LAEAEU_text.split()
                B_LAEAEU_coor = B_LAEAEU_text.split()
                C_LAEAEU_coor = C_LAEAEU_text.split()
                D_LAEAEU_coor = D_LAEAEU_text.split()
                #coordinates A,B,C,D in LAEA EU system:
                A_LAEAEU_xll  = float(A_LAEAEU_coor[0])
                A_LAEAEU_yll  = float(A_LAEAEU_coor[1])
                B_LAEAEU_xll  = float(B_LAEAEU_coor[0])
                B_LAEAEU_yll  = float(B_LAEAEU_coor[1])
                C_LAEAEU_xll  = float(C_LAEAEU_coor[0])
                C_LAEAEU_yll  = float(C_LAEAEU_coor[1])
                D_LAEAEU_xll  = float(D_LAEAEU_coor[0])
                D_LAEAEU_yll  = float(D_LAEAEU_coor[1])
                #calculate length, width, and area:
                width1  = (((B_LAEAEU_xll-A_LAEAEU_xll)**(2))+((B_LAEAEU_yll-A_LAEAEU_yll)**(2)))**(0.5)
                width2  = (((C_LAEAEU_xll-D_LAEAEU_xll)**(2))+((C_LAEAEU_yll-D_LAEAEU_yll)**(2)))**(0.5)
                width   = 0.5*(width1+width2)
                length1 = (((C_LAEAEU_xll-B_LAEAEU_xll)**(2))+((C_LAEAEU_yll-B_LAEAEU_yll)**(2)))**(0.5)
                length2 = (((D_LAEAEU_xll-A_LAEAEU_xll)**(2))+((D_LAEAEU_yll-A_LAEAEU_yll)**(2)))**(0.5)
                length  = 0.5*(length1+length2)
                area    = width*length
                #writing the result in the file:                       #%12f #%12.8f      #%12.8f      #%15.3f       #%15.3f       #%15.3f       #%15.3f       #%15.3f       #%15.3f       #%15.3f       #%15.3f     #%15.3f #%15.3f #%18.6f
                outline = "%12i %12.8f %12.8f %15.3f %15.3f %15.3f %15.3f %15.3f %15.3f %15.3f %15.3f %15.3f %15.3f %18.6f \n" % (i, Axll_latlon, Ayll_latlon, A_LAEAEU_xll, A_LAEAEU_yll, B_LAEAEU_xll, B_LAEAEU_yll, C_LAEAEU_xll, C_LAEAEU_yll, D_LAEAEU_xll, D_LAEAEU_yll, width, length, area)
                area_summary.write(str(outline))
                checkpunt = i%10
                if checkpunt > 0:
                	print ""
                	print irow_latlon, icol_latlon
                else:
                	print ""
                	area_summary.close()
                	area_summary = open('SUM_CELLAREA.ASC','a')
                	print irow_latlon, icol_latlon
                if i == 10:
                	#os.system('pause')
area_summary.close()