/[escript]/branches/doubleplusgood/doc/examples/inversion/create_netcdf.py
ViewVC logotype

Diff of /branches/doubleplusgood/doc/examples/inversion/create_netcdf.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4190 by caltinay, Thu Jan 24 03:53:19 2013 UTC revision 4191 by gross, Mon Feb 11 08:22:47 2013 UTC
# Line 12  Line 12 
12  #  #
13  ##############################################################################  ##############################################################################
14    
15  """This example shows how to create compatible netCDF input files for inversion"""  """
16    This example shows how to create a netCDF input files for inversion with
17    inversion in esys.downunder.
18    """
19    
20  from datetime import datetime  from datetime import datetime
21  import numpy as np  import numpy as np
# Line 40  NY=10 Line 43  NY=10
43  # Dummy value (for unset areas)  # Dummy value (for unset areas)
44  MISSING=-9999  MISSING=-9999
45    
46  # Units of the data array  # Data error
47  UNITS='mGal'  SIGMA = 3.
48    
49    
50    
51  # The actual data array, must have shape (NY, NX)  # The actual data array, must have shape (NY, NX)
52  DATA=10*np.random.normal(size=(NY, NX))  # these are just some random numbers.
53    DATA=10*np.random.normal(size=(NY, NX), scale=SIGMA)
54    ERROR=np.ones((NY, NX)) * SIGMA
55    
56  ##############################################################################  ##############################################################################
57  ###################### Keep everything below this line #######################  ###################### Keep everything below this line #######################
58  ##############################################################################  ##############################################################################
   
 # conventions used in the file  
 conventions="CF-1.0, COARDS, Unidata Dataset Discovery v1.0"  
 # projection string  
 esri_pe_string="GEOGCS[\\\"GCS_WGS_1984\\\",DATUM[\\\"D_WGS_1984\\\",SPHEROID[\\\"WGS_1984\\\",6378137.0,298.257223563]],PRIMEM[\\\"Greenwich\\\",0.0],UNIT[\\\"Degree\\\",0.0174532925199433]]"  
59  # file log  # file log
60  history=datetime.now().strftime("%d-%m-%Y")+" created using python script"  history=datetime.now().strftime("%d-%m-%Y")+" created using python script"
61  # license  # license
# Line 63  longitude=np.linspace(ORIGIN_X, ORIGIN_X Line 65  longitude=np.linspace(ORIGIN_X, ORIGIN_X
65  latitude=np.linspace(ORIGIN_Y, ORIGIN_Y-NY*DELTA_Y, NY, endpoint=False)  latitude=np.linspace(ORIGIN_Y, ORIGIN_Y-NY*DELTA_Y, NY, endpoint=False)
66    
67  o=netcdf_file(FILENAME,'w')  o=netcdf_file(FILENAME,'w')
 o.cdm_data_type="Grid"  
 o.Conventions=conventions  
68  o.history=history  o.history=history
69  o.license=license  o.license=license
 o.Metadata_Conventions=conventions  
 o.summary=SUMMARY  
70  o.title=TITLE  o.title=TITLE
71    o.summary=SUMMARY
72  o.createDimension("latitude", NY)  o.createDimension("latitude", NY)
73  o.createDimension("longitude", NX)  o.createDimension("longitude", NX)
74    
75  v=o.createVariable("latitude", latitude.dtype, ["latitude"])  v=o.createVariable("latitude", latitude.dtype, ["latitude"])
 v.axis = "Y"  
 v.long_name="Latitude"  
 v.standard_name="latitude"  
 v.units="degrees_north"  
76  v.data[:]=latitude  v.data[:]=latitude
77    
78  v=o.createVariable("longitude", longitude.dtype, ["longitude"])  v=o.createVariable("longitude", longitude.dtype, ["longitude"])
 v.axis = "X"  
 v.long_name="Longitude"  
 v.standard_name="longitude"  
 v.units="degrees_east"  
79  v.data[:]=longitude  v.data[:]=longitude
80    
81  v=o.createVariable("data", DATA.dtype, ["latitude","longitude"])  v=o.createVariable("data", DATA.dtype, ["latitude","longitude"])
 v.coordinates="lon lat"  
 v.esri_pe_string=esri_pe_string  
 v.long_name="Bouguer_anomaly"  
82  v.missing_value=MISSING  v.missing_value=MISSING
 v.units=UNITS  
83  v.data[:]=DATA  v.data[:]=DATA
84    
85    # can be omitted.
86    v=o.createVariable("error", DATA.dtype, ["latitude","longitude"])
87    v.missing_value=MISSING
88    v.data[:]=ERROR
89    
90  o.close()  o.close()
91    

Legend:
Removed from v.4190  
changed lines
  Added in v.4191

  ViewVC Help
Powered by ViewVC 1.1.26