I have an array from Euro-Cordex data which has a rotated pole projection from a Netcdf file:
grid_mapping_name: rotated_latitude_longitude
grid_north_pole_latitude: 39.25
grid_north_pole_longitude: -162.0
float64 rlon(rlon)
standard_name: grid_longitude
long_name: longitude in rotated pole grid
units: degrees
axis: X
unlimited dimensions:
current shape = (424,)
filling on, default _FillValue of 9.969209968386869e+36 used),
('rlat', <class 'netCDF4._netCDF4.Variable'>
float64 rlat(rlat)
standard_name: grid_latitude
long_name: latitude in rotated pole grid
units: degrees
axis: Y
unlimited dimensions:
current shape = (412,)
The dimensions are rlon (424) and rlat (412). I used some codes to convert these rotated lat lons into normal lat/lons. Now, I have two matrices with shape of (424, 412). The first one shows the longitude coordinates, and the second one shows the latitude coordinates.
Now, I want to convert the initial image (424, 412) to a image with the extents that I want:Min lon : 25, Max lon: 45, Min Lat: 35, Max lat: 43
lats = np.empty((len(rlat), len(rlon)))
lons = np.empty((len(rlat), len(rlon)))
for j in range (len(rlon)):
for i in range(len(rlat)):
lons[i, j] = unrot_lon(rlat[i],rlon[j],39.25,-162.0)
lats[i, j] = unrot_lat(rlat[i],rlon[j],39.25,-162.0)
a = lons<=45
aa = lons>=25
aaa = a*aa
b = lats<=43
bb = lats>=35
bbb = b*bb
c = bbb*aaa
The last matrix (c) is a boolean matrix which shows the pixels that I am interested according to the extents that I defined:
Now, I want to do two things that I fail in both:
First I would like to plot this image with the boundries on a basemap. For that I located the llcrnlon, llcrnlat, urcrnlon and urcrnlon based on the boolean matrix and by using some imagination:
llcrlon = 25.02#ok
llcrlat = np.nanmin(lats[c])# ok
urcrlon = np.nanmax(lons[c])#ok
urcrlat = np.nanmax(lats[np.where(lons==urcrlon)])#ok
Then I used the following codes to plot the image on a basemap:
lonss = np.linspace(np.min(lons[c]), np.max(lons[c]), (424-306+1))
latss = np.linspace(np.min(lats[c]), np.max(lats[c]), (170-73+1))
pl.figure(dpi = 250)
map = Basemap(projection='rotpole',llcrnrlon=llcrlon,llcrnrlat=llcrlat,urcrnrlon=urcrlon,urcrnrlat=urcrlat,resolution='i', o_lat_p = 39.25, o_lon_p =-162., lon_0=35, lat_0=45)
map.drawcoastlines()
map.drawstates()
parallels = np.arange(35,43,2.) #
meridians = np.arange(25,45,2.) #
map.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)
map.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)
lons, lats = np.meshgrid(lonss, latss)
x, y = map(lons, lats)
mapp = map.pcolormesh(x,y,WTD[73:170, 306:])
So, the map is not well-fit to the basemap projection. I would like to find out what is wrong.
Second, I would like to reproject this map to normal lat/lon. For that, I use the following codes to define a new grid:
targ_lons = np.linspace(25, 45, 170)
targ_lats = np.linspace(43, 35, 70)
T_Map = np.empty((len(targ_lats), len(targ_lons)))
T_Map[:] = np.nan
Then, I am trying to figure out the differences between the lon/lat matrices I produced in the beginning and my newly defined grids. Then, using the indices which represent the minimum/less than a specific threshold, fill in the new gridded image.
for i in range(len(targ_lons)):
for j in range(len(targ_lats)):
lon_extr = np.where(abs(lons-targ_lons[i])<0.01)
lat_extr = np.where(abs(lats-targ_lats[j])<0.01)
So here, if we have i=0 and j=0,
then:
lon_extr = (array([ 7, 16, 25, 34, 35, 43, 44, 53, 63, 72, 73, 82, 83, 92, 93, 102, 103, 112, 113, 122, 123, 133, 143, 153, 154, 164,
174, 175, 185, 195, 196, 206, 217, 227, 238, 248, 259, 269, 280,
290, 300, 321, 331, 341, 360, 370, 389], dtype=int64),
array([320, 319, 318, 317, 317, 316, 316, 315, 314, 313, 313, 312, 312,
311, 311, 310, 310, 309, 309, 308, 308, 307, 306, 305, 305, 304,
303, 303, 302, 301, 301, 300, 299, 298, 297, 296, 295, 294, 293,
292, 291, 289, 288, 287, 285, 284, 282], dtype=int64))
and
lat_extr=(array([143, 143, 143, 143, 143, 143, 143, 143, 143, 143, 143, 143, 143,
143, 143, 143, 143, 144, 144, 144, 144, 144, 144, 145, 145, 145,
145, 146, 146, 146, 146, 147, 147, 147, 148, 148, 149, 149, 150,
150, 151, 151, 152, 152, 153, 153, 154, 154, 155, 156, 156, 157,
157, 158, 158, 159, 159, 160, 160, 161, 162, 162, 163, 164, 164,
165, 167, 168, 168, 169, 169, 170, 170, 171, 174, 175, 177, 178,
180, 181, 183, 186, 190, 191, 192, 204, 205, 210, 214], dtype=int64),
array([251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 263,
264, 265, 266, 267, 227, 228, 229, 289, 290, 291, 214, 215, 303,
304, 204, 205, 313, 314, 196, 321, 322, 189, 329, 182, 336, 176,
342, 170, 348, 165, 353, 160, 358, 155, 363, 150, 146, 372, 142,
376, 138, 380, 134, 384, 130, 388, 126, 123, 395, 119, 116, 402,
405, 106, 103, 415, 100, 418, 97, 421, 94, 86, 83, 78, 75,
70, 68, 63, 56, 47, 45, 43, 19, 17, 8, 1], dtype=int64))
Now, I need to be able to pull the common coordinates and fill in the T_Map. I'm confused at this point. Is there a function for easy way to pull out the common lat/lon from these two arrays?
