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
Post a Comment