Working with data in grids and making maps

Warning: This tutorial is still part of the ongoing implementation of the fatiando.gridder.Grid class (PR 128). It will not work with the current development version. You'll need the grid branch in order for this to work.

Projections rely on cartopy. To install cartopy using conda:

conda install cartopy -c rsignell

License

All content can be freely used and adapted under the terms of the Creative Commons Attribution 4.0 International License.



In [1]:

    
%matplotlib inline
from fatiando.gridder import Grid
from fatiando.vis import mpl
from fatiando import utils
import fatiando
import numpy as np



In [2]:

    
import cartopy.crs as crs
import cartopy.feature as cf



In [3]:

    
fatiando.__version__









    Out[3]:





'0.2-206-gfc5f461'

Regular grids



In [4]:

    
g = Grid.regular((-1.5e3, 1.5e3, -1e3, 1e3), (51, 51), z=0)



In [5]:

    
g.add_attribute('data', utils.gaussian2d(g.x, g.y, 1000, 500) - utils.gaussian2d(g.x, g.y, 500, 1000))









    Out[5]:





<fatiando.gridder.Grid at 0x7f8153d2ff90>



In [6]:

    
g.contourf('data')









    Out[6]:





<matplotlib.contour.QuadContourSet instance at 0x7f8151982830>



In [7]:

    
g[g.y > 0].interp((50, 50)).contourf('data')









    Out[7]:





<matplotlib.contour.QuadContourSet instance at 0x7f815182dab8>



In [8]:

    
g.cut([0, 1000, 0, 500]).contourf('data')









    Out[8]:





<matplotlib.contour.QuadContourSet instance at 0x7f81516a65a8>

.x is the North-South coordinate.



In [9]:

    
g.x









    Out[9]:





array([-1500., -1440., -1380., ...,  1380.,  1440.,  1500.])

.y is the East-West coordinate.



In [10]:

    
g.y









    Out[10]:





array([-1000., -1000., -1000., ...,  1000.,  1000.,  1000.])



In [11]:

    
g.y.reshape(g.shape)









    Out[11]:





array([[-1000., -1000., -1000., ..., -1000., -1000., -1000.],
       [ -960.,  -960.,  -960., ...,  -960.,  -960.,  -960.],
       [ -920.,  -920.,  -920., ...,  -920.,  -920.,  -920.],
       ..., 
       [  920.,   920.,   920., ...,   920.,   920.,   920.],
       [  960.,   960.,   960., ...,   960.,   960.,   960.],
       [ 1000.,  1000.,  1000., ...,  1000.,  1000.,  1000.]])



In [12]:

    
g.z









    Out[12]:





array([ 0.,  0.,  0., ...,  0.,  0.,  0.])



In [13]:

    
g.area









    Out[13]:





(-1500.0, 1500.0, -1000.0, 1000.0)



In [14]:

    
g.spacing









    Out[14]:





(60.0, 40.0)



In [15]:

    
g.add_attribute('data', utils.gaussian2d(g.x, g.y, 1000, 500) - utils.gaussian2d(g.x, g.y, 500, 1000))









    Out[15]:





<fatiando.gridder.Grid at 0x7f8153d2ff90>



In [16]:

    
mpl.figure(figsize=(5, 5))
mpl.subplot(221)
mpl.axis('scaled')
g.contourf('data')

mpl.subplot(222)
mpl.axis('scaled')
g.contour('data')

mpl.subplot(223)
mpl.axis('scaled')
g.contour('data', colors='k')

mpl.subplot(224)
mpl.axis('scaled')
g.pcolor('data')
mpl.tight_layout(pad=0.0)









    



/home/leo/bin/anaconda/lib/python2.7/site-packages/matplotlib/text.py:52: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal
  if rotation in ('horizontal', None):
/home/leo/bin/anaconda/lib/python2.7/site-packages/matplotlib/text.py:54: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal
  elif rotation == 'vertical':



In [17]:

    
g.add_attribute('positive_data', 10*utils.gaussian2d(g.x, g.y, 1000, 500)).contourf('positive_data')
mpl.axis('scaled')









    Out[17]:





(-1000.0, 1000.0, -1500.0, 1500.0)



In [18]:

    
g.add_attribute('negative_data', -10*utils.gaussian2d(g.x, g.y, 1000, 500)).contourf('negative_data')
mpl.axis('scaled')









    Out[18]:





(-1000.0, 1000.0, -1500.0, 1500.0)



In [19]:

    
p = g.profile([-1400, 0], [1400, 0], 100)



In [20]:

    
mpl.plot(p.distance, p.positive_data, '-r')
mpl.plot(p.distance, p.negative_data, '-b')
mpl.xlim(0, p.distance.max())









    Out[20]:





(0, 2800.0)

Scattered points



In [21]:

    
scat = Grid.scatter((-180, 180, -90, 90), 200, seed=0)
scat.add_attribute('data', 10 + utils.gaussian2d(scat.x, scat.y, 100, 20))









    Out[21]:





<fatiando.gridder.Grid at 0x7f81504a2cd0>



In [22]:

    
scat.contour('data')









    



---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-22-8747f569fb29> in <module>()
----> 1 scat.contour('data')

/home/leo/src/fatiando/fatiando/gridder.pyc in contour(self, attribute, levels, autoranges, ax, basemap, label, clabel, style, linewidth, **kwargs)
    571         assert self.shape is not None, \
    572             "'shape' not set. Can only plot regular grids. " + \
--> 573             "Try method 'interp' to produce a regular grid"
    574         if isinstance(attribute, str):
    575             attribute = getattr(self, attribute)

AssertionError: 'shape' not set. Can only plot regular grids. Try method 'interp' to produce a regular grid



In [24]:

    
g = scat.interp((50, 50))



In [25]:

    
mpl.figure(figsize=(8, 5))
mpl.subplot(221)
g.contourf(g.data)
scat.points()

mpl.subplot(222)
g.pcolor(g.data)
scat.points('xk')

mpl.subplot(223)
g.contour(g.data)

mpl.subplot(224)
g.contour(g.data, colors='k')
mpl.tight_layout(pad=0)









    



/home/leo/bin/anaconda/lib/python2.7/site-packages/matplotlib/contour.py:381: RuntimeWarning: invalid value encountered in true_divide
  dist = np.add.reduce(([(abs(s)[i] / L[i]) for i in range(xsize)]), -1)

With projections



In [26]:

    
grid = Grid.regular((-90, -35, -180, 180), (101, 101), z=250e3, lonlat=True, projection=crs.Geodetic())
grid.add_attribute('data', 500 + utils.gaussian2d(grid.lon, grid.lat, 40, 10, y0=-65))









    Out[26]:





<fatiando.gridder.Grid at 0x7f815046b250>



In [27]:

    
grid.area









    Out[27]:





(-90.0, -35.0, -180.0, 180.0)



In [28]:

    
grid.contourf('data', ax=mpl.axes(), transform=None)









    Out[28]:





<matplotlib.contour.QuadContourSet instance at 0x7f815054dab8>



In [29]:

    
proj = grid.reproject(crs.SouthPolarStereo())



In [30]:

    
fig, axes = mpl.subplots(1, 3, subplot_kw=dict(projection=proj.projection), figsize=(10, 3))
ax = axes[0]
proj.contourf('data', ax=ax)
ax.coastlines()
ax = axes[1]
proj.pcolor('data', ax=ax)
ax.coastlines()
ax = axes[2]
proj.contour('data', ax=ax, colors='k', levels=8)
ax.coastlines(color='b')
ax.gridlines()
mpl.tight_layout()



In [31]:

    
proj = grid.reproject(crs.Mollweide())
proj.contourf('data', colorbar='horizontal')
mpl.gca().coastlines()
mpl.gca().set_global()



In [32]:

    
proj = grid.reproject(crs.Miller())
proj.contourf('data', colorbar='horizontal')
mpl.gca().coastlines()









    Out[32]:





<cartopy.mpl.feature_artist.FeatureArtist at 0x7f81502d2dd0>

Input the data in the projected coordinates.



In [33]:

    
grid = Grid.scatter((-90, -15, -180, 180), 100, z=250e3, seed=0, projection=crs.Geodetic()).reproject(crs.SouthPolarStereo())
#  Add the attribute in the projected coordinates
grid.add_attribute('data', utils.gaussian2d(grid.x, grid.y, 5000e3, 2000e3) - utils.gaussian2d(grid.x, grid.y, 2000e3, 5000e3))









    Out[33]:





<fatiando.gridder.Grid at 0x7f8151006410>



In [34]:

    
grid.interp((50, 50)).pcolor('data')
mpl.gca().coastlines()
grid.points()









    Out[34]:





[<matplotlib.lines.Line2D at 0x7f814affc910>]



In [35]:

    
proj = grid.reproject(crs.Mollweide()).interp((50, 50))
proj.contourf('data', colorbar='horizontal')
grid.reproject(proj.projection).points('xk')
mpl.gca().coastlines()
mpl.gca().set_global()



In [36]:

    
proj = grid.reproject(crs.Geodetic()).cut((-90, -40, -160, 20)).reproject(crs.SouthPolarStereo()).interp((50, 50))
proj.contour('data')
mpl.gca().coastlines()









    Out[36]:





<cartopy.mpl.feature_artist.FeatureArtist at 0x7f814ad65450>



In [37]:

    
proj = grid.reproject(crs.Miller()).interp((50, 50))
proj.contour('data', colorbar='horizontal', cmap=mpl.cm.seismic)
grid.reproject(proj.projection).points('xk')
mpl.gca().coastlines()
mpl.gca().set_global()

Load data from files

Text and numpy binary files

Save some data to a text file.



In [38]:

    
scat = Grid.scatter((-180, 180, -90, 90), 200, seed=0)
scat.add_attribute('data', 10 + utils.gaussian2d(scat.x, scat.y, 100, 20))
scat.add_attribute('other_data', -10 - utils.gaussian2d(scat.x, scat.y, 100, 20))
np.savetxt('tmp.txt', np.transpose([scat.x, scat.y, scat.data, scat.other_data]))

Not giving column names will assume columns are: x y data1 data2 data3 ...



In [39]:

    
g = Grid.load('tmp.txt')



In [40]:

    
g.get_attributes()









    Out[40]:





['data1', 'data2']



In [41]:

    
g.interp((50, 50)).pcolor('data1')









    Out[41]:





<matplotlib.collections.QuadMesh at 0x7f814ac7a090>



In [42]:

    
g.interp((50, 50)).pcolor('data2')









    Out[42]:





<matplotlib.collections.QuadMesh at 0x7f814abf4190>

If you give column names, they will be used to name the attributes.



In [43]:

    
g = Grid.load('tmp.txt', column_names='x,y,positive,negative'.split(','))
g.get_attributes()









    Out[43]:





['positive', 'negative']



In [44]:

    
g.interp((50, 50)).contourf('positive')









    Out[44]:





<matplotlib.contour.QuadContourSet instance at 0x7f814a96a368>



In [45]:

    
g.interp((50, 50)).contourf('negative')









    Out[45]:





<matplotlib.contour.QuadContourSet instance at 0x7f814a7d1ab8>

Loading numpy binary files works the same way.



In [46]:

    
np.save('tmp', np.transpose([scat.x, scat.y, scat.data, scat.other_data]))



In [47]:

    
g = Grid.load('tmp.npy', column_names=['x', 'y', 'positive', 'negative']).interp((50, 50))



In [48]:

    
g.contour(g.positive)









    Out[48]:





<matplotlib.contour.QuadContourSet instance at 0x7f814a3b9ef0>



In [49]:

    
g.contour(g.negative)









    Out[49]:





<matplotlib.contour.QuadContourSet instance at 0x7f814a2de758>

CSV files



In [50]:

    
scat = Grid.scatter((-180, 180, -90, 90), 200, seed=0)
scat.add_attribute('data', 10 + utils.gaussian2d(scat.x, scat.y, 100, 20))
scat.add_attribute('other_data', -10 - utils.gaussian2d(scat.x, scat.y, 100, 20))









    Out[50]:





<fatiando.gridder.Grid at 0x7f814a3f1ad0>



In [51]:

    
with open('tmp.csv', 'w') as f:
    f.write('x; y; positive_data; negative_data;   \n')
    np.savetxt(f, np.transpose([scat.x, scat.y, scat.data, scat.other_data]), delimiter='; ')



In [52]:

    
!head tmp.csv









    



x; y; positive_data; negative_data;   
1.757286141383690392e+01; -3.387674124106153783e+01; 1.005502517129078122e+01; -1.005502517129078122e+01
7.746817189407101978e+01; 3.534182798678270387e+01; 1.002416756632407946e+01; -1.002416756632407946e+01
3.699481538579178164e+01; -2.200466892735343549e+01; 1.025992158930449705e+01; -1.025992158930449705e+01
1.615794587888288447e+01; -5.767133803926573421e+01; 1.000023851794733076e+01; -1.000023851794733076e+01
-2.748427223799430408e+01; -8.555782888956038335e+01; 1.000000001045805043e+01; -1.000000001045805043e+01
5.252188070399620301e+01; -7.789506633661525825e+01; 1.000000019603782775e+01; -1.000000019603782775e+01
-2.246860394543068651e+01; 3.229069922974211693e+01; 1.007014386314265231e+01; -1.007014386314265231e+01
1.410382802815487366e+02; -8.334567979911838620e+00; 1.011499780830624040e+01; -1.011499780830624040e+01
1.669185937803705428e+02; 6.584257999569999242e+00; 1.005532320553900050e+01; -1.005532320553900050e+01

Will try to get the column names from the first line of the CSV file.



In [53]:

    
g = Grid.load('tmp.csv')
g.get_attributes()









    Out[53]:





['positive_data', 'negative_data']



In [54]:

    
g.interp((50, 50)).contourf('positive_data')









    Out[54]:





<matplotlib.contour.QuadContourSet instance at 0x7f8149947878>

You can specify the column names using the column_names argument. This will overwrite the values read from the file.



In [55]:

    
g = Grid.load('tmp.csv', column_names='lat lon height gravity'.split())
g.get_attributes()









    Out[55]:





['gravity', 'height']



In [56]:

    
g.interp((50, 50)).pcolor('gravity')









    Out[56]:





<matplotlib.collections.QuadMesh at 0x7f81497bb390>

If the column names are not present in the file, will assume that they are: x y data1 data2 data3 ...



In [57]:

    
np.savetxt('tmp.csv', np.transpose([scat.x, scat.y, scat.data, scat.other_data]), delimiter='; ')



In [58]:

    
!head tmp.csv









    



1.757286141383690392e+01; -3.387674124106153783e+01; 1.005502517129078122e+01; -1.005502517129078122e+01
7.746817189407101978e+01; 3.534182798678270387e+01; 1.002416756632407946e+01; -1.002416756632407946e+01
3.699481538579178164e+01; -2.200466892735343549e+01; 1.025992158930449705e+01; -1.025992158930449705e+01
1.615794587888288447e+01; -5.767133803926573421e+01; 1.000023851794733076e+01; -1.000023851794733076e+01
-2.748427223799430408e+01; -8.555782888956038335e+01; 1.000000001045805043e+01; -1.000000001045805043e+01
5.252188070399620301e+01; -7.789506633661525825e+01; 1.000000019603782775e+01; -1.000000019603782775e+01
-2.246860394543068651e+01; 3.229069922974211693e+01; 1.007014386314265231e+01; -1.007014386314265231e+01
1.410382802815487366e+02; -8.334567979911838620e+00; 1.011499780830624040e+01; -1.011499780830624040e+01
1.669185937803705428e+02; 6.584257999569999242e+00; 1.005532320553900050e+01; -1.005532320553900050e+01
-4.196105322272003946e+01; 7.140083274726157470e+01; 1.000000244547717365e+01; -1.000000244547717365e+01



In [59]:

    
g = Grid.load('tmp.csv')
g.get_attributes()









    Out[59]:





['data1', 'data2']



In [60]:

    
g.interp((50, 50)).contourf('data2')









    Out[60]:





<matplotlib.contour.QuadContourSet instance at 0x7f81495ac5f0>

.gdf files from ICGEM

.gdf files downloaded from http://icgem.gfz-potsdam.de/ICGEM/potato/Service.html.



In [61]:

    
g = Grid.load('go_cons_gcf_2_tim_r5-antartica.gdf')
print(g.metadata)









    



generating_institute     gfz-potsdam
     generating_date     2014/10/01
        product_type     gravity_field
                body     earth
           modelname     go_cons_gcf_2_tim_r5
     max_used_degree           280
         tide_system     tide_free
          functional     gravity_ell  (centrifugal term included)
                unit     mgal
          refsysname     WGS84
            gmrefpot      3.98600441800E+14 m**3/s**2
        radiusrefpot     6378137.000 m
          flatrefpot      3.352810664747480E-03   (1/298.25722356300)
         omegarefpot      7.29211500000E-05 1/s
       long_lat_unit     degree
      latlimit_north      -60.000000000000    
      latlimit_south      -90.000000000000    
      longlimit_west       0.0000000000000    
      longlimit_east       360.00000000000    
            gridstep      0.50000000000000    
     height_over_ell      10000.0000 m
  latitude_parallels            61
 longitude_parallels           721
number_of_gridpoints         43981
            gapvalue        9999999.0000
       weighted_mean      9.7945572E+05 mgal
            maxvalue      9.8016875E+05 mgal
            minvalue      9.7874044E+05 mgal
         signal_wrms      3.7810914E+02 mgal
         grid_format     long_lat_value
 
          longitude    latitude    gravity_ell
            [deg.]      [deg.]         [mgal]



In [62]:

    
g.get_attributes()









    Out[62]:





['gravity_ell', 'height']



In [63]:

    
proj = g.set_projection(crs.Geodetic()).reproject(crs.SouthPolarStereo())



In [64]:

    
proj.pcolor(proj.gravity_ell)
mpl.gca().coastlines()
mpl.title('Gravity (mGal) at {:g} km height'.format(0.001*g.height[0]))









    Out[64]:





<matplotlib.text.Text at 0x7f8149869a90>



In [65]:

    
g = Grid.load('eigen-6c3stat-antartica.gdf')
print(g.metadata)









    



generating_institute     gfz-potsdam
     generating_date     2014/10/01
        product_type     gravity_field
                body     earth
           modelname     eigen-6c3stat
     max_used_degree          1949
         tide_system     tide_free
          functional     gravity_earth  (centrifugal term included)
                unit     mgal
          refsysname     WGS84
            gmrefpot      3.98600441800E+14 m**3/s**2
        radiusrefpot     6378137.000 m
          flatrefpot      3.352810664747480E-03   (1/298.25722356300)
         omegarefpot      7.29211500000E-05 1/s
    normal_potential      6.263685171456948E+07 m**2/s**2
       long_lat_unit     degree
      latlimit_north      -60.000000000000    
      latlimit_south      -90.000000000000    
      longlimit_west       0.0000000000000    
      longlimit_east       360.00000000000    
            gridstep      0.50000000000000    
  latitude_parallels            61
 longitude_parallels           721
number_of_gridpoints         43981
            gapvalue        9999999.0000
       weighted_mean      9.8229713E+05 mgal
            maxvalue      9.8314427E+05 mgal
            minvalue      9.8171161E+05 mgal
         signal_wrms      3.1568614E+02 mgal
         grid_format     long_lat_height_value
 
          longitude    latitude  h_over_geoid  gravity_earth
            [deg.]      [deg.]     [meter]         [mgal]



In [66]:

    
g.get_attributes()









    Out[66]:





['gravity_earth', 'h_over_geoid']



In [67]:

    
proj = g.set_projection(crs.Geodetic()).reproject(crs.SouthPolarStereo())



In [68]:

    
fig, axes = mpl.subplots(1, 2, figsize=(11, 4), subplot_kw=dict(projection=proj.projection))
ax = axes[0]
proj.pcolor(proj.gravity_earth, ax=ax)
ax.coastlines()
ax.set_title('Gravity (mGal) at topography')
ax = axes[1]
proj.pcolor(proj.h_over_geoid, ax=ax, cmap=mpl.cm.Greens)
ax.coastlines()
ax.set_title('Topography (m)')
ax.gridlines()
mpl.tight_layout(pad=0)

Saving to CSV format



In [69]:

    
scat = Grid.scatter((-180, 180, -90, 90), 200, seed=0)
scat.add_attribute('my_data', 10 + utils.gaussian2d(scat.x, scat.y, 100, 20))
scat.dump_csv('tmp.csv')



In [70]:

    
!head tmp.csv









    



x; y; my_data
17.57286;-33.87674;10.05503
77.46817;35.34183;10.02417
36.99482;-22.00467;10.25992
16.15795;-57.67134;10.00024
-27.48427;-85.55783;10.00000
52.52188;-77.89507;10.00000
-22.46860;32.29070;10.07014
141.03828;-8.33457;10.11500
166.91859;6.58426;10.05532



In [71]:

    
Grid.load('tmp.csv').interp((50, 50)).contourf('my_data')









    Out[71]:





<matplotlib.contour.QuadContourSet instance at 0x7f8151239320>