matplotlib - Python Basemap :Cant get map of specific coordinates -


i need show aod range in map of kathmandu valley . downloaded modis hdf file , tried show in map.but,i able plot within nepal , when try plot within kathmandu valley (coordinates latitude 27 28 , longitude 85 86 ) map of other regions not of valley.the plot shown hereplot of aod how can show within whole kathmandu valley? code follows

from matplotlib.path import path import matplotlib.patches patches file_name='mod04_l2.a2016092.0455.006.2016102170619.hdf' file = sd(file_name, sdc.read) # lat , lon info lat = file.select('latitude') latitude = lat[:] min_lat=latitude.min() max_lat=latitude.max() lon = file.select('longitude') longitude = lon[:] min_lon=longitude.min() max_lon=longitude.max()  sds=file.select(59) attributes=sds.attributes() scale_factor=attributes['scale_factor']         #get valid range aod sds range=sds.getrange() min_range=min(range) max_range=max(range)          #get sds data data=sds.get()         #get data within valid range valid_data=data.ravel() valid_data=[x x in valid_data if x>=min_range] valid_data=[x x in valid_data if x<=max_range] valid_data=np.asarray(valid_data)         #scale valid data valid_data=valid_data*scale_factor         #find average average=sum(valid_data)/len(valid_data)         #find standard deviation stdev=np.std(valid_data)         #print information print('\nthe valid range of values is: ',round(min_range*scale_factor,3), ' ',round(max_range*scale_factor,3),'\nthe average is: ',round(average,3),'\nthe standard deviation is: ',round(stdev,3)) print('the range of latitude in file is: ',min_lat,' ',max_lat, 'degrees \nthe range of longitude in file is: ',min_lon, ' ',max_lon,' degrees')   attrs = sds.attributes(full=1) fillvalue=attrs['_fillvalue']             # fillvalue[0] attribute value (-9999) fv = fillvalue[0]             #turn fillvalues nan data=data.astype(float) data[data == fv] = np.nan             #create map data = np.ma.masked_array(data, np.isnan(data)) fig = plt.figure(figsize=(20,20))  ax1 = plt.subplot2grid((2,2), (0,0)) ax2 = plt.subplot2grid((2,2), (1,0)) ax3 = plt.subplot2grid((2,2), (0,1), rowspan=2)  map1 = basemap(projection='merc', resolution='h',lat_0=27.7, lon_0=85.3, llcrnrlat=min_lat, urcrnrlat=max_lat, llcrnrlon=min_lon, urcrnrlon = max_lon, ax=ax1)#whole region map1.drawcoastlines(linewidth=1) map1.drawcountries(linewidth=1) map1.drawparallels(np.arange(-90., 120., 5), labels=[1, 0, 0, 0]) map1.drawmeridians(np.arange(-180., 181., 5), labels=[0, 0, 0, 1]) x, y = map1(longitude, latitude) p=map1.pcolormesh(x, y, data*scale_factor) plt.autoscale() fig.colorbar(p, ax=ax1)   map2 = basemap(projection='cyl',resolution='h',llcrnrlon=79,llcrnrlat=26,urcrnrlon=88.3,urcrnrlat=31, ax=ax2) #nepal map2.drawcoastlines(linewidth=1) map2.drawcountries(linewidth=1) map2.drawparallels(np.arange(-90., 120., 1), labels=[1, 0, 0, 0]) map2.drawmeridians(np.arange(-180., 181., 1), labels=[0, 0, 0, 1]) x, y = map2(longitude, latitude)  p=map2.pcolormesh(x, y, data*scale_factor) plt.autoscale() fig.colorbar(p, ax=ax2)  map3 = basemap(projection='cyl',resolution='h',llcrnrlon=85,llcrnrlat=27,urcrnrlon=86,urcrnrlat=28, ax=ax3)   #kathmandu domain map3.drawcoastlines(linewidth=1) map3.drawcountries(linewidth=1) map3.drawparallels(np.arange(-90., 120., 1), labels=[1, 0, 0, 0]) map3.drawmeridians(np.arange(-180., 181., 1), labels=[0, 0, 0, 1]) x, y = map3(longitude, latitude) map3.pcolormesh(x, y, data*scale_factor) plt.autoscale()  lbx1, lby1 = map1(*map2(map2.xmin, map2.ymin, inverse= true)) ltx1, lty1 = map1(*map2(map2.xmin, map2.ymax, inverse= true)) rtx1, rty1 = map1(*map2(map2.xmax, map2.ymax, inverse= true)) rbx1, rby1 = map1(*map2(map2.xmax, map2.ymin, inverse= true))  verts1 = [     (lbx1, lby1), # left, bottom     (ltx1, lty1), # left, top     (rtx1, rty1), # right, top     (rbx1, rby1), # right, bottom     (lbx1, lby1), # ignored     ]  codes2 = [path.moveto,          path.lineto,          path.lineto,          path.lineto,          path.closepoly,          ]  path = path(verts1, codes2) patch = patches.pathpatch(path, lw=2) ax1.add_patch(patch)    fig = plt.gcf()             # show plot window. plt.show() plt.savefig('plot.png')  print('\nall valid files have been processed.') 


Comments

Popular posts from this blog

php - Vagrant up error - Uncaught Reflection Exception: Class DOMDocument does not exist -

vue.js - Create hooks for automated testing -

Add new key value to json node in java -