Changeset 30

Timestamp:

06/19/14 18:22:12 (10 years ago)

Author:

lahlod

Message:

scrit nouvelle grille

Location:

trunk/src/scripts_Laura

Files:

• cartesian_grid_test.py (modified) (1 diff)
• read_emis_glace.py (modified) (2 diffs)

trunk/src/scripts_Laura/cartesian_grid_test.py

@@ -11 +11 @@
 
 
-##############################
-# parameters of the new grid #
-##############################
 
-# apres avoir fait un ffgrid ==> lon et lat regulieres #
+###############
+# test values #
+###############
+longi = array([-109.208, -117.39 , -117.928, -123.927, -120.325, -128.82 , -125.947, -122.596, -133.173, -124.778, -137.312, -134.984, -132.538, -129.883, -126.898, -136.826, -134.414, -131.833, -128.978, -141.05 , -138.707, -136.313, -133.787, -131.036, -127.931, -140.634, -138.243, -135.757, -133.087, -130.122, -126.704, -110.352, -142.613, -140.213, -137.751, -135.144, -132.291, -129.058, -125.246, -114.38 , -147.118, -144.65 , -142.229, -139.777, -137.216, -134.453, -131.372, -127.803, -118.005, -146.753])
 
-Rt = 6371
-alpha = (pi*Rt)/180
-beta = pi/180
-# definition of theta angle and coord x,y #
-x = np.zeros([len(lat_mesh), len(lon_mesh)], float)
-y = np.zeros([len(lat_mesh), len(lon_mesh)], float)
-for ilat in range (0, len(lat_mesh)):
-    for ilon in range (0, len(lon_mesh)):
-        if ((lon_mesh[ilon] >= 0.) & (lon_mesh[ilon] <= 90.)): # 4eme quadrant
-            theta = - (90 - lon_mesh[ilon]) * beta
-            x[ilat, ilon] = (90 - lat_mesh[ilat]) * alpha * cos(theta)
-            y[ilat, ilon] = (90 - lat_mesh[ilat]) * alpha * sin(theta)
-        else:
-            if ((lon_mesh[ilon] > 90.) & (lon_mesh[ilon] <= 180.)): # 1er quadrant
-                theta = (lon_mesh[ilon] - 90) * beta
-                x[ilat, ilon] = (90 - lat_mesh[ilat]) * alpha * cos(theta)
-                y[ilat, ilon] = (90 - lat_mesh[ilat]) * alpha * sin(theta)
-            else: 
-                if ((lon_mesh[ilon] >= -180.) & (lon_mesh[ilon] < 0.)): # 2eme et 3eme quadrants
-                    theta = (270 + lon_mesh[ilon]) * beta
-                    x[ilat, ilon] = (90 - lat_mesh[ilat]) * alpha * cos(theta)
-                    y[ilat, ilon] = (90 - lat_mesh[ilat]) * alpha * sin(theta)
+lati = array([ 68.164,  69.205,  70.121,  70.087,  70.422,  69.994,  70.359, 70.709,  69.846,  70.983,  69.619,  70.061,  70.474,  70.867, 71.246,  70.259,  70.693,  71.103,  71.496,  69.936,  70.437, 70.896,  71.325,  71.735,  72.129,  70.593,  71.08 ,  71.532, 71.959,  72.369,  72.764,  73.708,  70.723,  71.244,  71.721, 72.168,  72.595,  73.007,  73.403,  74.075,  70.188,  70.824, 71.384,  71.891,  72.361,  72.807,  73.235,  73.648,  74.394, 70.892])
+
+z = array([ 0.938 ,  0.899 ,  0.958 ,  0.944 ,  0.991 ,  0.867 ,  0.948 , 0.963 ,  0.911 ,  0.978 ,  0.935 ,  0.932 ,  0.957 ,  0.972 , 0.983 ,  0.945 ,  0.941 ,  0.965 ,  0.986 ,  0.88  ,  0.942 , 0.945 ,  0.966 ,  0.989 ,  0.9999,  0.914 ,  0.908 ,  0.93  , 0.96  ,  0.996 ,  0.9999,  0.952 ,  0.898 ,  0.88  ,  0.852 , 0.86  ,  0.906 ,  0.984 ,  0.976 ,  0.888 ,  0.794 ,  0.89  , 0.846 ,  0.849 ,  0.869 ,  0.812 ,  0.764 ,  0.904 ,  0.853 , 0.934 ])
+z0 = min(z)
+z1 = max(z)
 
 
-# definition of nex grid cells # 
+####################################################
+# from polar (lon, lat) to cartesian coords (x, y) #
+####################################################
+# associates a (x, y) couple to each point of (longitude, latitude) coordinates 
+# origin of the (x, y) grid : (x=0, y=0) <=> (lon=0, lat=0) 
+Rt = 6371.
+alpha = (pi*Rt)/180.
+beta = pi/180.
+# definition of theta angle and coord x,y #
+x = np.zeros([len(lati), len(longi)], float)
+y = np.zeros([len(lati), len(longi)], float)
+for ilat in range (0, len(lati)):
+    for ilon in range (0, len(longi)):
+        if ((longi[ilon] >= 0.) & (longi[ilon] <= 90.)): # 4eme quadrant
+            theta = (90. - longi[ilon]) * beta
+            x[ilat, ilon] = (90. - lati[ilat]) * alpha * cos(theta)
+            y[ilat, ilon] = (90. - lati[ilat]) * alpha * sin(-theta)
+        else:
+            if ((longi[ilon] > 90.) & (longi[ilon] <= 180.)): # 1er quadrant
+                theta = (longi[ilon] - 90.) * beta
+                x[ilat, ilon] = (90. - lati[ilat]) * alpha * cos(theta)
+                y[ilat, ilon] = (90. - lati[ilat]) * alpha * sin(theta)
+            else: 
+                if ((longi[ilon] >= -180.) & (longi[ilon] < 0.)): # 2eme et 3eme quadrants
+                    theta = (270. + longi[ilon]) * beta
+                    x[ilat, ilon] = (90. - lati[ilat]) * alpha * cos(theta)
+                    y[ilat, ilon] = (90. - lati[ilat]) * alpha * sin(theta)
+
+
+# definition of the new cartesian grid (x, y) #
+Sx = size(x)
+Sy = size(y)
 xx = reshape(x, size(x))
-yy = np.append(y[:,180], y[:,360])
+yy = reshape(y, size(y))
+mx = min(xx)
+Mx = max(xx)
+my = min(yy)
+My = max(yy)
 
-x0 = -2780.
-x1 = 2780.
-y0 = -2780.
-y1 = 2780.
-dx = 2
-dy = 2
-xvec = arange(x0, x1 + dx, dx)
-yvec = arange(y0, y1 + dy, dy)
+x0 = round(mx) - 50. # xmin
+x1 = round(Mx) + 50. # xmax
+dx = 25. # delta(x)
+xvec = np.arange(x0, x1+dx, dx)
 nx = size(xvec)
+y0 = round(my) - 50. # ymin
+y1 = round(My) + 50. # ymax
+dy = 25. # delta(y)
+yvec = np.arange(y0, y1+dy, dy)
 ny = size(yvec)
-n = size(xx)
+xgrid_cart, ygrid_cart = np.meshgrid(xvec, yvec) # new cartesian grid (centered on North pole)
 
-ix = np.zeros([n], int)
-for kk in range(0, n):
-    if xx[kk] == x0:
-        ix[kk] = 0
+
+# counting each grid cell #
+ix = np.zeros([Sx], int)
+for i in range (0, Sx): # boucle sur les points M (abscisses)
+    if xx[i] == x0:
+        ix[i] = 0
     else:
-        ix[kk] = math.ceil((xx[kk] - x0)/dx)
+        ix[i] = math.ceil((xx[i] - x0) / dx) - 1
 
-iy = np.zeros([n], int)
-for kk in range(0, n):
-    if yy[kk] == min(y):
-        iy[kk] = 0
+iy = np.zeros([Sy], int) 
+for j in range (0, Sy): # boucle sur les points M (ordonnees)
+    if yy[j] == y0:
+        iy[j] = 0
     else:
-        iy[kk] = math.ceil((yy[kk] - y0)/dy)
+        iy[j] = math.ceil((yy[j] - y0) / dy) - 1 
 
-zz = reshape(emis_mesh, size(emis_mesh))
-zgrid = np.zeros([len(ixx), len(iyy)], float)
-for i in range (0, size(x)):
-    zgrid[ixx[i], iyy[i]] = zgrid[ixx[i], iyy[i]] + zz[i]
+
+#inx = (ix >= 0) & (ix < nx)
+#iny = (iy >= 0) & (iy < ny)
+#inz = (z > z0) & (z < z1)
+#inn = inx & iny & inz
+#iix = ix[inn]
+#iiy = iy[inn]
+#zz = z[inn]
 
 
 
-xii,yii = np.meshgrid(xx ,yy)
+
+# calculation of distances between point M(x,y) and 4 grid points of its belonging cell #
+dist1 = np.zeros([Sx], float)
+dist2 = np.zeros([Sx], float)
+dist3 = np.zeros([Sx], float)
+dist4 = np.zeros([Sx], float)
+min_dist = np.zeros([Sx], float)
+for k in range (0, Sx): # boucle sur les points M (x et y)
+    dist1[k] = sqrt(((xx[k] - xvec[ix[k]]) ** 2) + ((yy[k] - yvec[iy[k]]) ** 2))
+    dist2[k] = sqrt(((xx[k] - xvec[ix[k] + 1]) ** 2) + ((yy[k] - yvec[iy[k]]) ** 2))
+    dist3[k] = sqrt(((xx[k] - xvec[ix[k]])**2) + ((yy[k] - yvec[iy[k] + 1]) ** 2))
+    dist4[k] = sqrt(((xx[k] - xvec[ix[k] + 1]) ** 2) + ((yy[k] - yvec[iy[k] + 1]) ** 2))
+    dist_vec = np.array([dist1[k], dist2[k], dist3[k], dist4[k]])
+    ind = np.where(dist_vec == min(dist_vec))
+    print 'grid cell no : ' + str(ind[0][0]+1)
+    
+
+##### probleme: z pas la même dimension (plus petite)!!#######
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+for i in range (0, 5):
+    np.where((xx - xvec[i]) <= dx)
+        
+    for j in range (0, len(yy)):
+        np.where ((xx[i] - xvec) <= dx):
+        print 'ok'
+        else:
+            print 'not ok'
+        np.where((yy[i] - yvec) <= dy)
+
+
+
+M = np.array([xx, yy])
+for k in range (0, len(xx)):
+    xcoord[k] = M[:,k][0]
+    ycoord[k] = M[:,k][1]
+
+
+
+
+ix = np.zeros([len(xx)], int)
+for i in range (0, len(xx)):
+    if xx[i] == x0:
+        ix[i] = 0
+    else:
+        ix[i] = math.ceil((xx[i] - x0) / dx) - 1
+
+iy = np.zeros([len(yy)], int)
+for j in range (0, len(yy)):
+    if yy[j] == y0:
+        iy[j] = 0
+    else:
+        iy[j] = math.ceil((yy[j] - y0) / dy) - 1 
+
+
+
+
+
+
 
 figure()

trunk/src/scripts_Laura/read_emis_glace.py

@@ -35 +35 @@
 tup=np.zeros([4,nbtotal],float)
 tdn=np.zeros([4,nbtotal],float)
+#trans=np.zeros([4,nbtotal],float)
 
 
@@ -94 +95 @@
 lat_mesh = outy
 
+plt.ion()
 figure()
 #m = Basemap(llcrnrlon=-180,urcrnrlon=180,llcrnrlat=65,urcrnrlat=90,projection='cyl',resolution='c',fix_aspect=True)
 m = Basemap(projection = 'npaeqd',boundinglat = 60, lon_0 = 0, resolution = 'l')
 xii, yii = m(*np.meshgrid(lon_mesh,lat_mesh))
-clevs=arange(zz0,zz1,0.01)
+#clevs=arange(-2400, -800, 10)
 m.drawcoastlines(linewidth=1)
 m.drawparallels(np.arange(60, 90, 20))
 m.drawmeridians(np.arange(-180, 180, 20))
-cs=m.contourf(xii,yii,emis_mesh, clevs, cmap=cm.s3pcpn_l_r) # emissivity
+cs=m.contourf(xii,yii,y, cmap=cm.s3pcpn_l_r) # emissivity
 cbar =colorbar(cs)
 #cbar.set_label(txt)