2013-10-22 6 views
2

NetCDF 파일의 데이터를 세계지도에 플롯하려고합니다. 그러나 내가 할 때마다 문제가있는 음모가 있습니다 (스크린 샷 참조). 최소값과 최대 값이 전체 데이터 범위를 포괄하므로 완전히 채색되어야하며 위도와 경도는 필요한대로 인쇄됩니다. 그것은 ncview에서 완벽하게 잘 나타내지 만 파이썬에서는 그렇지 않습니다. 데이터를지도에 올바르게 표시하려면 어떻게해야합니까? (AKA I가 청소 만 가지고는 어디?)데이터가 단편화되지 않도록하려면 어떻게합니까?

코드는 다음과 같습니다 :

import matplotlib.pyplot as plt 
import matplotlib.cm as cm 
import matplotlib as mpl 
import numpy as np 
import scipy.io.netcdf as S 
from mpl_toolkits.basemap import * 
import pyproj as py 
#import various libraries needed 

q = input("Enter month of data wanted (1 = Jan, 12 = Dec): ") 
#make sure that input is valid (integer 1 - 12) - Warn if not and close program 
if q>12: 
    print "You have not entered a valid month!" 
elif q<1: 
    print "You have not entered a valid month!" 
elif isinstance(q, (int, long)) == False: 
    print "You're just trying to annoy me now, aren't you?" 
    #if all is good, proceed to getting data 
else: 


qq = q-1 
f = S.netcdf_file('aerocom.HadGEM3-A-GLOMAP.AEROCOM-A2-CTRL.monthly.sconcbc.2006.nc','r') 
#read in the file as a netcdf file into a filehandle with data 
bc = f.variables['sconcbc'][qq,:,:] 
lon1 = f.variables['lon'][:] 
lat = f.variables['lat'][:] 
bc, lon = addcyclic(bc, lon1) 

print lon 
for x in np.nditer(lon, op_flags=['readwrite']): 
    if x>180 : 
     x[...] = x-360 
print lon 
print lat  

# use low resolution coastlines. 
map = Basemap(projection='cyl',lat_0=0.,lon_0=0,resolution='l') 
map.fillcontinents #(color='grey',lake_color='white') 
# draw the edge of the map projection region 
map.drawmapboundary(fill_color='white') 
#draw lat/lon grid lines every 30 degrees. 
map.drawmeridians(np.arange(-180,180,30)) 
map.drawparallels(np.arange(-90,90,30)) 

#make grid on world map from the long/lat in the data (shape of bc) 
lon, lat = np.meshgrid(lon, lat) 
    #create array for colourbar and contour pattern 
    Y =[1e-15,1e-14,1e-13,5e-13,1e-12,3e-12,5e-12,7e-12,9e-12,1e-11,3e-11,5e-11,6e-11,7e-11,9e-11,1e-10,3e-10,5e-10,7e-10,9e-10,1e-9,2e-9,3e-9,5e-9,7e-9,9e-9,1e-8,2e-8,4e-8,6e-8,8e-8] 
#make those the x and y axis, but in map form 
x, y =map(lon,lat) 


plt.clf() 
my_cmap = cm.get_cmap('CMRmap') 


cs = map.contourf(x,y,bc,levels = Y,cmap=my_cmap,locator=mpl.ticker.LogLocator()) 


# set colourbar with location and size, with labels. 


plt.colorbar(cmap=my_cmap,orientation="horizontal",shrink=0.5) 

font = {'family' : 'serif', 
    'color' : 'black', 
    'weight' : 'bold', 
    'size' : 24, 
    } 
#add plot details 
plt.title('Black Carbon monthly average surface concentrations for %s/2006 in kgm-3'%(q) ,fontdict=font) 
map.drawcoastlines(linewidth=0.75) 
map.drawcountries(linewidth=0.25) 
#show plot 
plt.show() 

스크린 샷 :

Black Carbon monthly average surface concentrations for 9/2006

dimensions: 
    time = UNLIMITED ; // (12 currently) 
    lat = 145 ; 
    lon = 192 ; 
    bnds = 2 ; 
variables: 
    double time(time) ; 
      time:bounds = "time_bnds" ; 
      time:units = "days since 1850-01-01" ; 
      time:calendar = "360_day" ; 
      time:axis = "T" ; 
      time:long_name = "time" ; 
      time:standard_name = "time" ; 
    double time_bnds(time, bnds) ; 
    double lat(lat) ; 
      lat:bounds = "lat_bnds" ; 
      lat:units = "degrees_north" ; 
      lat:axis = "Y" ; 
      lat:long_name = "latitude" ; 
      lat:standard_name = "latitude" ; 
    double lat_bnds(lat, bnds) ; 
    double lon(lon) ; 
      lon:bounds = "lon_bnds" ; 
      lon:units = "degrees_east" ; 
      lon:axis = "X" ; 
      lon:long_name = "longitude" ; 
      lon:standard_name = "longitude" ; 
    double lon_bnds(lon, bnds) ; 
    float sconcbc(time, lat, lon) ; 
      sconcbc:standard_name = "mass_concentration_of_black_carbon_dry_ 
aerosol_in_air" ; 
      sconcbc:long_name = "Surface concentration BC" ; 
      sconcbc:units = "kg m-3" ; 
      sconcbc:original_name = "mo: m01s77i008" ; 
      sconcbc:cell_methods = "time: mean" ; 
      sconcbc:history = "2012-04-18T10:44:22Z altered by CMOR: replace 
d missing value flag (-1.07374e+09) with standard missing value (1e+20)." ; 
      sconcbc:missing_value = 1.e+20f ; 
      sconcbc:_FillValue = 1.e+20f ; 
      sconcbc:associated_files = "gridspecFile: gridspec_REALM_fx_HadG 
EM3-A-GLOMAP_AEROCOM-A2-CTRL_r0i0p0.nc" ; 
+0

데이터를 가지고 있지 않은 곳에 'contourf' 또는'imshow' 보간과 같은 알고리즘을 사용할 수 있습니다 ... –

+0

dateline을 가로 지르는 다각형이 360도 회전하고 그들 자신과 교차하는 것처럼 보입니다. 그러나 실제로 데이터가 어떻게 보이는지 알지 못해도 확신하기가 어렵습니다. 공개 데이터입니까? –

+0

사이드 노트 : 스크린 샷을 찍는 대신 "파일 저장"버튼을 사용할 수 있습니다. –

답변

0

솔루션을 찾을 수 :

문제는 함께 있었다 0-360 위도 설정에서 -180에서 180까지 데이터를 가져 오는 루프.이 작업을 수행함에 따라 t는 위도를 올바르게 유지하지만 단조롭게 증가하지는 않습니다 (180에서 1.25까지는 -178.75에서 0으로 건너 뜁니다)

이것은 주기적 특성과 더불어 점을 그려 플롯 한 다음 등고선을 0 (최종 위도)과 180 (첫 번째 위도) 데이터를 서로 맞 춥니 다. 이것은 물론 서로의 중간에 위치합니다.

이 문제를 해결하려면 위도 줄을 정렬하고 데이터를 다시 정렬하여 새로 정렬 된 위도와 일치 시키십시오.

lon = sorted(lon) 


bc1 = bc[:,0:97] 
bc2 = bc[:,97:] 

bc = np.hstack((bc2,bc1)) 

이렇게하면 데이터와 위도가 올바른 순서로 재생성되고 플로팅 오류가 방지됩니다.

관련 문제