import numpy as np from datetime import datetime, timedelta import sys __author__ = 'Trond Kristiansen' __email__ = 'trond.kristiansen@imr.no' __created__ = datetime(2010, 1, 7) __modified__ = datetime(2010, 1, 7) __version__ = "0.1" __status__ = "Development, modified on 07.01.2010" def findSubsetIndices(grdMODEL,min_lat,max_lat,min_lon,max_lon): """ Get the indices that covers the new grid, and enables us to only store a subset of the large input grid. This routine checks wether you cross the 0 longitude line, so that start is in the West (negative longitudes) and end point is in the East (positive longitudes). if that is the case, the extraction of data is split at the 0 line so that we first get all data in the Western region before all the data in the Eastern region is taken. Then we combine the two arrays so that West comes before East. This is not the case in an array where the longitude values goes from 0 to 360. In such an array (as SODA is), the Westernmost data comes after the Easternmost. Therefore we need to concentate the Westernomst data into an array with the Easternomst data in the correct order. Input: Input longitudes has to range from -180 to 180 Returns: list of indices, boolean value that is True of splitting occurs (crosses 0 long) Indices are returned so that the Westernomst comes first Trond Kristiansen, 22.07.2009, 08.01.2010 """ if min_lon<0 and max_lon>0: splitExtract = True; Turns=2 grdMODEL.splitExtract=splitExtract else: splitExtract = False; Turns=1 grdMODEL.splitExtract=splitExtract grdMODEL.lon = np.where(grdMODEL.lon>180,grdMODEL.lon-360,grdMODEL.lon) """Array to store the results returned from the function""" res=np.zeros((Turns,4),dtype=np.float64) lats=grdMODEL.lat[:,0] lons=grdMODEL.lon[0,:] for k in range(Turns): if k==0 and splitExtract == True: minLon=min_lon; maxLon=0 minLon=minLon+360 maxLon=maxLon+360 elif k==1 and splitExtract == True: minLon=0; maxLon=max_lon else: minLon=min_lon; maxLon=max_lon distances1 = [] distances2 = [] indices=[] index=1 for point in lats: s1 = max_lat-point # (vector subtract) s2 = min_lat-point # (vector subtract) distances1.append((np.dot(s1, s1), point, index)) distances2.append((np.dot(s2, s2), point, index-1)) index=index+1 distances1.sort() distances2.sort() indices.append(distances1[0]) indices.append(distances2[0]) distances1 = [] distances2 = [] index=1 for point in lons: s1 = maxLon-point # (vector subtract) s2 = minLon-point # (vector subtract) distances1.append((np.dot(s1, s1), point, index)) distances2.append((np.dot(s2, s2), point, index-1)) index=index+1 distances1.sort() distances2.sort() indices.append(distances1[0]) indices.append(distances2[0]) """ Save final product: max_lat_indices,min_lat_indices,max_lon_indices,min_lon_indices""" minJ=indices[1][2] maxJ=indices[0][2] minI=indices[3][2] maxI=indices[2][2] res[k,0]=minI; res[k,1]=maxI; res[k,2]=minJ; res[k,3]=maxJ; """ Save final product: max_lat_indices,min_lat_indices,max_lon_indices,min_lon_indices""" grdMODEL.indices=res def checkDomain(grdMODEL,grdROMS): lonCHECK=False latCHECK=False if (grdMODEL.lon.min() < grdROMS.lon_rho.min() and grdMODEL.lon.max() > grdROMS.lon_rho.max()): lonCHECK=True if (grdMODEL.lat.min() < grdROMS.lat_rho.min() and grdMODEL.lat.max() > grdROMS.lat_rho.max()): latCHECK=True print "\n--------------------------" print "---> Area output: (longitude=%3.2f,latitude=%3.2f) to (longitude=%3.2f,latitude=%3.2f)"%(grdROMS.lon_rho.min(),grdROMS.lat_rho.min(),grdROMS.lon_rho.max(),grdROMS.lat_rho.max()) print "---> Area input : (longitude=%3.2f,latitude=%3.2f) to (longitude=%3.2f,latitude=%3.2f)"%(grdMODEL.lon.min(),grdMODEL.lat.min(),grdMODEL.lon.max(),grdMODEL.lat.max()) if latCHECK is True and lonCHECK is True: print "Domain check passed: Input domain data covers output domain" else: print "WARNING: Your input domain is smaller or not overlaying your output domain" print "IOsubset.py: EXIT" sys.exit() print "\n--------------------------" def organizeSplit(grdMODEL,grdROMS,type,vars,cdf): if grdMODEL.splitExtract is True: print "\nThe subset crosses the 0 degrees longitude line." print "Need to split the extraction of data into Western and" print "Eastern region. The two regions will be concatenated again.\n" """Save longitude and latitude values so that only variables are needed to be extracted after first loop""" lon1=grdMODEL.lon[int(grdMODEL.indices[0,2]):int(grdMODEL.indices[0,3]),int(grdMODEL.indices[0,0]):int(grdMODEL.indices[0,1])] lat1=grdMODEL.lat[int(grdMODEL.indices[0,2]):int(grdMODEL.indices[0,3]),int(grdMODEL.indices[0,0]):int(grdMODEL.indices[0,1])] lon2=grdMODEL.lon[int(grdMODEL.indices[1,2]):int(grdMODEL.indices[1,3]),int(grdMODEL.indices[1,0]):int(grdMODEL.indices[1,1])] lat2=grdMODEL.lat[int(grdMODEL.indices[1,2]):int(grdMODEL.indices[1,3]),int(grdMODEL.indices[1,0]):int(grdMODEL.indices[1,1])] """Convert the positive numbers above 180 to negative values""" lon1=np.where(lon1>180,lon1-360,lon1) grdMODEL.lon = np.concatenate((lon1,lon2),axis=1) grdMODEL.lat = np.concatenate((lat1,lat2),axis=1) print "Subset extracted for domain from West (%s) to East (%s)\n"%(grdMODEL.lon.min(),grdMODEL.lon.max()) else: grdMODEL.lon=grdMODEL.lon[int(grdMODEL.indices[0,2]):int(grdMODEL.indices[0,3]),int(grdMODEL.indices[0,0]):int(grdMODEL.indices[0,1])] grdMODEL.lat=grdMODEL.lat[int(grdMODEL.indices[0,2]):int(grdMODEL.indices[0,3]),int(grdMODEL.indices[0,0]):int(grdMODEL.indices[0,1])] checkDomain(grdMODEL,grdROMS)