Menu
▾
▴
matplotlib-checkins
|
From: <js...@us...> - 2008-08-06 22:04:10
|
Revision: 5985
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=5985&view=rev
Author: jswhit
Date: 2008-08-06 22:04:08 +0000 (Wed, 06 Aug 2008)
Log Message:
-----------
update warpimage, include plotsst
Modified Paths:
--------------
trunk/toolkits/basemap/examples/run_all.py
trunk/toolkits/basemap/examples/warpimage.py
Modified: trunk/toolkits/basemap/examples/run_all.py
===================================================================
--- trunk/toolkits/basemap/examples/run_all.py 2008-08-06 21:57:19 UTC (rev 5984)
+++ trunk/toolkits/basemap/examples/run_all.py 2008-08-06 22:04:08 UTC (rev 5985)
@@ -2,7 +2,6 @@
test_files = glob.glob('*.py')
test_files.remove('run_all.py')
test_files.remove('fcstmaps.py')
-test_files.remove('plotsst.py')
test_files.remove('testgdal.py')
test_files.remove('pnganim.py')
test_files.remove('geos_demo_2.py')
Modified: trunk/toolkits/basemap/examples/warpimage.py
===================================================================
--- trunk/toolkits/basemap/examples/warpimage.py 2008-08-06 21:57:19 UTC (rev 5984)
+++ trunk/toolkits/basemap/examples/warpimage.py 2008-08-06 22:04:08 UTC (rev 5985)
@@ -42,8 +42,8 @@
fig=plt.figure()
# define orthographic projection centered on Europe.
m = Basemap(projection='ortho',lat_0=40,lon_0=40,resolution='l')
-# plot warped rgba image.
-im = m.bluemarble()
+# plot a gray-scale image specified from a URL.
+im = m.warpimage("http://earthobservatory.nasa.gov/Newsroom/BlueMarble/Images/gebco_bathy.5400x2700.jpg")
# draw coastlines.
m.drawcoastlines(linewidth=0.5,color='0.5')
# draw lat/lon grid lines every 30 degrees.
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
![http://earthobservatory.nasa.gov/Newsroom/BlueMarble/Images/gebco_bathy.5400x2700.jpg")](http://earthobservatory.nasa.gov/Newsroom/BlueMarble/Images/gebco_bathy.5400x2700.jpg")){kind=link}
|
From: <js...@us...> - 2008-08-07 18:02:35
|
Revision: 5991
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=5991&view=rev
Author: jswhit
Date: 2008-08-07 18:02:32 +0000 (Thu, 07 Aug 2008)
Log Message:
-----------
update testgdal example (now exercises more gdal features)
Modified Paths:
--------------
trunk/toolkits/basemap/examples/README
trunk/toolkits/basemap/examples/testgdal.py
Added Paths:
-----------
trunk/toolkits/basemap/examples/us_25m.dem
Modified: trunk/toolkits/basemap/examples/README
===================================================================
--- trunk/toolkits/basemap/examples/README 2008-08-07 18:00:25 UTC (rev 5990)
+++ trunk/toolkits/basemap/examples/README 2008-08-07 18:02:32 UTC (rev 5991)
@@ -73,9 +73,8 @@
hurrtrack.py plots hurricane tracks from shapefile data
(from nationalatlas.gov).
-testgdal.py is an example that shows how to plot raster geospatial data
-using the gdal module (http://gdal.maptools.org). Requires two data
-files that must be downloaded manually, see docstrings for URLs.
+testgdal.py is an example that shows how to plot GIS data
+using the gdal module (http://gdal.maptools.org).
panelplot.py shows how to make multi-panel plots, taking special care to
make sure that map aspect ratios are preserved (so you don't get squished
Modified: trunk/toolkits/basemap/examples/testgdal.py
===================================================================
--- trunk/toolkits/basemap/examples/testgdal.py 2008-08-07 18:00:25 UTC (rev 5990)
+++ trunk/toolkits/basemap/examples/testgdal.py 2008-08-07 18:02:32 UTC (rev 5991)
@@ -1,67 +1,58 @@
"""
-example showing how to plot a USGS DEM file using
+example showing how to plot data from a DEM file and an ESRI shape file using
gdal (http://pypi.python.org/pypi/GDAL).
-
-Data files must be downloaded manually from USGS:
-http://edcftp.cr.usgs.gov/pub/data/DEM/250/D/denver-w.gz
-http://edcftp.cr.usgs.gov/pub/data/nationalatlas/countyp020.tar.gz
"""
from osgeo import gdal, ogr
-from mpl_toolkits.basemap import Basemap
+from mpl_toolkits.basemap import Basemap, cm
import numpy as np
import matplotlib.pyplot as plt
+from numpy import ma
-# download from
-# http://edcftp.cr.usgs.gov/pub/data/DEM/250/D/denver-w.gz
-gd = gdal.Open('denver-w')
+# read DEM file using gdal.
+gd = gdal.Open('us_25m.dem')
# get data from DEM file
array = gd.ReadAsArray()
# get lat/lon coordinates from DEM file.
coords = gd.GetGeoTransform()
-llcrnrlon = coords[0]
-urcrnrlon = llcrnrlon+(array.shape[1]-1)*coords[1]
-urcrnrlat = coords[3]
-llcrnrlat = urcrnrlat+(array.shape[0]-1)*coords[5]
-# create Basemap instance.
-m = Basemap(llcrnrlon=llcrnrlon,llcrnrlat=llcrnrlat,urcrnrlon=urcrnrlon,urcrnrlat=urcrnrlat,projection='cyl')
-# create a figure, add an axes
-# (leaving room for a colorbar).
-fig = plt.figure()
-ax = fig.add_axes([0.1,0.1,0.75,0.75])
-# plot image from DEM over map.
-im = m.imshow(array,origin='upper')
-# make a colorbar.
-cax = plt.axes([0.875, 0.1, 0.05, 0.75]) # setup colorbar axes.
-plt.colorbar(cax=cax) # draw colorbar
-plt.axes(ax) # make the original axes current again
+nlons = array.shape[1]
+nlats = array.shape[0]
+delon = coords[1]
+delat = coords[5]
+lons = coords[0] + delon*np.arange(nlons)
+lats = coords[3] + delat*np.arange(nlats)[::-1]
+# setup basemap instance.
+m = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,
+ projection='lcc',lat_1=33,lat_2=45,lon_0=-95)
+# transform to nx x ny regularly spaced (4 km) native projection grid
+nx = int((m.xmax-m.xmin)/4000.)+1; ny = int((m.ymax-m.ymin)/4000.)+1
+topoin = ma.masked_values(array[::-1,:],-999.)
+topodat = m.transform_scalar(topoin,lons,lats,nx,ny,masked=True)
+im = m.imshow(topodat,cmap=cm.GMT_haxby_r)
# draw meridians and parallels.
-m.drawmeridians(np.linspace(llcrnrlon+0.1,urcrnrlon-0.1,5),labels=[0,0,0,1],fmt='%4.2f')
-m.drawparallels(np.linspace(llcrnrlat+0.1,urcrnrlat-0.1,5),labels=[1,0,0,0],fmt='%4.2f')
-# plot county boundaries from
-# http://edcftp.cr.usgs.gov/pub/data/nationalatlas/countyp020.tar.gz
-g = ogr.Open ("countyp020.shp")
-L = g.GetLayer(0)
-for feat in L:
- field_count = L.GetLayerDefn().GetFieldCount()
- geo = feat.GetGeometryRef()
- if geo.GetGeometryCount()<2:
- g1 = geo.GetGeometryRef( 0 )
- x =[g1.GetX(i) for i in range(g1.GetPointCount()) ]
- y =[g1.GetY(i) for i in range(g1.GetPointCount()) ]
- m.plot(x,y,'k')
- for count in range( geo.GetGeometryCount()):
- geom = geo.GetGeometryRef ( count )
- for cnt in range( geom.GetGeometryCount()):
- g1 = geom.GetGeometryRef( cnt )
- x =[g1.GetX(i) for i in range(g1.GetPointCount()) ]
- y =[g1.GetY(i) for i in range(g1.GetPointCount()) ]
- m.plot(x,y,'k')
-# plot some cities.
-lons = [-105.22,-105.513,-105.316,-105.47]; lats = [39.76,39.801,39.633,39.41]
-names = ['Golden','Central City','Evergreen','Bailey']
-x,y = m(lons,lats)
-m.plot(x,y,'ko')
-for name,xx,yy in zip(names,x,y):
- plt.text(xx+0.01,yy+0.01,name)
-plt.title(gd.GetDescription()+' USGS DEM with county boundaries')
+m.drawparallels(np.arange(25,65,20),labels=[1,0,0,0])
+m.drawmeridians(np.arange(-120,-40,20),labels=[0,0,0,1])
+# plot state boundaries from shapefile using ogr.
+g = ogr.Open ("st99_d00.shp")
+L = g.GetLayer(0) # data is in 1st layer.
+for feat in L: # iterate over features in layer
+ geo = feat.GetGeometryRef()
+ # iterate over geometries.
+ for count in range(geo.GetGeometryCount()):
+ geom = geo.GetGeometryRef(count)
+ if not geom.GetGeometryCount(): # just on geometry.
+ # get lon,lat points
+ lons = [geom.GetX(i) for i in range(geom.GetPointCount())]
+ lats = [geom.GetY(i) for i in range(geom.GetPointCount())]
+ # convert to map projection coords.
+ x, y = m(lons,lats)
+ # plot on map.
+ m.plot(x,y,'k')
+ else: # iterate over nested geometries.
+ for cnt in range( geom.GetGeometryCount()):
+ g1 = geom.GetGeometryRef( cnt )
+ lons = [g1.GetX(i) for i in range(g1.GetPointCount())]
+ lats = [g1.GetY(i) for i in range(g1.GetPointCount())]
+ x, y = m(lons,lats)
+ m.plot(x,y,'k')
+plt.title(gd.GetDescription()+' with state boundaries from '+g.GetName())
plt.show()
Added: trunk/toolkits/basemap/examples/us_25m.dem
===================================================================
--- trunk/toolkits/basemap/examples/us_25m.dem (rev 0)
+++ trunk/toolkits/basemap/examples/us_25m.dem 2008-08-07 18:02:32 UTC (rev 5991)
@@ -0,0 +1,627 @@
+ncols 1405
+nrows 621
+xllcorner -125.020833333337
+yllcorner 24.062499999999676
+cellsize 0.04166667
+NODATA_value -9999
+0 0 0 0 0 0 0 0 3 35 153 114 174 496 309 627 251 218 453 858 1165 1025 352 108 590 35 75 251 668 919 1421 1302 1204 1553 1459 1368 1761 1709 1229 1425 425 241 960 616 433 1162 1429 1715 1624 1949 1925 1668 1656 1712 1771 1304 1289 1340 1626 1659 1424 1051 505 754 282 1085 2019 1877 1656 1792 1335 913 917 1050 903 1332 1211 1716 1673 1187 1701 1638 1667 1501 1164 278 822 1170 1618 1569 1325 1543 1255 1108 1233 993 1338 1227 898 1024 1032 1038 1054 1235 1316 1075 1110 1244 1189 1288 1389 1556 1628 1528 1567 1579 1601 1659 1620 1695 1583 1497 1460 1177 1343 1322 1357 1172 1453 1156 982 482 382 484 446 449 779 1210 1336 1452 1275 1108 1287 1192 1389 1645 1818 1805 1895 1810 1610 1210 1288 1622 1689 1707 1765 1875 2053 2043 1685 1689 1262 638 833 668 661 627 1229 1999 1441 944 1599 1615 1761 2164 2075 2309 2072 1903 1783 1296 601 829 1248 1509 1502 1803 2054 2094 1423 1263 1105 907 632 685 1486 1619 1866 1967 2078 2276 1946 1932 1715 1769 2376 2419 2073 2342 1830 2057 2232 2163 1716 2232 2026 2037 1368 1178 1121 858 796 862 1120 1213 2297 2094 1774 2070 1904 2063 2259 1763 2221 2297 1711 1881 1905 1948 2062 1910 1363 1262 1689 1685 1590 2033 1894 2100 1781 1732 1594 1640 1582 1666 1708 1532 1515 1614 1410 1356 1541 1583 1456 1460 1385 1305 1242 1117 1070 1029 1002 992 983 992 993 992 986 971 963 963 957 957 963 967 972 977 981 984 952 941 940 917 887 871 868 863 855 831 815 844 861 860 836 816 805 783 761 783 762 762 767 741 753 753 749 729 724 747 751 782 744 757 725 714 748 770 780 795 809 830 823 828 838 826 798 781 766 748 736 734 739 778 761 778 799 828 826 835 832 826 803 846 864 843 827 795 778 755 788 781 794 799 812 832 851 818 813 787 769 761 762 769 776 778 778 795 785 788 799 790 831 821 848 853 854 854 860 905 972 985 957 929 915 920 906 870 882 900 912 881 873 901 902 902 915 920 936 920 866 892 883 879 873 876 876 872 857 818 765 764 753 732 725 735 750 724 717 722 724 724 720 729 747 736 734 721 713 713 712 703 701 697 695 708 716 716 718 731 745 744 732 712 697 679 690 699 704 706 710 696 679 693 700 692 730 777 734 724 746 780 825 836 796 730 697 716 710 659 636 621 613 619 630 634 615 603 589 585 585 581 579 576 579 579 579 574 572 566 569 573 574 583 585 585 585 591 592 597 597 597 597 599 603 603 608 617 626 633 646 654 650 640 631 625 627 633 643 670 688 707 729 730 735 734 739 745 744 756 740 723 717 712 708 703 690 676 672 667 660 662 647 643 623 612 603 597 591 587 587 590 586 578 558 540 532 526 522 514 503 498 490 492 478 469 462 448 444 416 419 451 452 453 452 436 416 404 410 388 375 368 363 386 412 416 398 406 413 418 413 406 400 396 394 393 393 385 383 384 384 387 381 385 384 384 382 383 377 365 366 376 366 353 336 334 329 311 302 299 293 292 292 287 286 286 284 281 279 276 276 263 256 254 253 249 249 249 244 240 237 237 237 237 237 237 237 237 234 233 237 237 237 237 232 236 232 231 231 224 226 229 231 234 237 237 237 237 237 238 247 252 256 260 262 262 262 260 264 274 280 285 291 290 283 280 280 276 281 286 286 291 296 309 311 314 314 311 328 330 334 338 334 335 360 363 364 353 337 344 333 341 342 355 337 351 355 341 324 327 326 327 333 358 358 358 339 368 375 377 402 374 383 406 407 408 416 425 407 419 422 416 406 399 390 384 386 372 365 353 352 364 372 369 380 371 388 373 365 368 382 375 385 387 398 409 388 377 362 401 413 395 384 386 405 390 370 379 356 374 397 389 380 380 368 375 374 377 396 393 410 415 406 418 427 415 413 410 436 436 427 422 448 422 422 408 409 409 426 429 427 426 444 449 437 451 436 445 446 459 451 442 437 441 445 439 444 459 478 480 459 459 480 467 459 461 455 448 436 443 418 399 384 365 335 374 355 352 313 302 319 292 287 310 335 334 340 287 257 256 256 265 262 257 256 256 256 256 256 256 260 256 256 257 274 322 279 300 312 324 329 326 320 329 340 333 329 329 350 359 354 332 347 345 334 359 352 357 361 369 361 360 356 374 365 366 372 342 325 353 333 325 332 332 326 333 315 298 299 296 314 334 302 312 293 277 296 291 291 278 276 270 271 269 271 259 253 246 242 238 230 227 225 226 209 206 209 211 221 231 220 213 214 213 210 209 210 211 203 201 204 190 172 157 178 186 186 185 182 195 211 214 211 212 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 0 0 0 0 0 0 0 0 0 0 0 0 0
+0 0 0 0 0 0 0 0 22 5 4 87 195 184 289 679 719 235 191 334 298 476 301 5 534 699 272 132 918 1248 1441 869 1119 727 1580 1445 1238 1358 1019 1563 1430 816 131 592 310 1005 1310 1678 2003 1673 1296 1370 1351 1446 1638 1867 1662 1659 1828 1754 1521 502 1228 1768 992 277 1206 1704 1156 1715 1376 1539 1308 1638 1566 1204 1378 1481 1355 1335 1693 1787 1482 1174 1276 342 924 1126 1378 1565 1467 1333 1340 1360 1283 1146 1689 1449 976 1298 1427 1316 1089 1213 1194 1047 1146 1254 1229 1275 1525 1633 1521 1595 1660 1720 1750 1837 1750 1694 1501 1389 1243 938 1155 1458 1111 1085 908 785 548 412 345 378 431 535 737 959 1015 1023 952 962 1076 1204 1460 1602 1603 1763 1888 1730 1463 1135 1496 1683 1571 1958 1812 1756 1808 2057 1777 1697 1476 818 581 595 471 1164 1950 2241 1676 1558 1494 1916 2175 2243 1979 1725 1592 1805 1670 724 720 1549 2103 1903 1916 2190 2095 1410 2091 2052 2024 1559 686 553 1084 1361 1641 2279 2076 1968 2046 2072 2090 1615 2348 2253 1853 2082 2116 2035 2224 2082 1845 2353 2185 1562 1350 1176 1065 867 804 834 977 1513 2351 2230 1966 1798 2136 2144 2151 1774 1947 1962 1669 1928 2228 2286 2204 1880 1345 1215 1452 2057 1922 2115 2155 1800 1830 1953 1958 1827 1638 1500 1409 1467 1412 1529 1479 1310 1422 1556 1661 1555 1386 1281 1260 1266 1220 1115 1022 997 997 982 993 991 986 974 964 958 947 958 966 975 979 972 977 974 946 916 903 894 888 872 864 842 830 820 800 795 817 832 821 813 802 788 768 791 783 776 772 778 785 789 785 777 777 791 807 805 798 789 776 775 785 789 803 825 852 857 856 853 836 833 795 770 757 745 751 766 788 818 824 817 849 868 881 896 898 903 888 867 931 916 916 894 834 808 836 839 861 902 897 857 857 848 831 801 779 785 784 786 793 796 823 833 821 873 875 838 840 878 902 892 891 923 1019 1040 1033 985 989 921 877 905 892 872 908 923 903 892 917 920 912 900 906 908 910 898 885 855 832 885 854 840 816 803 806 777 755 747 736 745 738 746 726 723 722 728 730 741 753 761 772 774 761 757 726 711 709 707 701 697 702 709 719 724 724 713 703 707 708 710 714 710 714 708 712 720 724 710 690 680 699 701 708 769 791 763 818 814 765 733 736 737 746 755 726 696 673 653 637 623 623 636 645 642 623 608 597 591 584 578 575 576 574 579 577 573 568 567 572 573 579 585 585 588 592 597 597 597 597 597 597 603 603 606 613 621 628 634 636 636 638 631 630 632 628 642 673 702 723 738 738 744 752 755 773 779 784 769 746 750 744 745 726 692 676 669 664 658 652 651 636 614 610 602 596 592 587 585 585 575 560 536 530 530 524 518 510 498 489 488 481 469 459 443 437 418 391 430 435 438 447 438 428 413 391 369 363 367 392 414 422 397 387 378 383 391 379 379 388 378 369 368 370 371 366 368 376 377 384 384 384 384 379 378 375 373 382 383 378 366 362 347 343 354 354 335 308 301 294 292 292 289 286 283 271 290 286 286 281 263 256 249 249 249 243 241 239 237 237 237 237 237 237 237 237 237 231 231 232 231 230 227 227 228 227 224 230 229 231 233 236 237 237 237 239 243 247 255 260 262 266 265 262 264 269 274 275 277 286 292 292 286 281 282 281 286 290 292 302 310 312 316 323 324 335 336 340 340 333 333 351 352 364 359 370 353 336 356 355 335 357 358 365 357 345 324 316 326 334 317 333 345 355 371 344 338 338 359 385 391 397 395 400 412 416 435 422 423 400 396 439 414 401 394 386 366 369 348 371 383 370 363 373 366 361 370 372 383 396 381 399 400 403 401 410 411 405 429 408 387 407 376 384 409 393 384 383 417 422 412 399 419 419 419 395 417 421 437 428 405 402 395 415 420 416 413 411 417 421 411 409 418 435 439 449 448 451 450 449 451 451 448 451 446 446 466 466 443 446 446 452 453 469 478 455 458 457 459 454 454 450 451 451 440 423 424 413 417 396 384 349 364 335 345 298 302 304 285 267 304 299 257 256 256 256 256 266 298 304 274 256 260 256 256 261 256 256 256 256 256 256 256 276 311 307 306 305 334 333 352 331 329 329 331 331 335 356 357 356 360 360 368 370 359 382 383 392 375 364 360 383 370 364 338 329 331 329 330 331 330 318 337 331 331 344 300 309 307 300 297 289 297 298 302 288 281 296 278 271 269 270 268 268 256 241 238 237 237 234 235 220 209 242 238 237 237 237 229 216 223 225 225 224 217 213 211 210 200 168 186 209 213 220 222 221 223 238 219 214 219 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 0 0 0 0 0 0 0 0 0 0 0 0 0
+0 0 0 0 0 0 0 0 0 0 0 2 120 163 164 307 454 206 218 879 786 463 22 1 101 127 4 360 972 1369 1224 274 791 1548 1272 1716 1361 1269 703 1065 1921 1228 559 209 258 641 1266 1235 1884 1559 1464 1961 1541 776 1186 1162 2035 2015 1472 1312 690 1175 1341 969 984 693 255 1102 1356 1264 1530 1651 1600 1569 1255 1163 1451 1226 1402 1629 1523 1205 904 1188 772 468 531 552 977 964 1485 1440 1644 1455 1282 1159 1526 1293 1021 1151 1312 1277 1137 1036 1119 1032 1125 1199 1199 1352 1473 1588 1640 1664 1589 1594 1604 1668 1606 1437 1301 1277 1183 1164 1014 734 822 871 637 390 366 341 353 462 565 645 848 1004 1188 1244 1330 1320 1193 1559 1497 1493 1650 1629 1534 1401 1300 1086 1346 1331 1561 1958 1803 1592 1705 1931 1701 1363 1417 1275 611 589 584 1185 1867 2158 1758 1879 2165 1931 2346 2349 2353 1896 1488 1264 888 581 1522 1701 1485 1960 2081 2005 1501 2151 2223 1943 2096 1768 931 530 990 1557 2128 1968 2002 1622 1885 1965 1702 1423 2099 2370 2138 1542 2052 2190 1940 1978 2236 2230 2016 1625 1313 1227 1183 836 785 843 1017 1017 1955 2345 2047 2272 2326 2276 1999 1768 1735 1865 1623 1789 2011 1833 1686 2061 1387 1195 1423 1822 1975 1890 2234 2015 1766 1725 1624 1623 1556 1909 1615 1593 1398 1326 1349 1275 1415 1519 1566 1495 1621 1478 1398 1294 1192 1126 1044 999 989 976 984 993 1006 997 963 925 943 929 958 964 946 941 939 939 936 916 897 843 858 847 855 859 859 854 852 840 828 819 816 776 800 795 773 797 795 792 794 795 796 792 791 792 796 805 801 798 792 781 787 798 803 806 809 836 852 847 839 825 814 811 785 767 758 764 771 797 882 898 858 881 908 927 916 956 966 971 954 939 967 955 956 952 905 878 914 915 911 954 932 904 900 869 869 863 833 842 846 837 862 865 878 906 897 959 946 928 917 950 995 971 1016 1074 1069 1048 1022 989 993 943 902 876 865 894 909 915 903 906 926 936 917 915 917 914 910 904 858 812 795 785 779 790 788 775 768 769 762 747 733 732 733 741 755 754 748 766 760 753 761 785 771 766 758 747 751 723 715 706 701 701 701 701 703 709 714 718 709 713 714 718 719 725 726 737 732 730 718 710 691 692 692 713 726 724 719 728 726 724 729 718 724 727 740 753 740 749 749 702 670 655 641 630 635 640 651 645 626 610 595 588 582 577 573 578 579 578 574 569 566 569 573 581 586 591 591 596 597 597 597 597 599 603 604 604 604 608 615 620 622 624 622 621 623 625 620 620 634 656 702 741 762 786 773 760 756 765 770 753 747 747 740 746 755 738 704 673 663 654 652 651 646 623 611 604 597 592 585 585 582 578 565 549 527 519 526 520 509 500 494 486 480 468 453 438 432 432 413 382 383 407 422 440 435 418 369 392 395 406 423 426 422 421 421 409 401 394 395 391 384 365 356 365 374 376 376 371 370 366 366 375 380 377 377 381 377 377 378 378 375 375 371 367 359 357 362 365 352 333 311 303 295 292 292 274 290 292 292 292 287 278 267 255 249 249 244 243 237 237 237 237 237 233 237 237 237 237 236 233 235 236 232 232 231 231 230 225 228 231 231 232 233 237 237 240 243 243 249 256 260 262 268 268 263 267 268 270 280 279 277 282 289 291 286 286 286 286 291 293 307 310 315 324 329 326 332 335 337 339 336 338 328 338 362 363 372 367 376 372 351 349 362 360 349 345 341 334 324 324 326 325 319 329 323 318 345 357 361 377 376 390 398 405 415 420 436 435 417 420 403 419 428 414 410 387 375 360 368 376 368 364 362 358 374 370 366 369 366 374 382 391 401 411 397 448 469 434 434 419 405 372 378 406 409 415 392 404 411 406 418 399 395 403 419 422 409 403 402 402 410 418 406 408 419 422 411 423 423 430 441 430 425 451 433 422 424 442 445 449 456 455 451 464 454 453 464 449 442 450 451 448 449 460 447 473 454 452 453 452 454 452 458 450 432 430 430 446 449 423 403 364 355 329 309 310 300 277 303 269 256 320 300 282 257 256 256 256 256 259 322 290 279 258 256 256 256 256 256 256 256 256 256 256 270 306 302 330 349 358 340 337 343 351 348 336 330 334 332 333 354 364 387 393 364 373 402 403 378 365 365 360 359 355 348 328 339 337 326 335 325 315 317 337 331 341 341 331 306 297 292 284 293 294 301 301 295 295 295 299 298 295 271 268 267 265 260 243 240 238 237 242 251 218 239 238 234 242 237 238 235 237 238 240 235 231 233 227 219 207 193 214 211 215 226 237 241 240 244 241 237 241 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+20 0 0 0 0 0 0 0 0 28 4 0 3 24 100 65 105 174 315 157 164 78 2 29 8 2 300 1136 1131 1245 574 610 1504 1559 821 738 1216 1013 269 1008 1525 1718 1211 388 52 241 635 1131 1370 983 1141 2027 2032 1436 550 828 1474 2051 2033 1736 1329 1659 1549 1705 1312 1198 456 259 651 860 1668 1178 1417 1765 1050 702 790 1255 1300 1580 1209 1350 1144 1520 1475 435 1088 496 855 1186 1595 1593 1472 1560 1611 1318 1166 1242 1394 1406 1143 1035 1059 996 988 1086 1326 1245 1208 1373 1500 1488 1328 1413 1450 1433 1372 1282 1413 1621 1579 1245 1181 1157 933 939 686 380 346 356 377 475 661 741 991 1250 1241 1392 1412 1298 1304 1361 1274 1383 1406 1638 1662 1871 1685 1260 1180 1061 1220 1506 1674 1695 1649 1438 1951 1822 1438 1340 1277 945 624 481 612 1198 1808 1963 2106 1764 1961 1960 1905 2202 2296 2184 1665 1288 872 1028 1413 1662 1852 1745 2073 1951 1697 2238 1824 1945 1733 1240 752 532 936 1468 1827 2061 2009 1690 1670 2077 2016 1675 1472 1808 2160 1426 1772 2266 2267 2376 2101 1870 1886 1540 1232 1057 950 853 783 889 946 1634 1939 1973 1888 2240 2002 1950 1990 1489 1974 1814 1298 1785 2376 2047 1755 1797 1474 1270 1924 1603 1923 1751 2190 2102 1735 2028 1926 1746 1703 2069 1601 1546 1536 1407 1330 1263 1282 1407 1531 1485 1445 1441 1371 1400 1333 1152 1079 1000 977 966 969 973 975 968 937 944 942 940 928 916 902 913 919 929 930 910 849 884 892 874 869 868 857 853 854 844 841 835 822 813 801 799 807 802 801 804 809 812 814 813 809 813 804 800 804 808 799 785 797 802 812 827 846 845 843 830 811 796 786 771 765 763 786 829 821 865 914 936 904 1001 1027 1013 989 1026 1070 1050 1026 1007 1009 1003 1017 1029 969 950 966 1023 955 949 942 914 931 893 881 867 867 893 885 919 953 970 975 1028 1030 1057 1100 1002 1028 1028 1090 1077 1078 1073 1063 1054 1039 1007 1015 937 892 908 888 893 905 914 907 936 944 929 883 870 866 854 849 865 853 828 816 794 785 769 759 765 766 765 769 751 751 760 759 752 764 759 757 760 787 801 791 780 769 777 780 776 758 757 742 725 710 712 720 721 720 719 717 716 714 713 713 713 713 723 727 714 704 703 718 703 686 675 696 710 705 696 690 695 698 708 720 734 730 723 729 717 718 729 737 745 715 681 666 652 643 635 642 650 649 636 617 606 593 585 580 579 574 578 579 579 570 566 566 572 577 585 590 597 599 598 597 597 602 603 604 609 609 606 607 610 616 621 621 620 615 615 615 615 610 613 631 703 734 758 780 792 781 780 782 790 792 784 764 739 733 746 720 680 666 659 652 650 645 636 616 609 600 595 587 580 578 576 566 558 536 524 501 517 506 499 490 478 469 465 449 436 432 432 432 429 426 428 416 393 385 387 416 429 425 423 432 431 425 420 419 414 408 402 399 393 387 377 368 354 363 369 368 371 376 377 376 365 369 371 369 371 377 373 371 374 372 374 370 369 371 371 369 358 347 345 339 330 317 306 299 288 297 296 298 297 292 286 277 263 255 249 248 243 242 237 237 237 237 236 234 236 237 237 237 237 237 237 237 235 231 231 230 225 227 228 231 232 231 232 237 237 243 243 248 254 260 265 271 269 268 265 270 276 280 280 280 280 280 286 286 286 286 286 291 294 308 312 316 328 329 336 340 340 341 339 338 360 351 340 345 358 369 383 398 375 368 360 362 349 343 332 332 332 327 326 327 342 367 390 351 370 359 367 384 391 402 412 392 421 426 448 437 418 418 382 387 380 366 366 379 380 365 360 362 371 371 365 372 372 365 362 379 380 383 377 368 383 390 392 381 419 439 434 419 400 384 382 406 441 436 429 402 408 416 418 411 416 410 415 422 425 425 409 398 419 415 412 398 401 400 425 433 423 426 431 451 451 451 448 428 433 426 447 455 452 445 450 452 452 462 483 471 454 456 460 452 449 452 453 441 452 453 451 459 452 456 456 452 450 450 453 450 450 447 421 351 325 325 326 329 327 327 282 280 278 256 259 274 259 256 256 256 256 256 256 256 256 258 256 256 256 257 256 256 256 256 256 256 304 307 302 325 337 358 390 344 334 367 338 327 331 341 332 330 337 351 374 373 355 384 382 362 357 359 360 358 352 346 337 332 335 334 331 340 326 318 316 331 357 341 341 341 336 329 321 331 334 301 298 299 301 294 295 300 294 299 292 273 276 274 269 250 237 240 237 237 239 240 229 236 237 239 245 254 264 250 249 253 253 248 267 253 246 239 219 214 227 238 240 239 254 270 267 245 243 252 235 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+53 38 0 0 0 0 0 0 0 94 80 68 6 1 0 0 0 3 2 20 16 370 25 58 78 116 34 262 1033 1285 456 519 1333 1064 1152 488 230 295 794 1293 1555 1304 1194 916 211 110 413 902 898 864 1274 1598 1574 1299 530 1074 1759 2006 1860 1753 1189 1787 1852 1234 735 281 301 415 148 71 690 1514 834 1948 1620 693 922 1194 1450 1353 1159 1256 1338 1232 1137 812 433 666 845 1134 1368 1576 1557 1532 1313 1198 1153 1552 1529 1366 1227 1214 1008 1190 1134 1092 1266 1454 1182 1511 1566 1465 1272 1341 1333 1396 1331 1416 1505 1389 1275 1590 1654 1050 840 624 349 503 818 1104 1084 1250 1268 1409 1319 1637 1562 1524 1467 1313 1280 1199 1348 1641 1670 1964 1931 1686 1311 1190 979 1277 1165 1486 1695 1977 1806 1511 2013 1666 1285 1363 1212 1227 766 437 1091 1561 1901 2087 1520 1654 2133 1594 2050 1986 1709 2126 1741 1252 736 739 1719 1837 1621 2042 1750 2261 2266 2369 2021 1981 1966 1433 674 538 1124 1741 1985 1888 1286 1822 1911 1701 1586 1342 1650 1170 1345 1930 2162 2313 1997 1869 1911 1547 1711 1460 1208 1064 954 909 812 843 962 1824 2165 1703 1759 2167 1958 1599 1450 1365 1669 1886 1304 2058 1783 1912 2048 1978 1456 1196 1749 1670 2009 1634 2186 1996 2101 1937 1898 1815 1764 2044 1640 1534 1489 1382 1319 1297 1253 1387 1543 1536 1359 1330 1336 1348 1287 1237 1256 1045 979 968 962 952 943 935 939 938 938 952 934 915 899 918 878 920 926 866 895 910 889 878 876 864 859 859 852 848 850 843 838 836 835 829 828 820 816 818 820 842 843 831 825 821 817 813 811 805 804 804 794 811 815 826 835 853 844 826 803 793 798 792 799 800 862 897 880 874 922 962 974 967 1000 1053 1086 1042 1084 1153 1109 1082 1092 1132 1149 1143 1057 1003 1007 1036 1041 1024 995 950 945 907 899 906 930 965 969 985 1002 1074 1149 1185 1162 1174 1144 1116 1130 1121 1118 1109 1111 1092 1055 1068 1057 1037 1003 922 907 908 888 890 908 919 939 904 896 895 889 874 857 837 821 820 830 840 826 817 803 791 765 781 786 775 747 746 743 757 777 785 796 792 789 785 791 783 768 787 767 754 752 766 756 734 740 728 729 754 770 776 765 747 733 728 727 728 722 722 722 727 720 714 713 706 700 694 697 685 710 717 713 701 695 694 692 696 696 709 728 741 735 720 727 726 730 742 746 749 757 727 669 648 640 641 646 650 642 626 614 600 590 583 579 577 579 579 578 573 566 569 573 578 585 591 600 609 606 599 600 603 608 608 605 606 609 608 605 613 616 620 621 616 615 615 612 604 623 634 648 700 740 760 780 779 779 783 780 771 759 741 739 732 717 693 674 666 652 644 645 638 623 612 602 595 590 583 569 573 566 560 545 530 527 512 483 488 483 474 467 461 446 434 432 432 432 432 432 432 432 428 423 419 415 416 426 432 432 432 428 423 420 417 416 411 404 398 390 386 373 367 356 343 357 369 371 372 377 369 361 359 358 359 359 359 359 362 364 366 364 355 364 369 348 311 299 294 322 300 289 306 297 305 304 298 297 292 292 288 281 270 259 252 249 245 243 243 242 237 237 237 235 233 235 235 237 235 234 237 235 232 231 231 227 225 229 231 231 231 231 232 236 237 243 243 246 252 261 265 270 278 270 268 269 279 280 282 286 286 286 287 291 291 288 288 292 304 310 311 316 325 329 336 342 346 342 346 345 357 353 344 352 356 359 360 384 391 360 361 352 361 363 356 338 331 331 323 325 357 374 371 382 385 382 364 359 373 373 401 399 405 424 436 425 427 389 383 385 380 395 403 396 379 382 371 367 364 377 367 361 362 380 372 367 386 384 382 365 365 365 365 376 395 420 411 427 411 420 437 440 441 422 425 395 421 422 411 408 419 425 430 439 417 413 399 414 417 443 446 421 414 408 413 417 415 419 428 446 445 445 447 426 433 435 451 439 446 436 437 450 451 450 458 453 447 458 450 451 459 456 451 452 451 452 453 455 468 459 461 457 451 451 446 443 439 395 355 363 409 385 342 340 348 360 327 299 284 277 258 256 271 272 256 256 256 256 256 256 256 256 256 256 256 256 256 257 257 256 256 256 257 294 323 335 333 335 323 316 326 334 338 330 326 335 336 327 344 338 344 341 355 349 346 347 348 358 357 350 358 340 331 327 328 329 345 334 326 311 321 322 327 330 330 329 329 330 328 322 319 300 292 303 310 300 300 314 302 287 268 260 239 221 207 214 209 208 210 217 222 231 245 262 249 241 260 268 269 270 268 267 265 268 270 259 248 242 238 236 243 246 255 271 282 266 267 248 243 261 246 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+19 60 29 1 0 0 0 0 0 0 35 123 117 286 313 132 1 0 0 10 58 240 52 179 76 405 424 8 560 113 604 1088 1321 949 625 129 1035 1292 1200 1027 1236 720 1281 769 34 197 547 788 1122 843 913 1315 1353 1431 960 348 1122 1368 1273 1697 1177 1548 1477 668 895 1460 1381 925 1308 647 89 40 930 1751 1569 795 900 1069 1302 1288 1160 948 993 1284 1469 1022 249 757 1185 1245 1338 1219 1526 1691 1778 1282 1314 1657 1722 1465 1393 1167 1022 1180 1370 1139 1340 1398 1359 1496 1440 1277 1133 1197 1160 1291 1328 1212 1189 1146 1092 1459 1546 1089 1221 741 495 372 774 1335 1249 1233 1422 1694 1517 1677 1936 1782 1569 1272 1358 1342 1235 1645 1575 1642 1524 1447 1384 941 1257 996 1376 1989 2058 1946 1610 1615 1953 1469 1192 1434 1270 1516 783 559 1432 1861 2035 1701 1259 1725 2141 1962 1404 1867 1527 1105 940 1091 629 865 981 1199 1472 1609 1904 2090 1717 2144 1845 1758 1352 1167 720 570 939 834 1118 1207 1795 1775 2119 2158 1914 1643 1656 1154 2053 2190 2065 1822 1702 1547 1788 1463 1295 1137 1061 985 924 886 868 788 880 1172 1691 1564 2138 2082 1934 2079 1571 1564 1244 1500 1296 1927 1881 1633 1950 1815 1424 1444 1349 1635 1952 1562 1802 1781 2085 1777 1917 1691 1860 1908 1837 1511 1453 1373 1295 1229 1218 1403 1556 1481 1495 1297 1212 1255 1283 1155 1129 1005 976 963 959 959 953 943 942 939 940 915 925 942 944 934 897 910 926 866 901 913 901 894 881 872 869 870 868 869 854 870 880 871 862 852 844 838 832 826 836 863 871 852 844 839 836 828 822 815 816 811 803 815 823 836 846 859 854 834 814 798 815 799 807 831 904 915 931 970 978 1052 1082 1070 1027 1084 1106 1091 1106 1157 1200 1189 1211 1230 1225 1185 1125 1108 1148 1173 1199 1152 1098 1057 1032 966 986 1005 1053 1119 1127 1057 1124 1211 1211 1214 1202 1189 1183 1160 1146 1154 1145 1123 1091 1102 1107 1080 1075 1060 984 924 906 894 893 895 915 929 919 898 911 913 883 874 854 842 826 822 848 854 799 793 795 787 788 799 785 761 776 782 774 768 751 742 748 753 753 754 752 763 775 752 751 744 740 754 738 732 719 723 731 731 728 733 763 754 745 755 748 741 721 721 721 717 728 729 728 721 699 675 691 702 700 696 701 702 699 698 696 698 700 719 728 727 731 731 730 729 740 730 722 729 745 772 746 678 653 642 642 648 654 650 636 621 610 595 588 581 579 579 579 579 575 568 567 573 578 580 591 600 606 613 607 604 603 603 605 604 608 609 609 608 609 614 619 621 618 620 613 605 603 609 622 627 635 652 683 705 719 727 731 735 724 698 688 696 686 669 665 666 652 646 635 633 626 613 603 595 589 582 568 567 563 556 547 528 523 515 498 475 467 470 464 455 445 435 432 432 432 432 432 432 432 432 432 432 432 432 432 432 432 433 438 437 427 420 417 416 422 432 439 455 406 384 367 357 348 359 364 370 370 366 370 367 361 358 321 311 308 324 324 313 306 307 317 311 315 363 361 354 341 329 323 322 317 315 310 303 296 292 292 288 283 277 262 257 249 249 247 243 243 238 237 237 237 237 237 236 232 234 234 234 233 231 227 227 227 225 226 231 231 233 232 231 237 242 243 243 249 252 259 262 274 280 282 281 281 283 286 287 289 293 292 292 292 292 292 292 292 295 309 315 316 328 332 334 341 340 338 340 339 337 331 334 341 346 362 363 367 365 352 350 339 331 330 326 324 328 326 328 324 324 326 334 334 356 362 360 352 359 359 356 379 396 384 424 442 412 404 414 410 396 403 398 419 404 378 366 364 363 359 362 359 359 359 364 393 396 404 415 388 370 388 382 368 365 385 402 390 386 416 421 458 459 432 413 403 424 423 444 408 406 420 414 420 419 406 406 413 425 440 427 406 430 419 433 426 425 442 444 453 452 444 436 441 439 442 439 442 469 469 449 441 444 457 462 447 447 451 471 456 465 462 474 467 481 460 458 480 462 457 456 449 452 450 449 434 426 413 394 373 388 372 418 392 341 357 369 320 279 265 279 268 276 256 256 258 256 256 256 256 256 256 256 256 256 256 257 257 256 256 256 256 257 292 318 314 316 334 328 324 326 335 327 329 326 352 335 333 348 347 341 353 352 349 351 353 359 357 344 332 329 327 324 339 345 347 332 324 327 362 328 331 329 334 344 327 368 351 346 302 286 302 332 308 280 290 294 295 253 231 238 256 249 250 242 261 263 252 258 256 270 274 294 292 279 280 276 274 275 275 274 272 271 268 268 254 244 247 247 248 252 251 263 274 299 303 271 261 246 256 260 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+40 3 3 0 0 0 0 0 0 0 0 0 2 55 233 447 236 1 0 1 15 25 33 50 108 307 1028 337 61 638 897 731 311 71 295 574 611 885 787 1047 1102 659 59 19 343 994 1130 1111 976 1258 849 1336 1639 1391 1053 318 674 591 1284 1825 1195 1433 985 1007 1476 1642 1464 1287 980 1195 1086 593 149 272 1015 601 443 887 1260 1445 1291 1125 812 1254 1127 648 209 634 1001 1286 1439 1469 1712 1363 1627 1381 1245 1557 1540 1469 1249 991 1227 1360 1260 992 1330 1351 1303 1458 1404 1265 1052 1191 1413 1405 1626 1789 1744 1548 1313 1043 1290 1391 1022 859 698 502 352 605 1110 1324 1569 1775 1788 1725 1571 1471 1333 1301 1340 1466 1089 1280 1293 1279 1224 1199 1134 1031 948 1261 1578 1883 2019 1808 1351 1699 1789 1366 1085 1435 1429 1022 544 684 1318 1905 2109 1871 1565 1186 2004 1709 1771 1213 830 1607 1559 804 712 1598 1844 1656 1666 1481 1796 1888 1623 1803 1832 1943 1920 1105 552 611 636 1042 1782 1801 2014 1809 1646 1954 1374 1380 1464 1636 1165 2082 2200 2012 2180 1883 1333 1493 1209 1073 1022 994 942 883 862 823 804 906 1097 1816 1987 2107 1635 1464 1459 1401 1334 1305 1379 1518 1668 1850 1872 1462 1182 1630 1669 1378 1416 1469 2041 1644 1483 1542 1645 1586 1545 1774 1759 1496 1452 1350 1232 1205 1177 1304 1349 1293 1328 1282 1188 1198 1197 1073 1039 1000 999 996 995 989 985 966 947 932 922 932 979 1015 997 967 937 865 914 913 867 907 912 911 911 906 906 900 899 899 899 898 882 871 878 878 871 868 870 864 859 865 873 877 866 856 843 838 835 839 832 827 830 829 824 842 849 854 851 838 824 817 827 810 803 826 899 942 1003 1026 1064 1113 1147 1136 1126 1090 1114 1109 1199 1223 1239 1241 1315 1394 1374 1288 1272 1232 1304 1307 1277 1226 1207 1166 1116 1066 1106 1151 1266 1257 1211 1157 1131 1185 1175 1208 1194 1173 1170 1144 1134 1141 1153 1120 1124 1129 1125 1102 1070 1001 931 906 904 905 905 906 916 922 913 919 931 937 911 900 898 881 863 863 849 856 838 833 814 809 795 783 770 784 801 803 790 776 784 787 780 772 762 742 733 735 735 728 726 725 727 728 734 726 725 746 743 738 749 765 748 748 739 755 756 752 767 752 730 719 720 724 720 702 669 685 692 688 699 725 726 719 724 734 741 728 733 740 719 704 724 727 725 718 720 713 715 725 736 740 730 703 673 657 647 646 646 652 651 644 630 613 603 590 586 580 576 578 580 574 567 563 567 574 584 591 600 608 609 609 606 603 603 604 606 606 608 608 609 613 609 616 617 619 615 607 603 598 603 609 617 620 613 616 621 623 630 633 638 643 650 658 658 650 646 644 648 646 631 622 627 619 603 594 591 582 566 563 561 554 532 529 514 505 499 490 470 455 465 454 444 436 432 432 432 432 432 432 432 432 432 429 426 421 426 431 432 432 433 438 443 445 442 453 457 468 461 443 409 395 379 366 362 349 338 347 358 359 360 360 357 351 331 359 355 360 363 360 353 339 347 364 366 364 358 353 345 331 326 322 316 316 311 310 301 292 292 289 283 279 264 257 255 249 249 244 242 237 237 237 237 237 237 237 237 235 231 231 231 231 231 231 231 226 229 231 235 233 231 233 237 243 243 245 249 253 260 272 280 281 286 287 291 289 294 299 303 308 309 308 307 303 298 293 295 298 310 311 316 325 329 328 329 330 334 334 329 329 325 336 337 334 335 336 348 340 329 329 324 327 324 325 326 337 330 326 323 327 325 326 324 334 332 342 356 359 369 388 393 411 395 430 431 426 422 436 399 416 417 410 400 404 378 396 375 360 359 361 370 370 362 359 380 388 404 402 387 377 391 375 381 370 375 367 369 394 418 429 437 443 439 432 420 414 420 419 410 432 424 420 428 417 420 426 428 443 437 416 398 409 416 436 429 423 450 449 434 438 451 460 439 450 461 471 462 469 476 475 472 460 476 463 481 488 489 465 477 475 489 483 478 479 486 478 467 458 458 461 450 449 426 422 420 429 418 399 381 356 420 421 384 324 308 339 339 294 270 314 341 296 263 258 258 256 256 256 256 281 256 256 256 265 301 271 256 256 256 256 256 257 281 298 306 319 357 348 336 325 327 325 328 329 328 347 360 368 362 374 354 368 370 368 360 362 357 359 339 327 324 329 344 347 354 338 319 350 344 325 327 328 328 329 324 330 338 307 308 311 323 296 273 265 275 275 279 276 267 273 292 279 265 265 272 277 278 259 281 297 299 331 319 303 283 293 295 294 295 287 285 274 273 264 261 260 268 275 268 267 255 257 263 268 288 274 259 245 243 252 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+130 74 17 0 0 0 0 0 0 0 0 0 0 0 70 354 504 347 4 0 0 0 0 10 137 294 810 766 37 2 130 790 953 1137 1014 1039 919 1110 456 1067 755 956 89 127 824 1123 1408 1159 757 1081 1083 1186 1314 1341 1181 580 159 513 1222 1253 1533 834 743 1455 1333 1441 1117 1068 1050 1103 1251 1351 760 92 10 548 261 782 978 964 1462 1415 1195 1501 1405 1103 463 773 1029 1229 1215 1363 1546 1632 1510 1222 1321 1603 1439 1521 1358 1032 1107 1355 1511 1263 989 1240 1042 1382 1281 959 1279 1442 1479 1641 1719 1717 1745 1612 1526 1290 1064 930 930 790 883 632 395 349 708 1198 1456 1757 1600 1891 1697 1403 1198 1204 1253 1576 1117 1204 1079 1340 1323 1260 1084 889 1204 1594 1745 1804 1981 1519 1376 1867 1618 1122 1097 1363 1644 1108 526 687 1292 1768 2054 1659 1414 1354 1323 925 1116 785 1382 1601 890 655 1286 1717 1944 1983 2055 1857 1256 1426 1266 1048 927 794 617 592 581 534 738 610 1151 1518 2085 1997 1681 1712 1689 2015 1952 1798 1260 1069 1168 1209 1457 1034 1025 1044 1010 943 948 942 937 904 867 809 780 803 911 1377 1897 2134 2188 1842 1414 1281 1214 1625 1731 1976 1509 1790 1994 1166 1224 1816 1688 1388 1757 1509 1519 1566 1434 1461 1415 1371 1694 1584 1592 1397 1380 1281 1198 1173 1158 1196 1199 1196 1203 1197 1169 1089 1064 1026 1028 1017 1029 1041 1075 1092 1061 1014 972 940 951 953 966 1019 1020 977 957 925 866 900 923 922 919 920 922 925 922 923 926 925 926 924 916 904 880 887 901 898 898 894 890 889 884 881 877 870 859 852 847 842 843 841 839 836 838 828 853 859 852 836 832 831 833 808 819 838 881 937 1003 966 994 1007 1027 1057 1104 1100 1214 1272 1148 1351 1428 1405 1431 1418 1412 1399 1355 1348 1249 1207 1275 1329 1249 1188 1172 1102 1093 1157 1255 1220 1219 1193 1110 1144 1129 1165 1165 1134 1125 1115 1153 1164 1147 1094 1118 1114 1128 1091 1045 934 915 924 915 909 908 913 948 942 954 964 960 919 915 917 901 874 857 851 841 837 872 904 894 890 893 854 815 798 811 817 801 810 783 770 765 758 750 745 741 734 733 735 745 754 775 759 735 731 721 730 757 746 734 756 768 754 739 761 760 766 755 748 745 745 746 732 716 703 677 697 692 713 710 717 735 746 751 756 777 765 759 748 720 711 726 719 715 713 713 713 713 717 724 725 720 715 701 677 666 652 646 643 648 654 645 636 618 604 593 587 580 579 576 573 573 566 568 562 573 580 591 596 603 603 604 608 603 597 600 603 603 609 609 605 603 609 615 607 605 611 614 614 607 599 595 593 593 598 602 610 620 626 618 620 624 629 636 630 630 624 631 638 624 621 621 621 604 591 591 581 565 560 560 545 531 521 513 498 498 491 478 458 456 455 442 442 437 432 432 432 432 432 432 426 426 427 424 432 431 431 429 423 422 436 447 454 457 457 466 469 458 430 416 405 393 378 365 355 355 356 341 333 332 344 348 342 369 368 365 368 368 365 346 362 373 371 371 375 378 375 385 374 349 330 317 316 312 310 305 297 292 289 283 280 269 261 256 251 249 246 243 238 237 237 237 237 237 237 237 237 236 232 234 231 231 231 231 226 229 231 236 237 237 237 239 243 246 249 253 256 261 274 280 282 287 293 301 298 303 310 310 313 315 313 310 310 310 310 310 309 310 310 311 311 312 316 316 324 329 329 329 329 327 324 324 324 327 326 324 336 351 338 330 329 326 335 335 331 325 326 324 324 324 329 329 353 345 356 370 366 378 370 365 377 381 354 388 397 375 381 389 416 410 384 402 427 395 397 367 366 370 368 388 392 381 365 373 390 397 396 378 400 399 397 402 385 366 367 370 384 407 416 411 416 412 408 412 426 426 420 422 421 418 441 461 425 444 430 415 427 426 429 429 444 430 413 423 444 453 449 425 436 448 458 450 452 455 460 475 481 480 468 476 472 462 483 481 480 483 475 474 483 485 488 483 485 482 479 466 473 476 472 460 448 438 435 430 422 420 403 362 389 401 433 376 299 272 305 334 296 272 296 346 319 282 256 256 260 267 257 271 286 314 256 258 283 308 259 258 256 256 256 256 258 279 320 352 322 347 355 362 369 370 374 376 374 374 386 392 388 366 372 388 391 391 390 397 376 358 344 337 339 336 346 352 327 311 331 330 334 331 330 327 330 323 325 328 329 328 334 347 337 333 324 315 312 316 343 344 343 303 312 295 293 304 304 310 327 316 312 329 339 321 330 322 312 304 314 308 303 301 320 301 287 275 276 292 325 329 324 275 273 268 262 268 273 300 291 273 270 260 254 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+247 203 92 6 16 15 0 0 0 0 0 0 0 0 0 1 91 376 585 108 0 0 0 0 90 267 584 728 268 48 758 797 1117 965 1140 856 611 1040 549 488 227 89 3 375 884 968 1139 989 1028 559 875 977 1459 1410 1096 967 656 197 1110 1375 1206 1189 818 684 820 1272 1383 820 906 897 1652 1435 1192 719 56 14 51 974 930 1051 1226 772 1417 1401 1260 745 377 678 850 1446 1461 1346 1271 1632 1135 1567 1499 1437 1321 1317 1423 1109 1045 1321 1530 1431 1190 972 975 1248 1003 1208 1493 1515 1686 1792 1844 1665 1604 1766 1555 1578 1368 1359 1158 870 731 693 531 350 462 992 1411 1464 1649 1816 1750 1477 1206 1115 1339 1301 1103 990 1278 1361 1252 1265 1153 841 1411 1614 1866 1915 1757 1280 1981 1533 1197 1167 1050 1371 1693 1050 570 571 1249 1942 1923 1811 1625 1951 1787 1558 804 1178 1162 661 938 1279 1155 1432 1767 1407 1733 1749 978 596 923 1169 1327 1301 1408 1647 1695 851 530 780 1703 1926 2039 2096 1654 2055 2136 2093 2127 2011 1483 1769 1788 1765 1315 1780 1897 1799 1519 1087 1018 995 978 866 890 929 843 806 789 881 1021 1807 2188 1547 1124 1373 1356 1635 1753 2262 2102 1415 1195 1236 1862 1962 1747 1426 1560 1628 1978 1970 2071 1731 1771 1488 1732 1309 1465 1408 1228 1220 1216 1199 1176 1162 1166 1147 1100 1113 1105 1072 1051 1042 1046 1081 1125 1114 1118 1106 1080 1030 990 967 953 979 1003 1006 1004 986 959 947 933 884 901 923 925 921 922 927 932 935 932 934 943 943 942 937 934 927 916 913 903 888 881 893 883 867 863 873 879 875 871 862 856 854 852 847 842 843 824 835 847 854 849 842 817 834 824 834 860 895 913 943 984 1024 1056 1057 1063 1095 1139 1195 1244 1175 1289 1300 1326 1383 1391 1405 1387 1365 1321 1289 1193 1184 1244 1268 1228 1148 1092 1082 1127 1152 1150 1122 1083 1092 1105 1117 1106 1083 1089 1156 1147 1128 1138 1108 1061 1089 1093 984 939 947 946 932 925 924 938 981 995 990 1000 978 949 936 944 923 911 904 897 857 845 852 896 915 917 943 956 934 881 860 847 808 787 784 775 781 785 757 745 750 740 762 744 740 763 789 764 735 725 725 727 726 742 756 730 753 762 760 758 764 774 781 785 771 753 737 738 750 756 730 702 707 710 702 702 708 752 754 769 763 758 768 786 761 734 730 729 733 730 726 722 721 720 725 725 727 731 734 722 702 680 664 655 646 644 646 646 645 640 625 610 597 590 584 579 578 568 569 572 560 560 571 575 587 595 603 609 606 603 603 597 597 597 604 604 608 608 603 608 608 604 603 602 604 606 605 610 609 608 605 600 594 590 590 592 598 606 609 612 613 613 620 620 621 616 612 608 615 603 589 582 580 563 558 556 544 530 512 510 499 489 488 477 467 451 446 441 438 435 434 432 432 432 432 432 432 432 432 432 429 424 427 431 435 433 431 424 432 446 466 469 452 439 434 419 403 386 376 383 378 374 370 368 365 364 370 367 369 371 373 371 371 371 366 373 376 383 389 389 425 437 439 436 458 435 399 345 318 312 310 305 299 292 287 283 280 268 262 256 255 249 249 243 242 237 237 237 237 232 235 237 237 237 236 235 234 232 232 227 227 231 231 236 237 237 242 247 249 250 256 258 262 262 270 279 285 291 295 308 310 310 313 316 319 320 316 311 315 316 316 322 322 323 324 327 317 326 328 329 329 329 329 329 329 334 329 332 325 323 323 326 332 347 335 341 345 340 336 326 323 327 323 325 323 325 336 343 355 358 361 358 363 362 361 375 380 354 382 371 384 375 365 392 392 399 398 433 416 402 397 396 397 412 383 419 420 366 362 368 380 401 392 411 416 426 426 401 372 380 390 399 379 375 375 384 397 416 414 441 456 442 446 429 436 422 431 447 437 416 404 412 437 428 427 425 437 449 439 431 426 424 450 422 441 467 477 455 448 459 477 473 485 476 482 491 492 472 470 471 473 485 482 480 484 492 493 477 489 483 481 476 481 471 459 465 453 452 456 453 434 423 379 370 382 405 420 380 357 310 294 279 278 274 256 289 299 280 256 256 272 279 270 258 264 268 256 256 256 256 256 257 257 256 256 268 310 295 338 324 339 377 375 386 404 417 403 386 411 420 396 383 373 361 361 376 390 391 391 397 381 363 356 348 352 348 354 342 315 334 369 382 374 335 338 318 332 331 329 359 366 340 330 342 356 357 337 332 332 351 365 365 322 327 335 338 306 331 326 328 336 323 330 344 335 332 334 329 329 332 331 335 316 299 299 301 298 281 320 332 343 322 300 296 275 271 266 267 270 292 280 272 271 268 268 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 ...
[truncated message content] |
|
From: <js...@us...> - 2009-01-26 20:27:22
|
Revision: 6834
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=6834&view=rev
Author: jswhit
Date: 2009-01-26 20:27:17 +0000 (Mon, 26 Jan 2009)
Log Message:
-----------
specify calendar
Modified Paths:
--------------
trunk/toolkits/basemap/examples/plotsst.py
trunk/toolkits/basemap/examples/pnganim.py
Modified: trunk/toolkits/basemap/examples/plotsst.py
===================================================================
--- trunk/toolkits/basemap/examples/plotsst.py 2009-01-26 20:01:58 UTC (rev 6833)
+++ trunk/toolkits/basemap/examples/plotsst.py 2009-01-26 20:27:17 UTC (rev 6834)
@@ -15,8 +15,8 @@
dataset = NetCDFFile('http://nomads.ncdc.noaa.gov/thredds/dodsC/oisst/totalAagg')
# find index of desired time.
time = dataset.variables['time']
-nt = date2index(date, time)
-print num2date(time[nt],time.units)
+nt = date2index(date, time, calendar='standard')
+print num2date(time[nt],time.units, calendar='standard')
# read sst. Will automatically create a masked array using
# missing_value variable attribute.
sst = dataset.variables['sst'][nt]
Modified: trunk/toolkits/basemap/examples/pnganim.py
===================================================================
--- trunk/toolkits/basemap/examples/pnganim.py 2009-01-26 20:01:58 UTC (rev 6833)
+++ trunk/toolkits/basemap/examples/pnganim.py 2009-01-26 20:27:17 UTC (rev 6834)
@@ -9,15 +9,13 @@
import matplotlib.mlab as mlab
import numpy.ma as ma
import datetime, sys, time, subprocess
-from mpl_toolkits.basemap import Basemap, shiftgrid, NetCDFFile, num2date
+from mpl_toolkits.basemap import Basemap, shiftgrid, NetCDFFile, date2index, num2date
# times for March 1993 'storm of the century'
-YYYYMMDDHH1 = '1993031000'
-YYYYMMDDHH2 = '1993031700'
+date1 = datetime.datetime(1993,3,10,0)
+date2 = datetime.datetime(1993,3,17,0)
+print date1, date2
-YYYY = YYYYMMDDHH1[0:4]
-if YYYY != YYYYMMDDHH2[0:4]:
- raise ValueError,'dates must be in same year'
# set OpenDAP server URL.
URL="http://nomad1.ncep.noaa.gov:9090/dods/reanalyses/reanalysis-2/6hr/pgb/pgb"
@@ -32,18 +30,10 @@
latitudes = data.variables['lat'][:]
longitudes = data.variables['lon'][:].tolist()
times = data.variables['time']
-# convert numeric time values to datetime objects.
-fdates = num2date(times[:],units=times.units,calendar='standard')
-# put times in YYYYMMDDHH format.
-dates = [fdate.strftime('%Y%m%d%H') for fdate in fdates]
-if YYYYMMDDHH1 not in dates or YYYYMMDDHH2 not in dates:
- raise ValueError, 'date1 or date2 not a valid date (must be in form YYYYMMDDHH, where HH is 00,06,12 or 18)'
-# find indices bounding desired times.
-ntime1 = dates.index(YYYYMMDDHH1)
-ntime2 = dates.index(YYYYMMDDHH2)
+ntime1 = date2index(date1,times,calendar='standard')
+ntime2 = date2index(date2,times,calendar='standard')
print 'ntime1,ntime2:',ntime1,ntime2
-if ntime1 >= ntime2:
- raise ValueError,'date2 must be greater than date1'
+print num2date(times[ntime1],times.units,calendar='standard'), num2date(times[ntime2],times.units,calendar='standard')
# get sea level pressure and 10-m wind data.
slpdata = data.variables['presmsl']
udata = data.variables['ugrdprs']
@@ -52,7 +42,7 @@
slpin = 0.01*slpdata[ntime1:ntime2+1,:,:]
uin = udata[ntime1:ntime2+1,0,:,:]
vin = vdata[ntime1:ntime2+1,0,:,:]
-datelabels = dates[ntime1:ntime2+1]
+dates = num2date(times[ntime1:ntime2+1], times.units, calendar='standard')
# add cyclic points
slp = np.zeros((slpin.shape[0],slpin.shape[1],slpin.shape[2]+1),np.float64)
slp[:,:,0:-1] = slpin; slp[:,:,-1] = slpin[:,:,0]
@@ -68,13 +58,12 @@
print uin.min(), uin.max()
print vin.min(), vin.max()
print 'dates'
-print datelabels
+print dates
# make orthographic basemaplt.
m = Basemap(resolution='c',projection='ortho',lat_0=60.,lon_0=-60.)
plt.ion() # interactive mode on.
uin = udata[ntime1:ntime2+1,0,:,:]
vin = vdata[ntime1:ntime2+1,0,:,:]
-datelabels = dates[ntime1:ntime2+1]
# make orthographic basemaplt.
m = Basemap(resolution='c',projection='ortho',lat_0=60.,lon_0=-60.)
plt.ion() # interactive mode on.
@@ -94,7 +83,7 @@
l, b, w, h = pos.bounds
# loop over times, make contour plots, draw coastlines,
# parallels, meridians and title.
-for nt,date in enumerate(datelabels[1:]):
+for nt,date in enumerate(dates):
CS = m.contour(x,y,slp[nt,:,:],clevs,linewidths=0.5,colors='k',animated=True)
CS = m.contourf(x,y,slp[nt,:,:],clevs,cmap=plt.cm.RdBu_r,animated=True)
# plot wind vectors on lat/lon grid.
@@ -117,7 +106,7 @@
m.drawcoastlines(linewidth=1.5)
m.drawparallels(parallels)
m.drawmeridians(meridians)
- plt.title('SLP and Wind Vectors '+date)
+ plt.title('SLP and Wind Vectors '+str(date))
if nt == 0: # plot colorbar on a separate axes (only for first frame)
cax = plt.axes([l+w-0.05, b, 0.03, h]) # setup colorbar axes
fig.colorbar(CS,drawedges=True, cax=cax) # draw colorbar
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2009-06-23 13:15:40
|
Revision: 7231
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=7231&view=rev
Author: jswhit
Date: 2009-06-23 12:39:21 +0000 (Tue, 23 Jun 2009)
Log Message:
-----------
change mlab.load to np.loadtxt
Modified Paths:
--------------
trunk/toolkits/basemap/examples/contour_demo.py
trunk/toolkits/basemap/examples/maskoceans.py
trunk/toolkits/basemap/examples/panelplot.py
trunk/toolkits/basemap/examples/plotmap.py
trunk/toolkits/basemap/examples/plotmap_masked.py
trunk/toolkits/basemap/examples/plotmap_shaded.py
trunk/toolkits/basemap/examples/polarmaps.py
trunk/toolkits/basemap/examples/simpletest.py
trunk/toolkits/basemap/examples/simpletest_oo.py
trunk/toolkits/basemap/examples/test.py
Modified: trunk/toolkits/basemap/examples/contour_demo.py
===================================================================
--- trunk/toolkits/basemap/examples/contour_demo.py 2009-06-21 18:53:49 UTC (rev 7230)
+++ trunk/toolkits/basemap/examples/contour_demo.py 2009-06-23 12:39:21 UTC (rev 7231)
@@ -1,14 +1,13 @@
from mpl_toolkits.basemap import Basemap, shiftgrid
import numpy as np
import matplotlib.pyplot as plt
-import matplotlib.mlab as mlab
# examples of filled contour plots on map projections.
# read in data on lat/lon grid.
-hgt = mlab.load('500hgtdata.gz')
-lons = mlab.load('500hgtlons.gz')
-lats = mlab.load('500hgtlats.gz')
+hgt = np.loadtxt('500hgtdata.gz')
+lons = np.loadtxt('500hgtlons.gz')
+lats = np.loadtxt('500hgtlats.gz')
# shift data so lons go from -180 to 180 instead of 0 to 360.
hgt,lons = shiftgrid(180.,hgt,lons,start=False)
lons, lats = np.meshgrid(lons, lats)
Modified: trunk/toolkits/basemap/examples/maskoceans.py
===================================================================
--- trunk/toolkits/basemap/examples/maskoceans.py 2009-06-21 18:53:49 UTC (rev 7230)
+++ trunk/toolkits/basemap/examples/maskoceans.py 2009-06-23 12:39:21 UTC (rev 7231)
@@ -1,13 +1,12 @@
from mpl_toolkits.basemap import Basemap, shiftgrid, maskoceans, interp
import numpy as np
import matplotlib.pyplot as plt
-import matplotlib.mlab as mlab
# example showing how to mask out 'wet' areas on a contour or pcolor plot.
-topodatin = mlab.load('etopo20data.gz')
-lonsin = mlab.load('etopo20lons.gz')
-latsin = mlab.load('etopo20lats.gz')
+topodatin = np.loadtxt('etopo20data.gz')
+lonsin = np.loadtxt('etopo20lons.gz')
+latsin = np.loadtxt('etopo20lats.gz')
# shift data so lons go from -180 to 180 instead of 20 to 380.
topoin,lons1 = shiftgrid(180.,topodatin,lonsin,start=False)
Modified: trunk/toolkits/basemap/examples/panelplot.py
===================================================================
--- trunk/toolkits/basemap/examples/panelplot.py 2009-06-21 18:53:49 UTC (rev 7230)
+++ trunk/toolkits/basemap/examples/panelplot.py 2009-06-23 12:39:21 UTC (rev 7231)
@@ -3,13 +3,12 @@
from matplotlib.ticker import MultipleLocator
import numpy as np
import matplotlib.pyplot as plt
-import matplotlib.mlab as mlab
# read in data on lat/lon grid.
-hgt = mlab.load('500hgtdata.gz')
-lons = mlab.load('500hgtlons.gz')
-lats = mlab.load('500hgtlats.gz')
+hgt = np.loadtxt('500hgtdata.gz')
+lons = np.loadtxt('500hgtlons.gz')
+lats = np.loadtxt('500hgtlats.gz')
lons, lats = np.meshgrid(lons, lats)
# Example to show how to make multi-panel plots.
Modified: trunk/toolkits/basemap/examples/plotmap.py
===================================================================
--- trunk/toolkits/basemap/examples/plotmap.py 2009-06-21 18:53:49 UTC (rev 7230)
+++ trunk/toolkits/basemap/examples/plotmap.py 2009-06-23 12:39:21 UTC (rev 7231)
@@ -7,13 +7,12 @@
from mpl_toolkits.basemap import Basemap, shiftgrid
import numpy as np
import matplotlib.pyplot as plt
-import matplotlib.mlab as mlab
# read in topo data (on a regular lat/lon grid)
# longitudes go from 20 to 380.
-topoin = mlab.load('etopo20data.gz')
-lons = mlab.load('etopo20lons.gz')
-lats = mlab.load('etopo20lats.gz')
+topoin = np.loadtxt('etopo20data.gz')
+lons = np.loadtxt('etopo20lons.gz')
+lats = np.loadtxt('etopo20lats.gz')
# shift data so lons go from -180 to 180 instead of 20 to 380.
topoin,lons = shiftgrid(180.,topoin,lons,start=False)
@@ -30,18 +29,19 @@
fig=plt.figure(figsize=(8,8))
# add an axes, leaving room for colorbar on the right.
ax = fig.add_axes([0.1,0.1,0.7,0.7])
+# associate this axes with the Basemap instance.
+m.ax = ax
# plot image over map with imshow.
im = m.imshow(topodat,plt.cm.jet)
# setup colorbar axes instance.
pos = ax.get_position()
l, b, w, h = pos.bounds
cax = plt.axes([l+w+0.075, b, 0.05, h])
-plt.colorbar(cax=cax) # draw colorbar
-plt.axes(ax) # make the original axes current again
+plt.colorbar(im,cax=cax) # draw colorbar
# plot blue dot on boulder, colorado and label it as such.
xpt,ypt = m(-104.237,40.125)
m.plot([xpt],[ypt],'bo')
-plt.text(xpt+100000,ypt+100000,'Boulder')
+ax.text(xpt+100000,ypt+100000,'Boulder')
# draw coastlines and political boundaries.
m.drawcoastlines()
m.drawcountries()
@@ -53,6 +53,5 @@
meridians = np.arange(10.,360.,30.)
m.drawmeridians(meridians,labels=[1,1,0,1])
# set title.
-plt.title('ETOPO Topography - Lambert Conformal Conic')
-#plt.savefig('plotmap.pdf')
+ax.set_title('ETOPO Topography - Lambert Conformal Conic')
plt.show()
Modified: trunk/toolkits/basemap/examples/plotmap_masked.py
===================================================================
--- trunk/toolkits/basemap/examples/plotmap_masked.py 2009-06-21 18:53:49 UTC (rev 7230)
+++ trunk/toolkits/basemap/examples/plotmap_masked.py 2009-06-23 12:39:21 UTC (rev 7231)
@@ -9,14 +9,13 @@
import numpy.ma as ma
import numpy as np
import matplotlib.pyplot as plt
-import matplotlib.mlab as mlab
import matplotlib.colors as colors
# read in topo data (on a regular lat/lon grid)
# longitudes go from 20 to 380.
-topoin = mlab.load('etopo20data.gz')
-lonsin = mlab.load('etopo20lons.gz')
-latsin = mlab.load('etopo20lats.gz')
+topoin = np.loadtxt('etopo20data.gz')
+lonsin = np.loadtxt('etopo20lons.gz')
+latsin = np.loadtxt('etopo20lats.gz')
# shift data so lons go from -180 to 180 instead of 20 to 380.
topoin,lonsin = shiftgrid(180.,topoin,lonsin,start=False)
@@ -33,6 +32,8 @@
fig=plt.figure(figsize=(8,8))
# add an axes, leaving room for colorbar on the right.
ax = fig.add_axes([0.1,0.1,0.7,0.7])
+# associate this axes with the Basemap instance.
+m.ax = ax
# make topodat a masked array, masking values lower than sea level.
topodat = np.where(topodat < 0.,1.e10,topodat)
topodatm = ma.masked_values(topodat, 1.e10)
@@ -44,12 +45,11 @@
pos = ax.get_position()
l, b, w, h = pos.bounds
cax = plt.axes([l+w+0.075, b, 0.05, h])
-plt.colorbar(cax=cax) # draw colorbar
-plt.axes(ax) # make the original axes current again
+plt.colorbar(im,cax=cax) # draw colorbar
# plot blue dot on boulder, colorado and label it as such.
xpt,ypt = m(-104.237,40.125)
m.plot([xpt],[ypt],'bo')
-plt.text(xpt+100000,ypt+100000,'Boulder')
+ax.text(xpt+100000,ypt+100000,'Boulder')
# draw coastlines and political boundaries.
m.drawcoastlines()
m.drawcountries()
@@ -61,5 +61,5 @@
meridians = np.arange(10.,360.,30.)
m.drawmeridians(meridians,labels=[1,1,0,1])
# set title.
-plt.title('Masked ETOPO Topography - via imshow')
+ax.set_title('Masked ETOPO Topography - via imshow')
plt.show()
Modified: trunk/toolkits/basemap/examples/plotmap_shaded.py
===================================================================
--- trunk/toolkits/basemap/examples/plotmap_shaded.py 2009-06-21 18:53:49 UTC (rev 7230)
+++ trunk/toolkits/basemap/examples/plotmap_shaded.py 2009-06-23 12:39:21 UTC (rev 7231)
@@ -6,7 +6,6 @@
from mpl_toolkits.basemap import Basemap, shiftgrid
import numpy as np
import matplotlib.pyplot as plt
-import matplotlib.mlab as mlab
try:
from matplotlib.colors import LightSource
except ImportError:
@@ -15,9 +14,9 @@
# read in topo data (on a regular lat/lon grid)
# longitudes go from 20 to 380.
-topoin = mlab.load('etopo20data.gz')
-lons = mlab.load('etopo20lons.gz')
-lats = mlab.load('etopo20lats.gz')
+topoin = np.loadtxt('etopo20data.gz')
+lons = np.loadtxt('etopo20lons.gz')
+lats = np.loadtxt('etopo20lats.gz')
# shift data so lons go from -180 to 180 instead of 20 to 380.
topoin,lons = shiftgrid(180.,topoin,lons,start=False)
Modified: trunk/toolkits/basemap/examples/polarmaps.py
===================================================================
--- trunk/toolkits/basemap/examples/polarmaps.py 2009-06-21 18:53:49 UTC (rev 7230)
+++ trunk/toolkits/basemap/examples/polarmaps.py 2009-06-23 12:39:21 UTC (rev 7231)
@@ -8,13 +8,12 @@
from mpl_toolkits.basemap import Basemap
import numpy as np
import matplotlib.pyplot as plt
-import matplotlib.mlab as mlab
# read in topo data (on a regular lat/lon grid)
# longitudes go from 20 to 380.
-etopo = mlab.load('etopo20data.gz')
-lons = mlab.load('etopo20lons.gz')
-lats = mlab.load('etopo20lats.gz')
+etopo = np.loadtxt('etopo20data.gz')
+lons = np.loadtxt('etopo20lons.gz')
+lats = np.loadtxt('etopo20lats.gz')
print 'min/max etopo20 data:'
print etopo.min(),etopo.max()
Modified: trunk/toolkits/basemap/examples/simpletest.py
===================================================================
--- trunk/toolkits/basemap/examples/simpletest.py 2009-06-21 18:53:49 UTC (rev 7230)
+++ trunk/toolkits/basemap/examples/simpletest.py 2009-06-23 12:39:21 UTC (rev 7231)
@@ -1,11 +1,10 @@
from mpl_toolkits.basemap import Basemap
import numpy as np
import matplotlib.pyplot as plt
-import matplotlib.mlab as mlab
# read in topo data (on a regular lat/lon grid)
-etopo=mlab.load('etopo20data.gz')
-lons=mlab.load('etopo20lons.gz')
-lats=mlab.load('etopo20lats.gz')
+etopo=np.loadtxt('etopo20data.gz')
+lons=np.loadtxt('etopo20lons.gz')
+lats=np.loadtxt('etopo20lats.gz')
# create Basemap instance for Robinson projection.
m = Basemap(projection='robin',lon_0=0.5*(lons[0]+lons[-1]))
# make filled contour plot.
Modified: trunk/toolkits/basemap/examples/simpletest_oo.py
===================================================================
--- trunk/toolkits/basemap/examples/simpletest_oo.py 2009-06-21 18:53:49 UTC (rev 7230)
+++ trunk/toolkits/basemap/examples/simpletest_oo.py 2009-06-23 12:39:21 UTC (rev 7231)
@@ -8,15 +8,13 @@
from mpl_toolkits.basemap import Basemap
from matplotlib.figure import Figure
import numpy as np
-import matplotlib.mlab as mlab
import matplotlib.cm as cm
-from matplotlib.mlab import load
# read in topo data (on a regular lat/lon grid)
# longitudes go from 20 to 380.
-etopo = mlab.load('etopo20data.gz')
-lons = mlab.load('etopo20lons.gz')
-lats = mlab.load('etopo20lats.gz')
+etopo = np.loadtxt('etopo20data.gz')
+lons = np.loadtxt('etopo20lons.gz')
+lats = np.loadtxt('etopo20lats.gz')
# create figure.
fig = Figure()
canvas = FigureCanvas(fig)
Modified: trunk/toolkits/basemap/examples/test.py
===================================================================
--- trunk/toolkits/basemap/examples/test.py 2009-06-21 18:53:49 UTC (rev 7230)
+++ trunk/toolkits/basemap/examples/test.py 2009-06-23 12:39:21 UTC (rev 7231)
@@ -6,14 +6,13 @@
from mpl_toolkits.basemap import Basemap, shiftgrid
import numpy as np
import matplotlib.pyplot as plt
-import matplotlib.mlab as mlab
import matplotlib.colors as colors
# read in topo data (on a regular lat/lon grid)
# longitudes go from 20 to 380.
-topodatin = mlab.load('etopo20data.gz')
-lonsin = mlab.load('etopo20lons.gz')
-latsin = mlab.load('etopo20lats.gz')
+topodatin = np.loadtxt('etopo20data.gz')
+lonsin = np.loadtxt('etopo20lons.gz')
+latsin = np.loadtxt('etopo20lats.gz')
# shift data so lons go from -180 to 180 instead of 20 to 380.
topoin,lons = shiftgrid(180.,topodatin,lonsin,start=False)
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2009-08-24 19:03:46
|
Revision: 7555
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=7555&view=rev
Author: jswhit
Date: 2009-08-24 19:03:32 +0000 (Mon, 24 Aug 2009)
Log Message:
-----------
rename example
Added Paths:
-----------
trunk/toolkits/basemap/examples/daynight.py
Removed Paths:
-------------
trunk/toolkits/basemap/examples/terminator.py
Copied: trunk/toolkits/basemap/examples/daynight.py (from rev 7554, trunk/toolkits/basemap/examples/terminator.py)
===================================================================
--- trunk/toolkits/basemap/examples/daynight.py (rev 0)
+++ trunk/toolkits/basemap/examples/daynight.py 2009-08-24 19:03:32 UTC (rev 7555)
@@ -0,0 +1,75 @@
+import numpy as np
+from mpl_toolkits.basemap import Basemap, netcdftime
+import matplotlib.pyplot as plt
+from datetime import datetime
+
+def epem(date):
+ """
+ input: date - datetime object (assumed UTC)
+ ouput: gha - Greenwich hour angle, the angle between the Greenwich
+ meridian and the meridian containing the subsolar point.
+ dec - solar declination.
+ """
+ dg2rad = np.pi/180.
+ rad2dg = 1./dg2rad
+ # compute julian day from UTC datetime object.
+ jday = netcdftime.JulianDayFromDate(date)
+ jd = np.floor(jday) # truncate to integer.
+ # utc hour.
+ ut = d.hour + d.minute/60. + d.second/3600.
+ # calculate number of centuries from J2000
+ t = (jd + (ut/24.) - 2451545.0) / 36525.
+ # mean longitude corrected for aberration
+ l = (280.460 + 36000.770 * t) % 360
+ # mean anomaly
+ g = 357.528 + 35999.050 * t
+ # ecliptic longitude
+ lm = l + 1.915 * np.sin(g*dg2rad) + 0.020 * np.sin(2*g*dg2rad)
+ # obliquity of the ecliptic
+ ep = 23.4393 - 0.01300 * t
+ # equation of time
+ eqtime = -1.915*np.sin(g*dg2rad) - 0.020*np.sin(2*g*dg2rad) \
+ + 2.466*np.sin(2*lm*dg2rad) - 0.053*np.sin(4*lm*dg2rad)
+ # Greenwich hour angle
+ gha = 15*ut - 180 + eqtime
+ # declination of sun
+ dec = np.arcsin(np.sin(ep*dg2rad) * np.sin(lm*dg2rad)) * rad2dg
+ return gha, dec
+
+def daynightgrid(date, nlons):
+ """
+ date is datetime object (assumed UTC).
+ nlons is # of longitudes used to compute terminator."""
+ nlats = ((nlons-1)/2)+1
+ dg2rad = np.pi/180.
+ lons = np.linspace(-180,180,nlons)
+ # compute greenwich hour angle and solar declination
+ # from datetime object (assumed UTC).
+ tau, dec = epem(date)
+ longitude = lons + tau
+ lats = np.arctan(-np.cos(longitude*dg2rad)/np.tan(dec*dg2rad))/dg2rad
+ lons2 = np.linspace(-180,180,nlons)
+ lats2 = np.linspace(-90,90,nlats)
+ lons2, lats2 = np.meshgrid(lons2,lats2)
+ daynight = np.ones(lons2.shape, np.float)
+ for nlon in range(nlons):
+ daynight[:,nlon] = np.where(lats2[:,nlon]>lats[nlon],0,daynight[:,nlon])
+ return lons2,lats2,daynight
+
+# now, in UTC time.
+d = datetime.utcnow()
+
+# miller projection
+map = Basemap(projection='mill',lon_0=0)
+# plot coastlines, draw label meridians and parallels.
+map.drawcoastlines()
+map.drawparallels(np.arange(-90,90,30),labels=[1,0,0,0])
+map.drawmeridians(np.arange(-180,180,60),labels=[0,0,0,1])
+# create grid of day=0, night=1
+lons,lats,daynight = daynightgrid(d,1441)
+x,y = map(lons, lats)
+# contour this grid with 1 contour level, specifying colors.
+# (gray for night, axis background for day)
+map.contourf(x,y,daynight,1,colors=[plt.gca().get_axis_bgcolor(),'0.7'])
+plt.title('Day/Night Map for %s (UTC)' %d )
+plt.show()
Deleted: trunk/toolkits/basemap/examples/terminator.py
===================================================================
--- trunk/toolkits/basemap/examples/terminator.py 2009-08-24 18:00:34 UTC (rev 7554)
+++ trunk/toolkits/basemap/examples/terminator.py 2009-08-24 19:03:32 UTC (rev 7555)
@@ -1,75 +0,0 @@
-import numpy as np
-from mpl_toolkits.basemap import Basemap, netcdftime
-import matplotlib.pyplot as plt
-from datetime import datetime
-
-def epem(date):
- """
- input: date - datetime object (assumed UTC)
- ouput: gha - Greenwich hour angle, the angle between the Greenwich
- meridian and the meridian containing the subsolar point.
- dec - solar declination.
- """
- dg2rad = np.pi/180.
- rad2dg = 1./dg2rad
- # compute julian day from UTC datetime object.
- jday = netcdftime.JulianDayFromDate(date)
- jd = np.floor(jday) # truncate to integer.
- # utc hour.
- ut = d.hour + d.minute/60. + d.second/3600.
- # calculate number of centuries from J2000
- t = (jd + (ut/24.) - 2451545.0) / 36525.
- # mean longitude corrected for aberration
- l = (280.460 + 36000.770 * t) % 360
- # mean anomaly
- g = 357.528 + 35999.050 * t
- # ecliptic longitude
- lm = l + 1.915 * np.sin(g*dg2rad) + 0.020 * np.sin(2*g*dg2rad)
- # obliquity of the ecliptic
- ep = 23.4393 - 0.01300 * t
- # equation of time
- eqtime = -1.915*np.sin(g*dg2rad) - 0.020*np.sin(2*g*dg2rad) \
- + 2.466*np.sin(2*lm*dg2rad) - 0.053*np.sin(4*lm*dg2rad)
- # Greenwich hour angle
- gha = 15*ut - 180 + eqtime
- # declination of sun
- dec = np.arcsin(np.sin(ep*dg2rad) * np.sin(lm*dg2rad)) * rad2dg
- return gha, dec
-
-def daynightgrid(date, nlons):
- """
- date is datetime object (assumed UTC).
- nlons is # of longitudes used to compute terminator."""
- nlats = ((nlons-1)/2)+1
- dg2rad = np.pi/180.
- lons = np.linspace(-180,180,nlons)
- # compute greenwich hour angle and solar declination
- # from datetime object (assumed UTC).
- tau, dec = epem(date)
- longitude = lons + tau
- lats = np.arctan(-np.cos(longitude*dg2rad)/np.tan(dec*dg2rad))/dg2rad
- lons2 = np.linspace(-180,180,nlons)
- lats2 = np.linspace(-90,90,nlats)
- lons2, lats2 = np.meshgrid(lons2,lats2)
- daynight = np.ones(lons2.shape, np.float)
- for nlon in range(nlons):
- daynight[:,nlon] = np.where(lats2[:,nlon]>lats[nlon],0,daynight[:,nlon])
- return lons2,lats2,daynight
-
-# now, in UTC time.
-d = datetime.utcnow()
-
-# miller projection
-map = Basemap(projection='mill',lon_0=0)
-# plot coastlines, draw label meridians and parallels.
-map.drawcoastlines()
-map.drawparallels(np.arange(-90,90,30),labels=[1,0,0,0])
-map.drawmeridians(np.arange(-180,180,60),labels=[0,0,0,1])
-# create grid of day=0, night=1
-lons,lats,daynight = daynightgrid(d,1441)
-x,y = map(lons, lats)
-# contour this grid with 1 contour level, specifying colors.
-# (gray for night, axis background for day)
-map.contourf(x,y,daynight,1,colors=[plt.gca().get_axis_bgcolor(),'0.7'])
-plt.title('Day/Night Map for %s (UTC)' %d )
-plt.show()
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2010-07-07 14:02:11
|
Revision: 8529
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=8529&view=rev
Author: jswhit
Date: 2010-07-07 14:02:04 +0000 (Wed, 07 Jul 2010)
Log Message:
-----------
update URLs
Modified Paths:
--------------
trunk/toolkits/basemap/examples/fcstmaps.py
trunk/toolkits/basemap/examples/fcstmaps_axesgrid.py
Modified: trunk/toolkits/basemap/examples/fcstmaps.py
===================================================================
--- trunk/toolkits/basemap/examples/fcstmaps.py 2010-07-07 13:59:29 UTC (rev 8528)
+++ trunk/toolkits/basemap/examples/fcstmaps.py 2010-07-07 14:02:04 UTC (rev 8529)
@@ -15,16 +15,20 @@
YYYYMMDD = datetime.datetime.today().strftime('%Y%m%d')
# set OpenDAP server URL.
-URLbase="http://nomad1.ncep.noaa.gov:9090/dods/mrf/mrf"
-URL=URLbase+YYYYMMDD+'/mrf'+YYYYMMDD
-print URL+'\n'
try:
+ URLbase="http://nomad1.ncep.noaa.gov:9090/dods/mrf/mrf"
+ URL=URLbase+YYYYMMDD+'/mrf'+YYYYMMDD
data = NetCDFFile(URL)
except:
- msg = """
+ try:
+ URLbase="http://nomad2.ncep.noaa.gov:9090/dods/mrf/mrf"
+ URL=URLbase+YYYYMMDD+'/mrf'+YYYYMMDD
+ data = NetCDFFile(URL)
+ except:
+ msg = """
opendap server not providing the requested data.
Try another date by providing YYYYMMDD on command line."""
- raise IOError, msg
+ raise IOError, msg
# read lats,lons,times.
Modified: trunk/toolkits/basemap/examples/fcstmaps_axesgrid.py
===================================================================
--- trunk/toolkits/basemap/examples/fcstmaps_axesgrid.py 2010-07-07 13:59:29 UTC (rev 8528)
+++ trunk/toolkits/basemap/examples/fcstmaps_axesgrid.py 2010-07-07 14:02:04 UTC (rev 8529)
@@ -17,16 +17,20 @@
YYYYMMDD = datetime.datetime.today().strftime('%Y%m%d')
# set OpenDAP server URL.
-URLbase="http://nomad1.ncep.noaa.gov:9090/dods/mrf/mrf"
-URL=URLbase+YYYYMMDD+'/mrf'+YYYYMMDD
-print URL+'\n'
try:
+ URLbase="http://nomad1.ncep.noaa.gov:9090/dods/mrf/mrf"
+ URL=URLbase+YYYYMMDD+'/mrf'+YYYYMMDD
data = NetCDFFile(URL)
except:
- msg = """
+ try:
+ URLbase="http://nomad2.ncep.noaa.gov:9090/dods/mrf/mrf"
+ URL=URLbase+YYYYMMDD+'/mrf'+YYYYMMDD
+ data = NetCDFFile(URL)
+ except:
+ msg = """
opendap server not providing the requested data.
Try another date by providing YYYYMMDD on command line."""
- raise IOError, msg
+ raise IOError, msg
# read lats,lons,times.
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2011-01-17 17:48:10
|
Revision: 8924
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=8924&view=rev
Author: jswhit
Date: 2011-01-17 17:48:03 +0000 (Mon, 17 Jan 2011)
Log Message:
-----------
try to use netCDF4 module first before falling back on included NetCDFFile function
Modified Paths:
--------------
trunk/toolkits/basemap/examples/ccsm_popgrid.py
trunk/toolkits/basemap/examples/fcstmaps.py
trunk/toolkits/basemap/examples/fcstmaps_axesgrid.py
trunk/toolkits/basemap/examples/plothighsandlows.py
trunk/toolkits/basemap/examples/ploticos.py
trunk/toolkits/basemap/examples/plotprecip.py
Modified: trunk/toolkits/basemap/examples/ccsm_popgrid.py
===================================================================
--- trunk/toolkits/basemap/examples/ccsm_popgrid.py 2011-01-17 17:28:07 UTC (rev 8923)
+++ trunk/toolkits/basemap/examples/ccsm_popgrid.py 2011-01-17 17:48:03 UTC (rev 8924)
@@ -24,7 +24,11 @@
import numpy.ma as ma
import numpy as np
import matplotlib.pyplot as plt
-from mpl_toolkits.basemap import Basemap, NetCDFFile
+from mpl_toolkits.basemap import Basemap
+try:
+ from netCDF4 import Dataset as NetCDFFile
+except ImportError:
+ from mpl_toolkits.basemap import NetCDFFile
# read in data from netCDF file.
infile = 'ccsm_popgrid.nc'
Modified: trunk/toolkits/basemap/examples/fcstmaps.py
===================================================================
--- trunk/toolkits/basemap/examples/fcstmaps.py 2011-01-17 17:28:07 UTC (rev 8923)
+++ trunk/toolkits/basemap/examples/fcstmaps.py 2011-01-17 17:48:03 UTC (rev 8924)
@@ -5,7 +5,11 @@
import sys
import numpy.ma as ma
import datetime
-from mpl_toolkits.basemap import Basemap, NetCDFFile, addcyclic, num2date
+from mpl_toolkits.basemap import Basemap, addcyclic, num2date
+try:
+ from netCDF4 import Dataset as NetCDFFile
+except ImportError:
+ from mpl_toolkits.basemap import NetCDFFile
# today's date is default.
Modified: trunk/toolkits/basemap/examples/fcstmaps_axesgrid.py
===================================================================
--- trunk/toolkits/basemap/examples/fcstmaps_axesgrid.py 2011-01-17 17:28:07 UTC (rev 8923)
+++ trunk/toolkits/basemap/examples/fcstmaps_axesgrid.py 2011-01-17 17:48:03 UTC (rev 8924)
@@ -6,8 +6,12 @@
import sys
import numpy.ma as ma
import datetime
-from mpl_toolkits.basemap import Basemap, NetCDFFile, addcyclic, num2date
+from mpl_toolkits.basemap import Basemap, addcyclic, num2date
from mpl_toolkits.axes_grid1 import AxesGrid
+try:
+ from netCDF4 import Dataset as NetCDFFile
+except ImportError:
+ from mpl_toolkits.basemap import NetCDFFile
# today's date is default.
Modified: trunk/toolkits/basemap/examples/plothighsandlows.py
===================================================================
--- trunk/toolkits/basemap/examples/plothighsandlows.py 2011-01-17 17:28:07 UTC (rev 8923)
+++ trunk/toolkits/basemap/examples/plothighsandlows.py 2011-01-17 17:48:03 UTC (rev 8924)
@@ -5,8 +5,12 @@
import numpy as np
import matplotlib.pyplot as plt
import sys
-from mpl_toolkits.basemap import Basemap, NetCDFFile, addcyclic
+from mpl_toolkits.basemap import Basemap, addcyclic
from scipy.ndimage.filters import minimum_filter, maximum_filter
+try:
+ from netCDF4 import Dataset as NetCDFFile
+except ImportError:
+ from mpl_toolkits.basemap import NetCDFFile
def extrema(mat,mode='wrap',window=10):
"""find the indices of local extrema (min and max)
Modified: trunk/toolkits/basemap/examples/ploticos.py
===================================================================
--- trunk/toolkits/basemap/examples/ploticos.py 2011-01-17 17:28:07 UTC (rev 8923)
+++ trunk/toolkits/basemap/examples/ploticos.py 2011-01-17 17:48:03 UTC (rev 8924)
@@ -1,7 +1,11 @@
-from mpl_toolkits.basemap import Basemap, NetCDFFile
+from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
from numpy import ma
+try:
+ from netCDF4 import Dataset as NetCDFFile
+except ImportError:
+ from mpl_toolkits.basemap import NetCDFFile
# read in orography of icosahedral global grid.
f = NetCDFFile('C02562.orog.nc')
lons = (180./np.pi)*f.variables['grid_center_lon'][:]
Modified: trunk/toolkits/basemap/examples/plotprecip.py
===================================================================
--- trunk/toolkits/basemap/examples/plotprecip.py 2011-01-17 17:28:07 UTC (rev 8923)
+++ trunk/toolkits/basemap/examples/plotprecip.py 2011-01-17 17:48:03 UTC (rev 8924)
@@ -1,4 +1,8 @@
-from mpl_toolkits.basemap import Basemap, cm, NetCDFFile
+from mpl_toolkits.basemap import Basemap, cm
+try:
+ from netCDF4 import Dataset as NetCDFFile
+except ImportError:
+ from mpl_toolkits.basemap import NetCDFFile
import numpy as np
import matplotlib.pyplot as plt
import copy
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <js...@us...> - 2011-02-08 12:58:58
|
Revision: 8959
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=8959&view=rev
Author: jswhit
Date: 2011-02-08 12:58:49 +0000 (Tue, 08 Feb 2011)
Log Message:
-----------
celestial allsky map example from Tom Loredo.
Added Paths:
-----------
trunk/toolkits/basemap/examples/allskymap.py
trunk/toolkits/basemap/examples/allskymap_cr_example.py
Added: trunk/toolkits/basemap/examples/allskymap.py
===================================================================
--- trunk/toolkits/basemap/examples/allskymap.py (rev 0)
+++ trunk/toolkits/basemap/examples/allskymap.py 2011-02-08 12:58:49 UTC (rev 8959)
@@ -0,0 +1,390 @@
+"""
+AllSkyMap is a subclass of Basemap, specialized for handling common plotting
+tasks for celestial data.
+
+It is essentially equivalent to using Basemap with full-sphere projections
+(e.g., 'hammer' or 'moll') and the `celestial` keyword set to `True`, but
+it adds a few new methods:
+
+* label_meridians for, well, labeling meridians with their longitude values;
+
+* geodesic, a replacement for Basemap.drawgreatcircle, that can correctly
+ handle geodesics that cross the limb of the map, and providing the user
+ easy control over clipping (which affects thick lines at or near the limb);
+
+* tissot, which overrides Basemap.tissot, correctly handling geodesics that
+ cross the limb of the map.
+
+Created Jan 2011 by Tom Loredo, based on Jeff Whitaker's code in Basemap's
+__init__.py module.
+"""
+
+from numpy import *
+import matplotlib.pyplot as pl
+from matplotlib.pyplot import *
+from mpl_toolkits.basemap import Basemap, pyproj
+from mpl_toolkits.basemap.pyproj import Geod
+
+__all__ = ['AllSkyMap']
+
+def angle_symbol(angle, round_to=1.0):
+ """
+ Return a string representing an angle, rounded and with a degree symbol.
+
+ This is adapted from code in mpl's projections.geo module.
+ """
+ value = np.round(angle / round_to) * round_to
+ if pl.rcParams['text.usetex'] and not pl.rcParams['text.latex.unicode']:
+ return r"$%0.0f^\circ$" % value
+ else:
+ return u"%0.0f\u00b0" % value
+
+
+class AllSkyMap(Basemap):
+ """
+ AllSkyMap is a subclass of Basemap, specialized for handling common plotting
+ tasks for celestial data.
+
+ It is essentially equivalent to using Basemap with full-sphere projections
+ (e.g., 'hammer' or 'moll') and the `celestial` keyword set to `True`, but
+ it adds a few new methods:
+
+ * label_meridians for, well, labeling meridians with their longitude values;
+
+ * geodesic, a replacement for Basemap.drawgreatcircle, that can correctly
+ handle geodesics that cross the limb of the map, and providing the user
+ easy control over clipping (which affects thick lines at or near the
+ limb);
+
+ * tissot, which overrides Basemap.tissot, correctly handling geodesics that
+ cross the limb of the map.
+ """
+
+ # Longitudes corresponding to east and west edges, reflecting the
+ # convention that 180 deg is the eastern edge, according to basemap's
+ # underlying projections:
+ east_lon = 180.
+ west_lon = 180.+1.e-10
+
+ def __init__(self,
+ projection='hammer',
+ lat_0=0., lon_0=0.,
+ suppress_ticks=True,
+ boundinglat=None,
+ fix_aspect=True,
+ anchor='C',
+ ax=None):
+
+ if projection != 'hammer' and projection !='moll':
+ raise ValueError('Only hammer and moll projections supported!')
+
+ # Use Basemap's init, enforcing the values of many parameters that
+ # aren't used or whose Basemap defaults would not be altered for all-sky
+ # celestial maps.
+ Basemap.__init__(self, llcrnrlon=None, llcrnrlat=None,
+ urcrnrlon=None, urcrnrlat=None,
+ llcrnrx=None, llcrnry=None,
+ urcrnrx=None, urcrnry=None,
+ width=None, height=None,
+ projection=projection, resolution=None,
+ area_thresh=None, rsphere=1.,
+ lat_ts=None,
+ lat_1=None, lat_2=None,
+ lat_0=lat_0, lon_0=lon_0,
+ suppress_ticks=suppress_ticks,
+ satellite_height=1.,
+ boundinglat=None,
+ fix_aspect=True,
+ anchor=anchor,
+ celestial=True,
+ ax=ax)
+
+ # Keep a local ref to lon_0 for hemisphere checking.
+ self._lon_0 = self.projparams['lon_0']
+ self._limb = None
+
+ def drawmapboundary(self,color='k',linewidth=1.0,fill_color=None,\
+ zorder=None,ax=None):
+ """
+ draw boundary around map projection region, optionally
+ filling interior of region.
+
+ .. tabularcolumns:: |l|L|
+
+ ============== ====================================================
+ Keyword Description
+ ============== ====================================================
+ linewidth line width for boundary (default 1.)
+ color color of boundary line (default black)
+ fill_color fill the map region background with this
+ color (default is no fill or fill with axis
+ background color).
+ zorder sets the zorder for filling map background
+ (default 0).
+ ax axes instance to use
+ (default None, use default axes instance).
+ ============== ====================================================
+
+ returns matplotlib.collections.PatchCollection representing map boundary.
+ """
+ # Just call the base class version, but keep a copy of the limb
+ # polygon for clipping.
+ self._limb = Basemap.drawmapboundary(self, color=color,
+ linewidth=linewidth, fill_color=fill_color, zorder=zorder, ax=ax)
+ return self._limb
+
+ def label_meridians(self, lons, fontsize=10, valign='bottom', vnudge=0,
+ halign='center', hnudge=0):
+ """
+ Label meridians with their longitude values in degrees.
+
+ This labels meridians with negative longitude l with the value 360-l;
+ for maps in celestial orientation, this means meridians to the right
+ of the central meridian are labeled from 360 to 180 (left to right).
+
+ `vnudge` and `hnudge` specify amounts in degress to nudge the labels
+ from their default placements, vertically and horizontally. This
+ values obey the map orientation, so to nudge to the right, use a
+ negative `hnudge` value.
+ """
+ # Run through (lon, lat) pairs, with lat=0 in each pair.
+ lats = len(lons)*[0.]
+ for lon,lat in zip(lons, lats):
+ x, y = self(lon+hnudge, lat+vnudge)
+ if lon < 0:
+ lon_lbl = 360 + lon
+ else:
+ lon_lbl = lon
+ pl.text(x, y, angle_symbol(lon_lbl), fontsize=fontsize,
+ verticalalignment=valign,
+ horizontalalignment=halign)
+
+ def east_hem(self, lon):
+ """
+ Return True if lon is in the eastern hemisphere of the map wrt lon_0.
+ """
+ if (lon-self._lon_0) % 360. <= self.east_lon:
+ return True
+ else:
+ return False
+
+ def geodesic(self, lon1, lat1, lon2, lat2, del_s=.01, clip=True, **kwargs):
+ """
+ Plot a geodesic curve from (lon1, lat1) to (lon2, lat2), with
+ points separated by arc length del_s. Return a list of Line2D
+ instances for the curves comprising the geodesic. If the geodesic does
+ not cross the map limb, there will be only a single curve; if it
+ crosses the limb, there will be two curves.
+ """
+
+ # TODO: Perhaps return a single Line2D instance when there is only a
+ # single segment, and a list of segments only when there are two segs?
+
+ # TODO: Check the units of del_s.
+
+ # This is based on Basemap.drawgreatcircle (which draws an *arc* of a
+ # great circle), but addresses a limitation of that method, supporting
+ # geodesics that cross the map boundary by breaking them into two
+ # segments, one in the eastern hemisphere and the other in the western.
+ gc = pyproj.Geod(a=self.rmajor,b=self.rminor)
+ az12,az21,dist = gc.inv(lon1,lat1,lon2,lat2)
+ npoints = int((dist+0.5**del_s)/del_s)
+ # Calculate lon & lat for points on the arc.
+ lonlats = gc.npts(lon1,lat1,lon2,lat2,npoints)
+ lons = [lon1]; lats = [lat1]
+ for lon, lat in lonlats:
+ lons.append(lon)
+ lats.append(lat)
+ lons.append(lon2); lats.append(lat2)
+ # Break the arc into segments as needed, when there is a longitudinal
+ # hemisphere crossing.
+ segs = []
+ seg_lons, seg_lats = [lon1], [lat1]
+ cur_hem = self.east_hem(lon1)
+ for lon, lat in zip(lons[1:], lats[1:]):
+ if self.east_hem(lon) == cur_hem:
+ seg_lons.append(lon)
+ seg_lats.append(lat)
+ else:
+ # We should interpolate a new pt at the boundary, but in
+ # the mean time just rely on the step size being small.
+ segs.append( (seg_lons, seg_lats) )
+ seg_lons, seg_lats = [lon], [lat]
+ cur_hem = not cur_hem
+ segs.append( (seg_lons, seg_lats) )
+ # Plot each segment; return a list of the mpl lines.
+ lines = []
+ for lons, lats in segs:
+ x, y = self(lons, lats)
+ if clip and self._limb:
+ line = plot(x, y, clip_path=self._limb, **kwargs)[0]
+ else:
+ line = plot(x, y, **kwargs)[0]
+ lines.append(line)
+ # If there are multiple segments and no color args, reconcile the
+ # colors, which mpl will have autoset to different values.
+ # *** Does this screw up mpl's color set sequence for later lines?
+ if not kwargs.has_key('c') or kwargs.has_key('color'):
+ if len(lines) > 1:
+ c1 = lines[0].get_color()
+ for line in lines[1:]:
+ line.set_color(c1)
+ return lines
+
+ def tissot(self,lon_0,lat_0,radius_deg,npts,ax=None,**kwargs):
+ """
+ Draw a polygon centered at ``lon_0,lat_0``. The polygon
+ approximates a circle on the surface of the earth with radius
+ ``radius_deg`` degrees latitude along longitude ``lon_0``,
+ made up of ``npts`` vertices.
+
+ The polygon represents a Tissot's indicatrix
+ (http://en.wikipedia.org/wiki/Tissot's_Indicatrix),
+ which when drawn on a map shows the distortion inherent in the map
+ projection. Tissots can be used to display azimuthally symmetric
+ directional uncertainties ("error circles").
+
+ Extra keyword ``ax`` can be used to override the default axis instance.
+
+ Other \**kwargs passed on to matplotlib.patches.Polygon.
+
+ returns a list of matplotlib.patches.Polygon objects, with two polygons
+ when the tissot crosses the limb, and just one polygon otherwise.
+ """
+
+ # TODO: Just return the polygon (not a list) when there is only one
+ # polygon? Or stick with the list for consistency?
+
+ # This is based on Basemap.tissot, but addresses a limitation of that
+ # method by handling tissots that cross the limb of the map by finding
+ # separate polygons in the eastern and western hemispheres comprising
+ # the tissot.
+ ax = kwargs.pop('ax', None) or self._check_ax()
+ g = pyproj.Geod(a=self.rmajor,b=self.rminor)
+ az12,az21,dist = g.inv(lon_0,lat_0,lon_0,lat_0+radius_deg)
+ start_hem = self.east_hem(lon_0)
+ segs1 = [self(lon_0,lat_0+radius_deg)]
+ over, segs2 = [], []
+ delaz = 360./npts
+ az = az12
+ last_lon = lon_0
+ # Note adjacent and opposite edge longitudes, in case the tissot
+ # runs over the edge.
+ if start_hem: # eastern case
+ adj_lon = self.east_lon
+ opp_lon = self.west_lon
+ else:
+ adj_lon = self.west_lon
+ opp_lon = self.east_lon
+ for n in range(npts):
+ az = az+delaz
+ # skip segments along equator (Geod can't handle equatorial arcs)
+ if np.allclose(0.,lat_0) and (np.allclose(90.,az) or np.allclose(270.,az)):
+ continue
+ else:
+ lon, lat, az21 = g.fwd(lon_0, lat_0, az, dist)
+ # If in the starting hemisphere, add to 1st polygon seg list.
+ if self.east_hem(lon) == start_hem:
+ x, y = self(lon, lat)
+ # Add segment if it is in the map projection region.
+ if x < 1.e20 and y < 1.e20:
+ segs1.append( (x, y) )
+ last_lon = lon
+ # Otherwise, we cross hemispheres.
+ else:
+ # Trace the edge of each hemisphere.
+ x, y = self(adj_lon, lat)
+ if x < 1.e20 and y < 1.e20:
+ segs1.append( (x, y) )
+ # We presume if adj projection is okay, opposite is.
+ segs2.append( self(opp_lon, lat) )
+ # Also store the overlap in the opposite hemisphere.
+ x, y = self(lon, lat)
+ if x < 1.e20 and y < 1.e20:
+ over.append( (x, y) )
+ last_lon = lon
+ poly1 = Polygon(segs1, **kwargs)
+ ax.add_patch(poly1)
+ if segs2:
+ over.reverse()
+ segs2.extend(over)
+ poly2 = Polygon(segs2, **kwargs)
+ ax.add_patch(poly2)
+ return [poly1, poly2]
+ else:
+ return [poly1]
+
+
+if __name__ == '__main__':
+
+ # Note that Hammer & Mollweide projections enforce a 2:1 aspect ratio.
+ # Use figure size good for a 2:1 plot.
+ fig = figure(figsize=(12,6))
+
+ # Set up the projection and draw a grid.
+ map = AllSkyMap(projection='hammer')
+ # Save the bounding limb to use as a clip path later.
+ limb = map.drawmapboundary(fill_color='white')
+ map.drawparallels(np.arange(-75,76,15), linewidth=0.5, dashes=[1,2],
+ labels=[1,0,0,0], fontsize=9)
+ map.drawmeridians(np.arange(-150,151,30), linewidth=0.5, dashes=[1,2])
+
+ # Label a subset of meridians.
+ lons = np.arange(-150,151,30)
+ map.label_meridians(lons, fontsize=9, vnudge=1,
+ halign='left', hnudge=-1) # hnudge<0 shifts to right
+
+ # x, y limits are [0, 4*rt2], [0, 2*rt2].
+ rt2 = sqrt(2)
+
+ # Draw a slanted green line crossing the map limb.
+ line = plot([rt2,0], [rt2,2*rt2], 'g-')
+
+ # Draw a slanted magenta line crossing the map limb but clipped.
+ line = plot([rt2+.1,0+.1], [rt2,2*rt2], 'm-', clip_path=limb)
+
+ # Draw some geodesics.
+ # First a transparent thick blue geodesic crossing the limb but not clipped,
+ # overlayed by a thinner red geodesic that is clipped (by default), to
+ # illustrate the effect of clipping.
+ lines = map.geodesic(120, 30, 240, 60, clip=False, c='b', lw=7, alpha=.5)
+ lines = map.geodesic(240, 60, 120, 30, c='r', lw=3, alpha=.5)
+
+ # Next two large limb-crossing geodesics with the same path, but rendered
+ # in opposite directions, one transparent blue, the other transparent
+ # yellow. They should be right on top of each other, giving a greenish
+ # brown hue.
+ lines = map.geodesic(240, -60, 120, 30, c='b', lw=2, alpha=.5)
+ lines = map.geodesic(120, 30, 240, -60, c='y', lw=2, alpha=.5)
+
+ # What happens if a geodesic is given coordinates spanning more than
+ # a single rotation? Not sure what to expect, but it shoots off the
+ # map (clipped here). Perhaps we should ensure lons are in [0, 360].
+ #lines = map.geodesic(120, 20, 240+360, 50, del_s=.2, c='g')
+
+ # Two tissots fully within the limb.
+ poly = map.tissot(60, -15, 10, 100)
+ poly = map.tissot(280, 60, 10, 100)
+ #poly = map.tissot(90, -85, 10, 100)
+
+ # Limb-spanning tissots in each quadrant.
+ # lower left:
+ poly = map.tissot(170, -60, 15, 100)
+ # upper left:
+ poly = map.tissot(175, 70, 15, 100)
+ # upper right (note negative longitude):
+ poly = map.tissot(-175, 30, 15, 100, color='r', alpha=.6)
+ # lower right:
+ poly = map.tissot(185, -40, 10, 100)
+
+ # Plot the tissot centers as "+" symbols. Note the top left symbol
+ # would cross the limb without the clip_path argument; this might be
+ # desired to enhance visibility.
+ lons = [170, 175, -175, 185]
+ lats = [-60, 70, 30, -40]
+ x, y = map(lons, lats)
+ map.scatter(x, y, s=40, marker='+', linewidths=1, edgecolors='g',
+ facecolors='none', clip_path=limb, zorder=10) # hi zorder -> top
+
+ title('AllSkyMap demo: Clipped lines, markers, geodesics, tissots')
+ show()
Added: trunk/toolkits/basemap/examples/allskymap_cr_example.py
===================================================================
--- trunk/toolkits/basemap/examples/allskymap_cr_example.py (rev 0)
+++ trunk/toolkits/basemap/examples/allskymap_cr_example.py 2011-02-08 12:58:49 UTC (rev 8959)
@@ -0,0 +1,225 @@
+"""
+Example of astronomical use of AllSkyMap class in allskymap.py module
+
+Plot an all-sky map showing locations of the 27 highest-energy ultra-high
+energy cosmic rays detected by the Auger (south) experiment as of Aug 2007,
+and locations of 18 (fictitious!) candidate sources. Indicate CR direction
+uncertainties and source scattering scales with tissots, and show the
+nearest candidate source to each CR with geodesics.
+
+Created 2011-02-07 by Tom Loredo
+"""
+
+from cStringIO import StringIO
+import numpy as np
+from numpy import cos, sin, arccos, deg2rad, rad2deg
+import csv, re
+
+import matplotlib.pyplot as plt
+from allskymap import AllSkyMap
+from matplotlib.colors import Normalize
+from matplotlib.colorbar import ColorbarBase
+import matplotlib.ticker as ticker
+
+
+class Source:
+ """
+ Parse and store data for a celestial source.
+ """
+
+ int_re = re.compile(r'^[-+]?[0-9]+$')
+ # float_re = re.compile(r'^[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?$')
+
+ def __init__(self, id, year, day, l, b, sig=None, **kwds):
+ self.id = int(id)
+ self.year = int(year)
+ self.day = int(day)
+ self.l = float(l)
+ self._l = deg2rad(self.l) # radians
+ self.b = float(b)
+ self._b = deg2rad(self.b) # radians
+ if sig is not None:
+ self.sig = float(sig)
+ self._sig = deg2rad(self.sig)
+ else:
+ self.sig, self._sig = None, None
+ # If extra values are specified as keywords, set them as
+ # attributes. The quick way is to use self.__dict__.update(kwds),
+ # but we want to convert to int or float values if possible.
+ # *** Consider using matplotlib.cbook.converter.
+ if kwds:
+ for key, val in kwds.items():
+ try:
+ nval = float(val)
+ # It's a number; but it may be an int.
+ if self.int_re.match(str(val)):
+ nval = int(val)
+ setattr(self, key, nval)
+ except ValueError: # non-numerical value
+ setattr(self, key, val)
+
+ def gcangle(self, src):
+ """
+ Calculate the great circle angle to another source.
+ """
+ # Use the law of cosines; note it is usually expressed with colattitude.
+ c = sin(self._b)*sin(src._b) + \
+ cos(self._b)*cos(src._b)*cos(self._l - src._l)
+ return rad2deg(arccos(c))
+
+
+# Auger UHE cosmic ray data, Jan 2004 to Aug 2007
+# From Appendix A of Abraham et al. (2008); "Correlation of the highest-energy
+# cosmic rays with the positions of nearby active galactic nuclei,"
+# Astropart.Phys.29:188-204,2008; Erratum-ibid.30:45,2008
+
+# Year day ang S(1000) E (EeV) RA Dec Longitude Latitude
+# * = w/i 3.2 deg of AGN
+AugerData = StringIO(
+"""2004 125 47.7 252 70 267.1 -11.4 15.4 8.4
+2004 142 59.2 212 84 199.7 -34.9 -50.8 27.6 *
+2004 282 26.5 328 66 208.0 -60.3 -49.6 1.7 *
+2004 339 44.7 316 83 268.5 -61.0 -27.7 -17.0 *
+2004 343 23.4 323 63 224.5 -44.2 -34.4 13.0 *
+2005 54 35.0 373 84 17.4 -37.9 -75.6 -78.6 *
+2005 63 54.5 214 71 331.2 -1.2 58.8 -42.4 *
+2005 81 17.2 308 58 199.1 -48.6 -52.8 14.1 *
+2005 295 15.4 311 57 332.9 -38.2 4.2 -54.9 *
+2005 306 40.1 248 59 315.3 -0.3 48.8 -28.7 *
+2005 306 14.2 445 84 114.6 -43.1 -103.7 -10.3
+2006 35 30.8 398 85 53.6 -7.8 -165.9 -46.9 *
+2006 55 37.9 255 59 267.7 -60.7 -27.6 -16.5 *
+2006 81 34.0 357 79 201.1 -55.3 -52.3 7.3
+2006 185 59.1 211 83 350.0 9.6 88.8 -47.1 *
+2006 296 54.0 208 69 52.8 -4.5 -170.6 -45.7 *
+2006 299 26.0 344 69 200.9 -45.3 -51.2 17.2 *
+2007 13 14.3 762 148 192.7 -21.0 -57.2 41.8
+2007 51 39.2 247 58 331.7 2.9 63.5 -40.2 *
+2007 69 30.4 332 70 200.2 -43.4 -51.4 19.2 **
+2007 84 17.3 340 64 143.2 -18.3 -109.4 23.8 *
+2007 145 23.9 392 78 47.7 -12.8 -163.8 -54.4 *
+2007 186 44.8 248 64 219.3 -53.8 -41.7 5.9
+2007 193 18.0 469 90 325.5 -33.5 12.1 -49.0 *
+2007 221 35.3 318 71 212.7 -3.3 -21.8 54.1 *
+2007 234 33.2 365 80 185.4 -27.9 -65.1 34.5
+2007 235 42.6 276 69 105.9 -22.9 -125.2 -7.7
+""")
+AugerTable = csv.reader(AugerData, dialect='excel-tab')
+CRs = {}
+for id, row in enumerate(AugerTable):
+ # Make an integer ID from Year+Day (presumes none on same day!).
+ src = Source(id, row[0], row[1], row[7], row[8], E=float(row[4]))
+ CRs[src.id] = src
+print 'Parsed data for', len(CRs), 'UHE CRs...'
+
+# Partly fictitious candidate source locations.
+# src.id src.l_deg src.b_deg src.xProj src.yProj
+# tab-delimited
+CandData = StringIO(
+"""1 270. -28.
+2 229. -80.
+3 141. -47.
+4 172. -51.
+5 251. -51.
+6 241. -36.
+7 281. 26.
+8 241. 64.
+9 240. 64.
+10 148. 70.
+11 305. 13.
+12 98. 79.
+13 309. 19.
+14 104. 68.
+15 104. 68.
+16 321. 15.
+17 328. -14.
+18 177.5 -35.
+""")
+# Add this line above to see a tissot overlapping the map limb.
+CandTable = csv.reader(CandData, dialect='excel-tab')
+cands = {}
+for row in CandTable:
+ src = Source(row[0], 0, 0, row[1], row[2])
+ cands[src.id] = src
+print 'Parsed data for', len(cands), 'candidate sources...'
+
+# Calculate the separation matrix; track the closest candidate to each CR.
+sepn = {}
+for cr in CRs.values():
+ id, sep = None, 181.
+ for cand in cands.values():
+ val = cr.gcangle(cand)
+ sepn[cr.id, cand.id] = val
+ if val < sep:
+ id, sep = cand.id, val
+ # Store the closest candidate id and separation as a CR attribute.
+ cr.near_id = id
+ cr.near_ang = sep
+
+
+# Note that Hammer & Mollweide projections enforce a 2:1 aspect ratio.
+# Choose figure size for a 2:1 plot, with room at bottom for colorbar.
+fig = plt.figure(figsize=(12,7))
+main_ax = plt.axes([0.05, .19, .9, .75]) # rect=L,B,W,H
+
+# Set up the projection and draw a grid.
+m = AllSkyMap(ax=main_ax, projection='hammer')
+m.drawmapboundary(fill_color='white')
+m.drawparallels(np.arange(-75,76,15), linewidth=0.5, dashes=[1,2],
+ labels=[1,0,0,0], fontsize=9)
+m.drawmeridians(np.arange(-150,151,30), linewidth=0.5, dashes=[1,2])
+
+# Label a subset of meridians.
+lons = np.arange(-150,151,30)
+m.label_meridians(lons, fontsize=9, vnudge=1,
+ halign='left', hnudge=-1) # nudge<0 shifts to right
+
+# Plot CR directions.
+lvals = [src.l for src in CRs.values()]
+bvals = [src.b for src in CRs.values()]
+x, y = m(lvals, bvals)
+# These symbols will be covered by opaque tissots; plot them anyway
+# so there is a collection for the legend.
+cr_pts = m.scatter(x, y, s=8, c='r', marker='o', linewidths=.5,
+ edgecolors='none')
+
+# Plot tissots showing uncertainties, colored by energy.
+# We use 1 deg uncertainties, which are probably ~2 sigma for most events.
+Evals = np.array([src.E for src in CRs.values()])
+norm_E = Normalize(Evals.min()-10, Evals.max()+20) # -+ for jet_r for brt clrs
+# jet_r color map is in spectral sequence, with a small unsaturated
+# green range, helping to distinguish CRs from candidates.
+cmap = plt.cm.get_cmap('jet_r')
+for cr in CRs.values():
+ color = cmap(norm_E(cr.E))[0:3] # ignore alpha
+ m.tissot(cr.l, cr.b, 2., 30, ec='none', color=color, alpha=1)
+
+# Plot candidate directions.
+lvals = [src.l for src in cands.values()]
+bvals = [src.b for src in cands.values()]
+x, y = m(lvals, bvals)
+cand_pts = m.scatter(x, y, marker='+', linewidths=.5,
+ edgecolors='k', facecolors='none', zorder=10) # hi zorder -> top
+
+# Plot tissots showing possible scale of candidate scatter.
+for l, b in zip(lvals, bvals):
+ m.tissot(l, b, 5., 30, ec='none', color='g', alpha=0.25)
+
+# Show the closest candidate to each CR.
+for cr in CRs.values():
+ cand = cands[cr.near_id]
+ m.geodesic(cr.l, cr.b, cand.l, cand.b, lw=0.5, ls='-', c='g')
+
+plt.title('UHE Cosmic Rays and Candidate Sources')
+plt.legend([cr_pts, cand_pts], ['UHE CR', 'Candidate'],
+ frameon=False, loc='lower right', scatterpoints=1)
+
+# Plot a colorbar for the CR energies.
+cb_ax = plt.axes([0.25, .1, .5, .03], frameon=False) # rect=L,B,W,H
+#bar = ColorbarBase(cb_ax, cmap=cmap, orientation='horizontal', drawedges=False)
+vals = np.linspace(Evals.min(), Evals.max(), 100)
+bar = ColorbarBase(cb_ax, values=vals, norm=norm_E, cmap=cmap,
+ orientation='horizontal', drawedges=False)
+bar.set_label('CR Energy (EeV)')
+
+plt.show()
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
×
Want the latest updates on software, tech news, and AI?
Get latest updates about software, tech news, and AI from SourceForge directly in your inbox once a month.
Thanks for helping keep SourceForge clean.
X