program rd_wr_binary implicit none integer :: istatus character (len=256) :: name_land,name_ufrac real, allocatable, dimension(:,:) :: xlat_o,xlong_o,ufrac_o real, allocatable, dimension(:) :: xlat_i,xlong_i,ufrac_i real, allocatable, dimension(:) :: landuse_i real, allocatable, dimension(:,:) :: landuse_o,rarray integer :: nx ! x-dimension of the array integer :: ny ! y-dimension of the array integer :: nz ! z-dimension of the array integer :: isigned ! 0=unsigned data, 1=signed data integer :: endian ! 0=big endian, 1=little endian real :: scalefactor ! value to divide array elements by before truncation to integers integer :: wordsize ! number of bytes to use for each array element real :: latmin,longmin,latmax,longmax,dx,dy,latk,lonk,laty,lonx,pi,umax,a1,a2,a3 real*8 xx,yy,xlat,xlong real :: dlat, dlon real :: conv double precision a,b,c,d,rt,dist,dmin real xlongmax,utot_i,utot_o integer :: np,ip,ipmin,ix,iy,i,j,iip,iizone,ispher,iutm integer nurbm,icmax,imax integer igsize ! The following must be set before compiling ! Change nx, ny (same (-1) as in namelist.wps) and np (number of values extracted) character*70 aa ! parameter (nx=240,ny=270,nz=1,np=500544,nurbm=10) parameter (nx=696,ny=674,nz=1,np=3342500,nurbm=10) integer icount(nurbm),itot,itot_urb name_land='landuse_urban' name_ufrac='ufrac' ! Implement coordinates of the domain regarding the Region Of Interest latmin=21.8524 longmin=112.82 latmax=23.538 longmax=114.56 ! Did you extract the data in UTM or in WGS84? ! flag input data. If iutm=1 input data in UTM after extraction, if iutm=0 input data in WGS84 (lat long) iutm=1 ! Chose your interpolation : "Looking for the closest" is great if you have a domain grid size between ! 100-300m and if the domain surface is below 10 000 km2. "Most frequent value" is great with grid sizes ! greater than 250m, there is no limit for the domain surface ! igzise flag interpolation - if igsize.eq.0 look for the closest, if igsize.eq.1 most frequent value igsize=0 !! isigned = 0 !! endian = 0 !! wordsize = 1 !! scalefactor = 1.0 ! isigned = 1 ! endian = 1 ! wordsize = 4 ! scalefactor = 0.01 ! Do not change anything here isigned = 1 endian = 1 wordsize = 1 scalefactor = 1 allocate(ufrac_o(nx,ny)) allocate(xlat_o(nx,ny)) allocate(xlong_o(nx,ny)) allocate(landuse_o(nx,ny)) allocate(rarray(nx,ny)) ufrac_o=-9999 landuse_o=0 allocate(xlat_i(np)) allocate(xlong_i(np)) allocate(ufrac_i(np)) allocate(landuse_i(np)) dx=(longmax-longmin)/nx dy=(latmax-latmin)/ny write(*,*)'dx=',dx,'dy=',dy rt=6371000. pi=3.1415926535 xlongmax=-100000000000000000. ! If iutm=1 change the iizone regarding the UTM zone of the domain iizone=50 ispher=4 ! Change the file=[name of the city].txt open(unit=20,file='shenzhen_clip4.txt',status='old') read(20,*)aa ! write(*,*)aa do ip=1,np read(20,*) iip,xx,yy,landuse_i(ip) if(iutm.eq.1)then call utm2ll(xx,yy,xlat,xlong,conv,iizone,ispher) xlong_i(ip)=-xlong/3600. xlat_i(ip)=xlat/3600. elseif(iutm.eq.0)then xlong_i(ip)=xx xlat_i(ip)=yy endif ! write(*,*)xx,yy,-xlong_i(ip)/3600.,xlat_i(ip)/3600. enddo close(20) do ix=1,nx do iy=1,ny xlat_o(ix,iy)=latmin+dy*(iy-1) xlong_o(ix,iy)=longmin+dx*(ix-1) enddo enddo landuse_o(:,:)=0 if(igsize.eq.0)then write(*,*)'start looking for the closest' ipmin=-1 do ix=1+5,nx-5 do iy=1+5,ny-5 dmin=1000000. do ip=1,np latk=xlat_i(ip) lonk=xlong_i(ip) laty=xlat_o(ix,iy) lonx=xlong_o(ix,iy) a=cos((latk+laty)*pi/360.) dlat=rt*(latk-laty)*pi/180. dlon=rt*a*(lonk-lonx)*pi/180. dist=(dlat**2.+dlon**2.)**.5 if(dist.lt.dmin)then dmin=dist ipmin=ip endif enddo if(landuse_i(ipmin).lt.11.and.landuse_i(ipmin).ge.1)then landuse_o(ix,iy)=landuse_i(ipmin)+30 else landuse_o(ix,iy)=0 endif if(landuse_o(ix,iy).eq.30)write(*,*)'Not Working at all',ix,iy,landuse_i(ipmin) ! endif enddo enddo ! elseif(igsize.eq.1)then write(*,*)'start calculating most frequent value' do ix=1+5,nx-5 do iy=1+5,ny-5 itot=0 itot_urb=0 icount(:)=0 do ip=1,np latk=xlat_i(ip) lonk=xlong_i(ip) laty=xlat_o(ix,iy) lonx=xlong_o(ix,iy) if((xlat_i(ip).gt.(xlat_o(ix,iy)-dx/2.)).and.(xlat_i(ip).le.(xlat_o(ix,iy)+dx/2.)))then if((xlong_i(ip).gt.(xlong_o(ix,iy)-dy/2.)).and.(xlong_i(ip).le.(xlong_o(ix,iy)+dy/2.)))then itot=itot+1 if((landuse_i(ip).ge.1).and.(landuse_i(ip).lt.11))then icount(landuse_i(ip))=icount(landuse_i(ip))+1 itot_urb=itot_urb+1 endif endif endif enddo icmax=0 do i=1,nurbm if(icount(i).gt.icmax)then icmax=icount(i) imax=i endif enddo write(66,*)'icount=',icount,imax ! If you change the urban fraction condition the urban surface cover will change. ! Default is set to "greater than 0.5" if(icmax.gt.0.and.(float(itot_urb)/float(itot).gt.0.5))then ! if(icmax.gt.0)then landuse_o(ix,iy)=imax+30 ! ufrac_o(ix,iy)=float(itot_urb)/float(itot) else landuse_o(ix,iy)=0 endif write(66,*)'landuse_o(ix,iy)=',landuse_o(ix,iy) enddo enddo endif ! CAREFUL : Don't change the routine and subroutine down below !!!!!!!!!!! utot_o=0. utot_i=0. do i=1,nx do j=1,ny utot_o=utot_o+ufrac_o(i,j)*1000.*1000./3./3. enddo enddo write(*,*)utot_o,5.e06/utot_o*1.e4 do i=1,np utot_i=utot_i+ufrac_i(i)*1000.*1000./3./3. enddo write(*,*)utot_i,5.e06/utot_i*1.e4 ! ! Modify the field as necessary ! ! ! Write data to geogrid binary format using write_geogrid() ! isigned = 1 endian = 1 wordsize = 1 scalefactor = 1 open(unit=44,file='prueba') do ix=1,nx do iy=1,ny write(44,*)ix,iy,landuse_o(ix,iy) enddo enddo call write_geogrid((name_land), len_trim(name_land), landuse_o, nx, ny, nz, isigned, endian, scalefactor, wordsize) isigned = 1 endian = 1 wordsize = 4 scalefactor = 0.01 call write_geogrid((name_ufrac), len_trim(name_ufrac), ufrac_o, nx, ny, nz, isigned, endian, scalefactor, wordsize) ! ! Read data from geogrid binary format using read_geogrid() ! ! call read_geogrid('lll', len_trim('lll'), rarray, nx, ny, nz, isigned, endian, scalefactor, wordsize, istatus) ! if (istatus /= 0) then ! write(0,*) 'Error while reading '//trim('lll')//'. Quitting.' ! end if do iy=ny,1,-5 write(44,'(49(f6.3,1x))')(landuse_o(ix,iy),ix=1,nx,5) enddo deallocate(ufrac_o) deallocate(landuse_o) ! do ix=1,nx ! do iy=1,ny ! if(rarray(ix,iy).eq.31)write(*,*)'31',ix,iy ! enddo ! enddo end program rd_wr_binary subroutine utm2ll(x,y,slat,slon,conv,iizone,ispher) ! ! universal transverse mercator conversion ! ! rewritten 6-16-83 by j.f. waananen using formulation ! by j.p. snyder in bulletin 1532, pages 68-69 ! ! convert utm.f to convert from utm to lat lon (idir > 0) ! 2/5/2001 j. klinck ! implicit none ! real*8 axis(19),bxis(19) real*8 radsec real*8 ak0,a,b,es real*8 x,y,slat,slon,conv,cm,phi,dlam,epri,en,t real*8 c,em,xx,yy,um,e1 real*8 phi1,r,d,alam ! integer iizone,ispher,izone integer iutz ! data axis/6378206.4d0,6378249.145d0,6377397.155d0, & 6378157.5d0,6378388.d0,6378135.d0,6377276.3452d0, & 6378145.d0,6378137.d0,6377563.396d0,6377304.063d0, & 6377341.89d0,6376896.0d0,6378155.0d0,6378160.d0, & 6378245.d0,6378270.d0,6378166.d0,6378150.d0/ ! data bxis/6356583.8d0,6356514.86955d0,6356078.96284d0, & 6356772.2d0,6356911.94613d0,6356750.519915d0,6356075.4133d0, & 6356759.769356d0,6356752.31414d0,6356256.91d0,6356103.039d0, & 6356036.143d0,6355834.8467d0,6356773.3205d0,6356774.719d0, & 6356863.0188d0,6356794.343479d0,6356784.283666d0, & 6356768.337303d0/ ! data ak0/0.9996d0/ ! data radsec/206264.8062470964d0/ ! a = axis(ispher) b = bxis(ispher) es= (a**2-b**2)/a**2 izone = iizone ! compute utm zone(izone) and central meridian in seconds for ! geodetic to utm conversion where zone is not input. ! if (izone.eq.0 ) then write(*,*) ' ************* error exit from utm. ' write(*,*) ' zone is not given for conversion from ' write(*,*) ' x,y to lat,lon. ' return endif ! if(iabs(izone).gt.30) then iutz = iabs(izone)-30 cm=((iutz*6.0d0)-3.0d0)*(-3600.0d0) else iutz = 30-iabs(izone) cm=((iutz*6.0d0)+3.0d0)*3600.0d0 endif ! !--------------------------------------------------------------- ! inverse computation !--------------------------------------------------------------- yy = y if (izone.lt.0) yy = yy - 1.0d7 xx = x - 5.0d5 em = yy/ak0 um = em/(a*(1.d0-(es/4.d0)-(3.d0*es*es/6.4d1)- & (5.d0*es*es*es/2.56d2))) e1 = (1.d0-dsqrt(1.d0-es))/(1.d0+dsqrt(1.d0-es)) ! phi1 = um+((3.d0*e1/2.d0)-(2.7d1*e1**3/3.2d1))*dsin(2.d0*um) + & ((2.1d1*e1*e1/1.6d1)-(5.5d1*e1**4/3.2d1))*dsin(4.d0*um) + & (1.51d2*e1**3/9.6d1)*dsin(6.d0*um) ! en = a/dsqrt(1.0d0-es*dsin(phi1)**2) t = dtan(phi1)**2 epri = es/(1.d0-es) c = epri*dcos(phi1)**2 r = (a*(1.d0-es))/((1.d0-es*dsin(phi1)**2)**1.5d0) d = xx/(en*ak0) ! phi = phi1 - (en*dtan(phi1)/r) * ((d*d/2.d0) - & (5.d0+3.d0*t+10.d0*c-4.d0*c*c-9.d0*epri)*d**4/2.4d1 & + (6.1d1+9.d1*t+2.98d2*c+4.5d1*t*t & -2.52d2*epri-3.d0*c*c)*d**6/7.2d2) ! alam = (cm/radsec)-(d-(1.d0+2.d0*t+c)*d**3/6.d0 + & (5.d0-2.d0*c+2.8d1*t -3.d0*c*c+8.d0*epri+2.4d1*t*t) & *d**5/1.2d2)/dcos(phi1) ! slat = phi*radsec slon = alam*radsec dlam = -(slon-cm)/radsec !--------------------------------------------------------------- ! ! test to see if within definition of utm projection ! ! ---- latitude within 84 degrees if (dabs(slat).gt.302400) then write(99,*) ' ************* error exit from utm. ' write(99,*) ' latitude value is poleward of 84 degrees' write(99,*) ' utm is not valid.' write(99,*) ' calculation has continued but values ' write(99,*) ' may not be valid.' write(99,*) slat/3600.d0, slon/3600.d0 endif ! ---- delta longitude within 0.16 radians of central meridian if (dabs(dlam).gt.1.6d-1) then write(99,*) ' ************* error exit from utm. ' write(99,*) ' d lon not within 0.16 radians of central merid.' write(99,*) ' calculation has continued but values ' write(99,*) ' may not be valid.' write(99,*) slat/3600.d0, slon/3600.d0 endif ! ! ! compute convergence ! conv = dlam*(dsin(phi)+1.9587d-12*(dlam**2) & *dsin(phi)*dcos(phi*phi)) conv = conv*radsec ! return end