Changeset 31

Timestamp:

06/20/14 18:55:02 (10 years ago)

Author:

lahlod

Message:

new_cartesian_grid and read_bathy

Location:

trunk/src/scripts_Laura

Files:

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

trunk/src/scripts_Laura/cartesian_grid_test.py

@@ -10 +10 @@
 import draw_map
 
-
-
-###############
-# 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])
-
-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)
-
-
 ####################################################
 # 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)
+
+def new_cartesian_grid(longitude, latitude, z, z0, z1, dx, dy):
+
+    # 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) 
+    L = len(z)
+    #z0 = min(z)
+    #z1 = max(z)
+    Rt = 6371.
+    alpha = (pi*Rt)/180.
+    beta = pi/180.
+    x = np.zeros([L], float)
+    y = np.zeros([L], float)
+    for k in range (0, L):
+        if ((longitude[k] >= 0.) & (longitude[k] <= 90.)): # 4eme quadrant
+            theta = (90. - longitude[k]) * beta
+            x[k] = (90. - latitude[k]) * alpha * cos(theta)
+            y[k] = (90. - latitude[k]) * 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)
+            if ((longitude[k] > 90.) & (longitude[k] <= 180.)): # 1er quadrant
+                theta = (longitude[k] - 90.) * beta
+                x[k] = (90. - latitude[k]) * alpha * cos(theta)
+                y[k] = (90. - latitude[k]) * 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)
+                if ((longitude[k] >= -180.) & (longitude[k] < 0.)): # 2eme et 3eme quadrants
+                    theta = (270. + longitude[k]) * beta
+                    x[k] = (90. - latitude[k]) * alpha * cos(theta)
+                    y[k] = (90. - latitude[k]) * alpha * sin(theta)
 
 
-# definition of the new cartesian grid (x, y) #
-Sx = size(x)
-Sy = size(y)
-xx = reshape(x, size(x))
-yy = reshape(y, size(y))
-mx = min(xx)
-Mx = max(xx)
-my = min(yy)
-My = max(yy)
-
-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)
-xgrid_cart, ygrid_cart = np.meshgrid(xvec, yvec) # new cartesian grid (centered on North pole)
+    # definition of the new cartesian grid (x, y) #
+    Mx = max(x)
+    mx = min(x)
+    x0 = round(mx) - 50. # xmin
+    x1 = round(Mx) + 50. # xmax
+    #dx = 50. # delta(x)
+    xvec = np.arange(x0, x1+dx, dx)
+    nx = len(xvec)
+    My = max(y)
+    my = min(y)
+    y0 = round(my) - 50. # ymin
+    y1 = round(My) + 50. # ymax
+    #dy = 50. # delta(y)
+    yvec = np.arange(y0, y1+dy, dy)
+    ny = len(yvec)
+    xgrid_cart, ygrid_cart= np.meshgrid(xvec, yvec) # new cartesian grid (centered on North pole)
 
 
-# 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[i] = math.ceil((xx[i] - x0) / dx) - 1
+    # counting each grid cell #
+    ix = np.zeros([L], int)
+    i = 0
+    for i in range (0, L): # boucle sur les points M (abscisses) 
+        if x[i] == x0:
+            ix[i] = 0
+        else:
+            ix[i] = math.ceil((x[i] - x0) / dx) - 1 # associates each x vakue to a cell number
 
-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[j] = math.ceil((yy[j] - y0) / dy) - 1 
+    i = 0
+    iy = np.zeros([L], int) 
+    for i in range (0, L):# boucle sur les points M (ordonnees) 
+        if y[i] == y0:
+            iy[i] = 0
+        else:
+            iy[i] = math.ceil((y[i] - y0) / dy) - 1 # associates each y vakue to a cell number
 
 
-#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]
+    # calculation of distances between point M(x,y) and 4 grid points of its belonging cell #
+    close_point = np.zeros([L, 2], int)
+    k = 0
+    for k in range (0, L): # boucle sur les points M (x et y)
+    #    print 'x', x[k]
+    #    print 'y', y[k]
+    #    print 'xvec', xvec[ix[k]]
+    #    print 'yvec', yvec[iy[k]]
+        d1 = sqrt(((x[k] - xvec[ix[k]]) ** 2) + ((y[k] - yvec[iy[k]]) ** 2))
+     #   print 'd1', d1
+        d2 = sqrt(((x[k] - xvec[ix[k] + 1]) ** 2) + ((y[k] - yvec[iy[k]]) ** 2))
+    #    print 'd2', d2
+        d3 = sqrt(((x[k] - xvec[ix[k]]) ** 2) + ((y[k] - yvec[iy[k] + 1]) ** 2))
+     #   print 'd3', d3
+        d4 = sqrt(((x[k] - xvec[ix[k] + 1]) ** 2) + ((y[k] - yvec[iy[k] + 1]) ** 2))
+    #    print 'd4', d4
+        d_vec = np.array([d1, d2, d3, d4])
+        ind = np.where(d_vec == min(d_vec)) # finds the smallest distance between the 4 points of the grid
+    #    print 'grid cell no : ' + str(ind[0][0]+1)
+        point_vec = np.array([(ix[k], iy[k]), (ix[k] + 1, iy[k]), (ix[k], iy[k] + 1), (ix[k] + 1, iy[k] + 1)])
+    #    print 'chosen point in grid', point_vec[ind[0][0]] # we have chosen which point of the grid is closer to M // point_vec[ind[0][0]] = (cell no of closest x, celle no of closest y)
+        close_point[k, :] = point_vec[ind[0][0]]# we have chosen which point of the grid is closer to M // point_vec[ind[0][0]] = (cell no of closest x, celle no of closest y)
 
 
-
-
-# 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)
+    # association of z value to closest grid point #
+    zgrid = np.zeros([ny, nx], float) # z in new grid
+    ngrid = np.zeros([ny, nx], int) # nb of obs per new grid cell
+    z2grid = np.zeros([ny, nx], float) # z2 in new grid
+    for k in range (0, L):
+        zgrid[close_point[k, 1], close_point[k, 0]] = zgrid[close_point[k, 1], close_point[k, 0]] + z[k]
+        ngrid[close_point[k, 1], close_point[k, 0]] = ngrid[close_point[k, 1], close_point[k, 0]] + 1
+        z2grid[close_point[k, 1], close_point[k, 0]] = z2grid[close_point[k, 1], close_point[k, 0]] + (z[k] * z[k])
     
 
-##### probleme: z pas la même dimension (plus petite)!!#######
+    # z in each point of new grid #
+    ZGRID = zgrid / ngrid
+    # variance of z in each grid cell #
+    sigmagrid = np.zeros([ny, nx], float)
+    for j in range (0, nx):
+        for i in range (0, ny):
+            if (ngrid[i, j] > 1): # take away the cells where no obs 
+                sigmagrid[i, j] = (1 / (ngrid[i, j] - 1) * sqrt(z2grid[i, j] - ngrid[i, j] * ZGRID[i, j] * ZGRID[i, j]))
+            else:
+                sigmagrid[i, j] = NaN
+        
+
+    return ZGRID, ngrid, sigmagrid, xvec, yvec, xgrid_cart, ygrid_cart 
 
 
@@ -127 +133 @@
 
 
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-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()
-m = Basemap(projection = 'npaeqd',boundinglat = 60, lon_0 = 0, resolution = 'l')
-clevs=arange(0.5, 1., 0.01)
-m.drawcoastlines(linewidth=1)
-m.drawparallels(np.arange(60, 90, 20))
-m.drawmeridians(np.arange(-180, 180, 20))
-xii,yii = m(*np.meshgrid(xx ,yy))
-cs=m.contourf(xii,yii, zgrid, clevs, cmap=cm.s3pcpn_l_r)
-cbar =colorbar(cs)
-
+#emisGRID, nemis, sigmaemis, xvec, yvec, xgrid_cart, ygrid_cart = cartesian_grid_test.new_cartesian_grid(lon_zon, lat_zon, emis_zon, emis_zon.min(), emis_zon.max(), 40, 40)

trunk/src/scripts_Laura/read_emis_glace.py

@@ -113 +113 @@
 
 
+import