Changeset 1173 in lmdz_wrf for trunk/tools
- Timestamp:
- Oct 11, 2016, 4:25:28 PM (9 years ago)
- Location:
- trunk/tools
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/tools/module_ForInterpolate.F90
r900 r1173 1 1 ! Module to interpolate values from a giving projection and pressure interpolation 2 2 ! To be included in a python 3 ! f2py -m module_ForInterpolate --f90exec=/usr/bin/gfortran-4.7 -c module_generic.F90 module_ForInterpolate.F90 >& run_f2py.log 3 ! Follow compilation instructions from Makefile 4 ! Content 5 ! LlInterpolateProjection: Subroutine which provides the indices for a given interpolation of a projection 6 ! var2D_IntProj: Subroutine to interpolate a 2D variable 7 ! var3D_IntProj: Subroutine to interpolate a 3D variable 8 ! var4D_IntProj: Subroutine to interpolate a 4D variable 9 ! var5D_IntProj: Subroutine to interpolate a 5D variable 4 10 MODULE module_ForInterpolate 5 11 … … 714 720 END SUBROUTINE CoarseInterpolateExact 715 721 722 SUBROUTINE LlInterpolateProjection(inlonv, inlatv, projlon, projlat, intkind, outLlw, idimx, idimy, & 723 pdimx, pdimy) 724 ! Subroutine which provides the indices for a given interpolation of a projection 725 726 USE module_generic 727 728 IMPLICIT NONE 729 730 INTEGER, PARAMETER :: r_k = KIND(1.d0) 731 INTEGER, INTENT(in) :: idimx, idimy, pdimx, pdimy 732 REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in) :: projlon, projlat 733 REAL(r_k), DIMENSION(idimx,idimy), INTENT(in) :: inlonv, inlatv 734 CHARACTER(LEN=50), INTENT(in) :: intkind 735 REAL(r_k), DIMENSION(3,16,pdimx,pdimy), INTENT(out) :: outLlw 736 737 ! Local 738 INTEGER :: i,j,iv, ix, iy 739 REAL(r_k) :: mindiffLl, dist 740 REAL(r_k), DIMENSION(idimx,idimy) :: difflonlat 741 REAL(r_k), DIMENSION(2) :: extremelon, extremelat 742 INTEGER, DIMENSION(2) :: iLl 743 CHARACTER(LEN=50) :: fname 744 745 !!!!!!! Variables 746 ! idimx, idimy: dimension length of the input projection 747 ! pdimx, pdimy: dimension length of the target projection 748 ! in[lon/lat]: longitudes and latitudes of the target projection 749 ! proj[lon/lat]: longitudes and latitudes of the target projection 750 ! intkind: kind of interpolation 751 ! 'npp': nearest neighbourgh 752 ! 'dis': weighted distance, grid-output for SW, NW, NE, SE 753 ! outLlw: output interpolation result 754 ! for point pi,pj; up to 16 different values of 755 ! 1st: i-index in input projection 756 ! 2nd: j-index in input projection 757 ! 3rd: weight for i-index, j-index to use for ponderation (0<1.) 758 fname = 'LlInterpolateProjection' 759 760 extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /) 761 extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /) 762 763 ! Using case outside loop to be more efficient 764 SELECT CASE(TRIM(intkind)) 765 CASE ('dis') 766 DO i=1, pdimx 767 DO j=1, pdimy 768 ! Find point 769 difflonlat = SQRT((projlon(i,j)-inlonv)**2. + (projlat(i,j)-inlatv)**2.) 770 mindiffLl = MINVAL(difflonlat) 771 ! Getting the four surrounding values 772 iLl = index2DArrayR(difflonlat, idimx, idimy, mindiffLl) 773 IF ( projlon(i,j) < inlonv(iLl(1),iLl(2)) ) THEN 774 outLlw(1,1,i,j) = MAX(iLl(1)-1,1) 775 outLlw(1,2,i,j) = MAX(iLl(1)-1,1) 776 outLlw(1,3,i,j) = MIN(iLl(1),idimx) 777 outLlw(1,4,i,j) = MIN(iLl(1),idimx) 778 ELSE 779 outLlw(1,1,i,j) = MIN(iLl(1),1) 780 outLlw(1,2,i,j) = MIN(iLl(1),1) 781 outLlw(1,3,i,j) = MAX(iLl(1)+1,idimx) 782 outLlw(1,4,i,j) = MAX(iLl(1)+1,idimx) 783 END IF 784 IF ( projlat(i,j) < inlatv(iLl(1),iLl(2)) ) THEN 785 outLlw(2,1,i,j) = MIN(iLl(1)-1,1) 786 outLlw(2,2,i,j) = MAX(iLl(1),idimy) 787 outLlw(2,3,i,j) = MAX(iLl(1),idimy) 788 outLlw(2,4,i,j) = MIN(iLl(1)-1,1) 789 ELSE 790 outLlw(2,1,i,j) = MIN(iLl(1),1) 791 outLlw(2,2,i,j) = MAX(iLl(1)+1,idimy) 792 outLlw(2,3,i,j) = MAX(iLl(1)+1,idimy) 793 outLlw(2,4,i,j) = MIN(iLl(1),1) 794 END IF 795 ! Computing distances 796 DO iv=1,4 797 ix = INT(outLlw(1,iv,i,j)) 798 iy = INT(outLlw(2,iv,i,j)) 799 dist = SQRT( (projlon(i,j)-inlonv(ix,iy))**2 + (projlat(i,j)-inlatv(ix,iy))**2 ) 800 IF ( dist /= 0.) THEN 801 outLlw(3,iv,i,j) = 1./dist 802 ELSE 803 outLlw(3,iv,i,j) = 1. 804 END IF 805 END DO 806 ! Normalizing 807 outLlw(3,:,i,j) = outLlw(3,:,i,j)/(SUM(outLlw(3,:,i,j))) 808 END DO 809 END DO 810 CASE ('npp') 811 DO i=1, pdimx 812 DO j=1, pdimy 813 ! Find point 814 difflonlat = SQRT((projlon(i,j)-inlonv)**2. + (projlat(i,j)-inlatv)**2.) 815 mindiffLl = MINVAL(difflonlat) 816 outLlw(1:2,1,i,j) = index2DArrayR(difflonlat, idimx, idimy, mindiffLl) 817 outLlw(3,1,i,j) = 1. 818 END DO 819 END DO 820 END SELECT 821 822 END SUBROUTINE LlInterpolateProjection 823 824 SUBROUTINE var2D_IntProj(var2Din, inlonv, inlatv, projlon, projlat, intkind, varout, idimx, idimy, & 825 pdimx, pdimy) 826 ! Subroutine to interpolate a 2D variable 827 828 USE module_generic 829 830 IMPLICIT NONE 831 832 INTEGER, PARAMETER :: r_k = KIND(1.d0) 833 INTEGER, INTENT(in) :: idimx, idimy, pdimx, pdimy 834 REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in) :: projlon, projlat 835 REAL(r_k), DIMENSION(idimx,idimy), INTENT(in) :: inlonv, inlatv 836 CHARACTER(LEN=50), INTENT(in) :: intkind 837 REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in) :: var2Din 838 REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(out) :: varout 839 840 ! Local 841 INTEGER :: i, j, iv, ix, iy 842 REAL(r_k) :: w 843 REAL(r_k), DIMENSION(3,16,pdimx,pdimy) :: outLlw 844 CHARACTER(LEN=50) :: fname 845 846 !!!!!!! Variables 847 ! idimx, idimy: dimension length of the input projection 848 ! pdimx, pdimy: dimension length of the target projection 849 ! in[lon/lat]: longitudes and latitudes of the target projection 850 ! proj[lon/lat]: longitudes and latitudes of the target projection 851 ! intkind: kind of interpolation 852 ! 'npp': nearest neighbourgh 853 ! 'dis': weighted distance, grid-output for SW, NW, NE, SE 854 ! outLlw: output interpolation result 855 ! for point pi,pj; up to 16 different values of 856 ! 1st: i-index in input projection 857 ! 2nd: j-index in input projection 858 ! 3rd: weight for i-index, j-index to use for ponderation (0<1.) 859 ! var2Din: 2D variable to interpolate 860 ! varout: variable interpolated on the target projection 861 fname = 'var2D_IntProj' 862 863 CALL LlInterpolateProjection(inlonv, inlatv, projlon, projlat, intkind, outLlw, idimx, idimy, pdimx,& 864 pdimy) 865 866 SELECT CASE (intkind) 867 CASE('dis') 868 DO i=1, pdimx 869 DO j=1, pdimx 870 varout(i,j) = 0. 871 DO iv=1, 4 872 ix = INT(outLlw(1,iv,i,j)) 873 iy = INT(outLlw(2,iv,i,j)) 874 w = outLlw(3,iv,i,j) 875 varout(i,j) = varout(i,j) + w*var2Din(ix,iy) 876 END DO 877 END DO 878 END DO 879 CASE('npp') 880 DO i=1, pdimx 881 DO j=1, pdimx 882 ix = INT(outLlw(1,1,i,j)) 883 iy = INT(outLlw(2,1,i,j)) 884 varout(i,j) = var2Din(ix,iy) 885 END DO 886 END DO 887 END SELECT 888 889 END SUBROUTINE var2D_IntProj 890 891 SUBROUTINE var3D_IntProj(var3Din, inlonv, inlatv, projlon, projlat, intkind, varout, idimx, idimy, & 892 pdimx, pdimy, d3) 893 ! Subroutine to interpolate a 3D variable 894 895 USE module_generic 896 897 IMPLICIT NONE 898 899 INTEGER, PARAMETER :: r_k = KIND(1.d0) 900 INTEGER, INTENT(in) :: idimx, idimy, pdimx, pdimy, d3 901 REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in) :: projlon, projlat 902 REAL(r_k), DIMENSION(idimx,idimy), INTENT(in) :: inlonv, inlatv 903 CHARACTER(LEN=50), INTENT(in) :: intkind 904 REAL(r_k), DIMENSION(pdimx,pdimy,d3), INTENT(in) :: var3Din 905 REAL(r_k), DIMENSION(pdimx,pdimy,d3), INTENT(out) :: varout 906 907 ! Local 908 INTEGER :: i, j, k, iv, ix, iy 909 REAL(r_k) :: w 910 REAL(r_k), DIMENSION(3,16,pdimx,pdimy) :: outLlw 911 CHARACTER(LEN=50) :: fname 912 913 !!!!!!! Variables 914 ! idimx, idimy: dimension length of the input projection 915 ! pdimx, pdimy: dimension length of the target projection 916 ! in[lon/lat]: longitudes and latitudes of the target projection 917 ! proj[lon/lat]: longitudes and latitudes of the target projection 918 ! intkind: kind of interpolation 919 ! 'npp': nearest neighbourgh 920 ! 'dis': weighted distance, grid-output for SW, NW, NE, SE 921 ! outLlw: output interpolation result 922 ! for point pi,pj; up to 16 different values of 923 ! 1st: i-index in input projection 924 ! 2nd: j-index in input projection 925 ! 3rd: weight for i-index, j-index to use for ponderation (0<1.) 926 ! var3Din: 3D variable to interpolate 927 ! varout: variable interpolated on the target projection 928 fname = 'var3D_IntProj' 929 930 CALL LlInterpolateProjection(inlonv, inlatv, projlon, projlat, intkind, outLlw, idimx, idimy, pdimx,& 931 pdimy) 932 933 SELECT CASE (intkind) 934 CASE('dis') 935 DO i=1, pdimx 936 DO j=1, pdimx 937 DO k=1, d3 938 varout(i,j,k) = 0. 939 DO iv=1, 4 940 ix = INT(outLlw(1,iv,i,j)) 941 iy = INT(outLlw(2,iv,i,j)) 942 w = outLlw(3,iv,i,j) 943 varout(i,j,k) = varout(i,j,k) + w*var3Din(ix,iy,k) 944 END DO 945 END DO 946 END DO 947 END DO 948 CASE('npp') 949 DO i=1, pdimx 950 DO j=1, pdimx 951 DO k=1, d3 952 ix = INT(outLlw(1,1,i,j)) 953 iy = INT(outLlw(2,1,i,j)) 954 varout(i,j,k) = var3Din(ix,iy,k) 955 END DO 956 END DO 957 END DO 958 END SELECT 959 960 END SUBROUTINE var3D_IntProj 961 962 SUBROUTINE var4D_IntProj(var4Din, inlonv, inlatv, projlon, projlat, intkind, varout, idimx, idimy, & 963 pdimx, pdimy, d3, d4) 964 ! Subroutine to interpolate a 4D variable 965 966 USE module_generic 967 968 IMPLICIT NONE 969 970 INTEGER, PARAMETER :: r_k = KIND(1.d0) 971 INTEGER, INTENT(in) :: idimx, idimy, pdimx, pdimy, d3, d4 972 REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in) :: projlon, projlat 973 REAL(r_k), DIMENSION(idimx,idimy), INTENT(in) :: inlonv, inlatv 974 CHARACTER(LEN=50), INTENT(in) :: intkind 975 REAL(r_k), DIMENSION(pdimx,pdimy,d3,d4), INTENT(in) :: var4Din 976 REAL(r_k), DIMENSION(pdimx,pdimy,d3,d4), INTENT(out) :: varout 977 978 ! Local 979 INTEGER :: i, j, k, l, iv, ix, iy 980 REAL(r_k) :: w 981 REAL(r_k), DIMENSION(3,16,pdimx,pdimy) :: outLlw 982 CHARACTER(LEN=50) :: fname 983 984 !!!!!!! Variables 985 ! idimx, idimy: dimension length of the input projection 986 ! pdimx, pdimy: dimension length of the target projection 987 ! in[lon/lat]: longitudes and latitudes of the target projection 988 ! proj[lon/lat]: longitudes and latitudes of the target projection 989 ! intkind: kind of interpolation 990 ! 'npp': nearest neighbourgh 991 ! 'dis': weighted distance, grid-output for SW, NW, NE, SE 992 ! outLlw: output interpolation result 993 ! for point pi,pj; up to 16 different values of 994 ! 1st: i-index in input projection 995 ! 2nd: j-index in input projection 996 ! 3rd: weight for i-index, j-index to use for ponderation (0<1.) 997 ! var4Din: 4D variable to interpolate 998 ! varout: variable interpolated on the target projection 999 fname = 'var4D_IntProj' 1000 1001 CALL LlInterpolateProjection(inlonv, inlatv, projlon, projlat, intkind, outLlw, idimx, idimy, pdimx,& 1002 pdimy) 1003 1004 SELECT CASE (intkind) 1005 CASE('dis') 1006 DO i=1, pdimx 1007 DO j=1, pdimx 1008 DO k=1, d3 1009 DO l=1, d4 1010 varout(i,j,k,l) = 0. 1011 DO iv=1, 4 1012 ix = INT(outLlw(1,iv,i,j)) 1013 iy = INT(outLlw(2,iv,i,j)) 1014 w = outLlw(3,iv,i,j) 1015 varout(i,j,k,l) = varout(i,j,k,l) + w*var4Din(ix,iy,k,l) 1016 END DO 1017 END DO 1018 END DO 1019 END DO 1020 END DO 1021 CASE('npp') 1022 DO i=1, pdimx 1023 DO j=1, pdimx 1024 ix = INT(outLlw(1,1,i,j)) 1025 iy = INT(outLlw(2,1,i,j)) 1026 DO k=1, d3 1027 DO l=1, d4 1028 varout(i,j,k,l) = var4Din(ix,iy,k,l) 1029 END DO 1030 END DO 1031 END DO 1032 END DO 1033 END SELECT 1034 1035 END SUBROUTINE var4D_IntProj 1036 1037 1038 SUBROUTINE var5D_IntProj(var5Din, inlonv, inlatv, projlon, projlat, intkind, varout, idimx, idimy, & 1039 pdimx, pdimy, d3, d4, d5) 1040 ! Subroutine to interpolate a 5D variable 1041 1042 USE module_generic 1043 1044 IMPLICIT NONE 1045 1046 INTEGER, PARAMETER :: r_k = KIND(1.d0) 1047 INTEGER, INTENT(in) :: idimx, idimy, pdimx, pdimy, d3, d4, d5 1048 REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in) :: projlon, projlat 1049 REAL(r_k), DIMENSION(idimx,idimy), INTENT(in) :: inlonv, inlatv 1050 CHARACTER(LEN=50), INTENT(in) :: intkind 1051 REAL(r_k), DIMENSION(pdimx,pdimy,d3,d4,d5), INTENT(in) :: var5Din 1052 REAL(r_k), DIMENSION(pdimx,pdimy,d3,d4,d5), INTENT(out):: varout 1053 1054 ! Local 1055 INTEGER :: i, j, k, l, m, iv, ix, iy 1056 REAL(r_k) :: w 1057 REAL(r_k), DIMENSION(3,16,pdimx,pdimy) :: outLlw 1058 CHARACTER(LEN=50) :: fname 1059 1060 !!!!!!! Variables 1061 ! idimx, idimy: dimension length of the input projection 1062 ! pdimx, pdimy: dimension length of the target projection 1063 ! in[lon/lat]: longitudes and latitudes of the target projection 1064 ! proj[lon/lat]: longitudes and latitudes of the target projection 1065 ! intkind: kind of interpolation 1066 ! 'npp': nearest neighbourgh 1067 ! 'dis': weighted distance, grid-output for SW, NW, NE, SE 1068 ! outLlw: output interpolation result 1069 ! for point pi,pj; up to 16 different values of 1070 ! 1st: i-index in input projection 1071 ! 2nd: j-index in input projection 1072 ! 3rd: weight for i-index, j-index to use for ponderation (0<1.) 1073 ! var5Din: 5D variable to interpolate 1074 ! varout: variable interpolated on the target projection 1075 fname = 'var5D_IntProj' 1076 1077 CALL LlInterpolateProjection(inlonv, inlatv, projlon, projlat, intkind, outLlw, idimx, idimy, pdimx,& 1078 pdimy) 1079 1080 SELECT CASE (intkind) 1081 CASE('dis') 1082 DO i=1, pdimx 1083 DO j=1, pdimx 1084 DO k=1, d3 1085 DO l=1, d4 1086 DO m=1, d5 1087 varout(i,j,k,l,m) = 0. 1088 DO iv=1, 4 1089 ix = INT(outLlw(1,iv,i,j)) 1090 iy = INT(outLlw(2,iv,i,j)) 1091 w = outLlw(3,iv,i,j) 1092 varout(i,j,k,l,m) = varout(i,j,k,l,m) + w*var5Din(ix,iy,k,l,m) 1093 END DO 1094 END DO 1095 END DO 1096 END DO 1097 END DO 1098 END DO 1099 CASE('npp') 1100 DO i=1, pdimx 1101 DO j=1, pdimx 1102 ix = INT(outLlw(1,1,i,j)) 1103 iy = INT(outLlw(2,1,i,j)) 1104 DO k=1, d3 1105 DO l=1, d4 1106 DO m=1, d5 1107 varout(i,j,k,l,m) = var5Din(ix,iy,k,l,m) 1108 END DO 1109 END DO 1110 END DO 1111 END DO 1112 END DO 1113 END SELECT 1114 1115 END SUBROUTINE var5D_IntProj 1116 716 1117 SUBROUTINE Interpolate(projlon, projlat, lonvs, latvs, mindiff, inpt, diffs, ilonlat, dimx, dimy, & 717 1118 Ninpts) -
trunk/tools/nc_var.py
r1092 r1173 44 44 'get_namelist_vars', 'get_Variables', \ 45 45 'getvalues_lonlat', 'grattr', \ 46 'grmattr', 'idims', 'igattrs', 'increaseDimvar', 'isgattrs', 'isvattrs', 'ivars',\47 'i vattrs', 'LMDZ_toCF', 'maskvar', 'model_characteristics',\46 'grmattr', 'idims', 'igattrs', 'increaseDimvar', 'isgattrs', \ 47 'isvattrs', 'ivars', 'ivattrs', 'LMDZ_toCF', 'maskvar', 'model_characteristics', \ 48 48 'ncreplace', 'ncstepdiff', 'netcdf_concatenation', 'netcdf_fold_concatenation', \ 49 'netcdf_fold_concatenation_HMT', ' Partialmap_Entiremap',\49 'netcdf_fold_concatenation_HMT', 'reproject', 'Partialmap_Entiremap', \ 50 50 'Partialmap_EntiremapFor', 'Partialmap_EntiremapForExact', \ 51 51 'pinterp', 'remapnn', \ … … 247 247 elif oper == 'Partialmap_EntiremapForExact': 248 248 ncvar.Partialmap_EntiremapForExact(opts.values, opts.ncfile, opts.varname) 249 elif oper == 'reproject': 250 ncvar.reproject(opts.values, opts.ncfile, opts.varname) 249 251 elif oper == 'seasmean': 250 252 ncvar.seasmean(timename, opts.ncfile, opts.varname) -
trunk/tools/nc_var_tools.py
r1170 r1173 90 90 # remapnn: Function to remap to the nearest neightbor a variable using projection from another file 91 91 # remapnn_old: Function to remap to the nearest neightbor a variable using projection from another fil 92 # reproject: Function to reproject values to another one 92 93 # seasmean: Function to compute the seasonal mean of a variable 93 94 # sellonlatbox: Function to select a lotlan box from a data-set … … 15782 15783 #Partialmap_Entiremap('longitude,latitude,std,5000.,lonlat', '/home/lluis/etudes/DYNAMICO/ORCHIDEE/interpolation/data/carteveg5km.nc', 'vegetation_map') 15783 15784 15785 def reproject(values, filen, variable): 15786 """ Function to reproject values to another one 15787 reproject(values, filen, variable) 15788 values=[inlonn],[inlatn],[projfile],[projlonn],[projlatn],[kind],[dimvars] 15789 inlonn= name of the longitude values in the file to interpolate 15790 inlatn= name of the latitude values in the file to interpolate 15791 projlonn= name of the longitude values in the file with the target projection 15792 projlatn= name of the latitude values in the file with the target projection 15793 kind= kind of interpolation 15794 'npp': nearest point 15795 'dis': distance weighted (4 closest) 15796 dimvars= ':', separated list of couples of dimension, and variable with its values [dimn]@[variabledim] 15797 filen= name of file to interpolate 15798 variable = ',' list of variable to interpolate ('all' for all variables) 15799 """ 15800 import module_ForInt as fin 15801 fname = 'reproject' 15802 15803 ofile = fname + '.nc' 15804 15805 if values == 'h': 15806 print fname + '_____________________________________________________________' 15807 print reproject.__doc__ 15808 quit() 15809 15810 arguments = '[inlonn],[inlatn],[projfile],[projlonn],[projlatn],[kind],[dimvars]' 15811 gen.check_arguments(fname, values, arguments, ',') 15812 15813 inlonn = values.split(',')[0] 15814 inlatn = values.split(',')[1] 15815 projfile = values.split(',')[2] 15816 projlonn = values.split(',')[3] 15817 projlatn = values.split(',')[4] 15818 kind = values.split(',')[5] 15819 dimvars = gen.Str_list(values.split(',')[6],':') 15820 15821 varns = gen.Str_list(variable,',') 15822 15823 inf = NetCDFFile(filen, 'r') 15824 if not os.path.isfile(projfile): 15825 print errormsg 15826 print ' ' + fname + ": file with target projection '" + projfile + \ 15827 "' does not exist !!" 15828 quit(-1) 15829 projf = NetCDFFile(projfile, 'r') 15830 15831 if not inf.variables.has_key(inlonn): 15832 print errormsg 15833 print ' ' +fname+ ": input file '" + filen + "' does not have longitude '" +\ 15834 inlonn + "' !!" 15835 quit(-1) 15836 if not inf.variables.has_key(inlatn): 15837 print errormsg 15838 print ' ' +fname+ ": input file '" + filen + "' does not have latitude '" + \ 15839 inlatn + "' !!" 15840 quit(-1) 15841 if not projf.variables.has_key(projlonn): 15842 print errormsg 15843 print ' ' + fname + ": file with target projection '" + projfile + \ 15844 "' does not have longitude '" + projlonn + "' !!" 15845 quit(-1) 15846 if not projf.variables.has_key(projlatn): 15847 print errormsg 15848 print ' ' + fname + ": file with target projection '" + projfile + \ 15849 "' does not have latitude '" + projlatn + "' !!" 15850 quit(-1) 15851 15852 oinlon = inf.variables[inlonn] 15853 oinlat = inf.variables[inlatn] 15854 oprojlon = projf.variables[projlonn] 15855 oprojlat = projf.variables[projlatn] 15856 15857 inlonvals0 = oinlon[:] 15858 inlatvals0 = oinlat[:] 15859 olonvals0 = oprojlon[:] 15860 olatvals0 = oprojlat[:] 15861 15862 # Getting 2D longitudes and latitudes 15863 ilonvals, ilatvals = gen.lonlat2D(ilonvals0, ilatvals0) 15864 olonvals, olatvals = gen.lonlat2D(olonvals0, olatvals0) 15865 15866 if len(ilonvals.shape) != 2: 15867 print errormsg 15868 print ' ' + fname + ': wrong input longitude rank:', ilonvals.shape, \ 15869 'it must be 2!!' 15870 quit(-1) 15871 15872 if len(ilatvals.shape) != 2: 15873 print errormsg 15874 print ' ' + fname + ': wrong input latitude rank:', ilatvals.shape, \ 15875 'it must be 2!!' 15876 quit(-1) 15877 15878 if len(olonvals.shape) != 2: 15879 print errormsg 15880 print ' ' + fname + ': wrong output longitude rank:', olonvals.shape, \ 15881 'it must be 2!!' 15882 quit(-1) 15883 15884 if len(olatvals.shape) != 2: 15885 print errormsg 15886 print ' ' + fname + ': wrong output latitude rank:', olatvals.shape, \ 15887 'it must be 2!!' 15888 quit(-1) 15889 15890 for i in range(2): 15891 if ilonvals[i] != ilatvals[i]: 15892 print errormsg 15893 print ' ' + fname + ': dimension',i,'of input longitude:', ilonvals[i], \ 15894 'differs from latitude:', ilatvals[i], 'they must be the same!!' 15895 quit(-1) 15896 if rlonvals[i] != rlatvals[i]: 15897 print errormsg 15898 print ' ' + fname + ': dimension',i,'of output longitude:', ilonvals[i],\ 15899 'differs from latitude:', ilatvals[i], 'they must be the same!!' 15900 quit(-1) 15901 15902 idx = ilonvals.shape[1] 15903 idy = ilonvals.shape[0] 15904 rdx = rlonvals.shape[1] 15905 rdy = rlonvals.shape[0] 15906 15907 ilon = ilonvals.astype('float64') 15908 ilat = ilatvals.astype('float64') 15909 rlon = rlonvals.astype('float64') 15910 rlat = rlatvals.astype('float64') 15911 15912 ilont = ilon.transpose() 15913 ilatt = ilat.transpose() 15914 rlont = rlon.transpose() 15915 rlatt = rlat.transpose() 15916 15917 # dimension-vardims 15918 vardimns = {} 15919 for vardn in dimvars: 15920 dimn = vardn.split('@')[0] 15921 dimvn = vardn.split('@')[1] 15922 vardimns[dimn] = dimvn 15923 15924 # File creation 15925 onewnc = NetCDFFile(ofile, 'w') 15926 15927 # Dimensions 15928 newdim = onewnc.createDimension('lon',rdx) 15929 newdim = onewnc.createDimension('lat',rdy) 15930 15931 # Which variables to re-project 15932 if varns[0] == 'all': 15933 varns = inf.variables.keys() 15934 # No interpolation of lon, lat 15935 varns.remove(inlonn) 15936 varns.remove(inlatn) 15937 15938 for varn in varns: 15939 if not inf.variables.has_key(inlatn): 15940 print errormsg 15941 print ' ' + fname + ": input file '" + filen + "' does not have '" + \ 15942 varn + "' !!" 15943 quit(-1) 15944 print ' ' + fname + ": re-projecting '" + varn + "' ..." 15945 ovarin = inf.variables[varn] 15946 varin = ovarin[:] 15947 15948 Nvarind = len(varin.shape) 15949 if Nvarind < 2: 15950 print errormsg 15951 print ' ' + fname + ': wrong input data rank:', varin.shape, \ 15952 'it must be at least 2!!' 15953 quit(-1) 15954 15955 varin = varin.astype('float64') 15956 varint = varin.transpose() 15957 15958 if Nvarind == 2: 15959 varout = fin.var2D_IntProj(var2Din=varint, inlonv=inlont, inlatv=inlatt, \ 15960 projlon=rlont, projlat=rlatt, intkind=kind, idimx=idx, idimy=idy, \ 15961 pdimx=rdx, pdimy=rdy) 15962 elif Nvarind == 3: 15963 varout = fin.var3D_IntProj(var3Din=varint, inlonv=inlont, inlatv=inlatt, \ 15964 projlon=rlont, projlat=rlatt, intkind=kind, idimx=idx, idimy=idy, \ 15965 pdimx=rdx, pdimy=rdy, d3=varin.shape[0]) 15966 elif Nvarind == 4: 15967 varout = fin.var4D_IntProj(var4Din=varint, inlonv=inlont, inlatv=inlatt, \ 15968 projlon=rlont, projlat=rlatt, intkind=kind, idimx=idx, idimy=idy, \ 15969 pdimx=rdx, pdimy=rdy, d3=varin.shape[0], d4=varin.shape[1]) 15970 elif Nvarind == 5: 15971 varout = fin.var4D_IntProj(var4Din=varint, inlonv=inlont, inlatv=inlatt, \ 15972 projlon=rlont, projlat=rlatt, intkind=kind, idimx=idx, idimy=idy, \ 15973 pdimx=rdx, pdimy=rdy, d3=varin.shape[0], d4=varin.shape[1], \ 15974 d5=varin.shape[2]) 15975 else: 15976 print errormsg 15977 print ' ' + fname + ': rank:', Nvarind, 'for input data not ready !!' 15978 quit(-1) 15979 15980 # Variable's dimensions and their related variables 15981 dimsfile = onewnc.dimensions() 15982 for dimn in ovarin.dimensions(): 15983 if not gen.searchInlist(dimsfile, dimn): 15984 odimn = inf.dimensions[dimn] 15985 if odimn.isunlimited(): 15986 newdim = onewnc.createDimension(dimn, None) 15987 else: 15988 newdim = onewnc.createDimension(dimn, len(odimn)) 15989 if not vardimsn.has_key(dimn): 15990 print errormsg 15991 print ' ' + fname + ": dimension '" + dimn + "' does not have "+\ 15992 'assigned any variable !!' 15993 print ' assigned pairs (via [dimvars]) _______' 15994 gen.print_dictionary(vardimsn) 15995 quit(-1) 15996 else: 15997 add_vars(inf,onewnc,vardimsn[dimn]) 15998 15999 # Interpolated variable 16000 vardims = ovarin.dimensions() 16001 newvardims = gen.replace_list(vardims,inlonn,'lon') 16002 newvardims = gen.replace_list(newvardims,inlatn,'lat') 16003 varattrs = ovarin.ncattrs() 16004 if gen.searchInlist(varattrs,'_FillValue'): 16005 fillv = ovarin.getncattr('_FillValue') 16006 newvar = onewnc.createVariable(varn, 'f4', tuple(newvardims), \ 16007 fill_value=fillv) 16008 varattrs.remove('_FillValue') 16009 else: 16010 newvar = onewnc.createVariable(varn, 'f4', tuple(newvardims)) 16011 16012 newvar[:] = varout[:] 16013 for attrn in varattrs: 16014 attrv = ovarin.getncattr(attrn) 16015 set_attribute(newvar, attrn, attrv) 16016 16017 onewnc.sync() 16018 16019 # Global attributes 16020 add_globalattrs(inf,onewnc,'all') 16021 onewnc.sync() 16022 16023 inf.close() 16024 projf.close() 16025 onewnc.close() 16026 16027 print fname + ": successful re-projection saved in '" + ofile + "' !!" 16028 16029 return 16030 15784 16031 def lonlatvarsFile(lonvals,latvals,dnx,dny,oLlnc): 15785 16032 """ Function to CF-define the longitudes and latitudes variables within a file
Note: See TracChangeset
for help on using the changeset viewer.