Changeset 7742


Ignore:
Timestamp:
May 24, 2010, 4:34:13 PM (15 years ago)
Author:
hudson
Message:

Further split up file conversion to modularise shallow_water module.

Location:
anuga_core/source/anuga
Files:
7 added
3 edited
1 copied

Legend:

Unmodified
Added
Removed

  • anuga_core/source/anuga/shallow_water/data_manager.py

    r7736 r7742  
    16561656    fid.close()
    16571657
    1658 
    1659 ##
    1660 # @brief Convert ASC file to DEM file.
    1661 # @param basename_in Stem of input filename.
    1662 # @param basename_out Stem of output filename.
    1663 # @param use_cache ??
    1664 # @param verbose True if this function is to be verbose.
    1665 # @return
    1666 # @note A PRJ file with same stem basename must exist and is used to fix the
    1667 #       UTM zone, datum, false northings and eastings.
    1668 def convert_dem_from_ascii2netcdf(basename_in, basename_out=None,
    1669                                   use_cache=False,
    1670                                   verbose=False):
    1671     """Read Digital Elevation model from the following ASCII format (.asc)
    1672 
    1673     Example:
    1674     ncols         3121
    1675     nrows         1800
    1676     xllcorner     722000
    1677     yllcorner     5893000
    1678     cellsize      25
    1679     NODATA_value  -9999
    1680     138.3698 137.4194 136.5062 135.5558 ..........
    1681 
    1682     Convert basename_in + '.asc' to NetCDF format (.dem)
    1683     mimicking the ASCII format closely.
    1684 
    1685     An accompanying file with same basename_in but extension .prj must exist
    1686     and is used to fix the UTM zone, datum, false northings and eastings.
    1687 
    1688     The prj format is assumed to be as
    1689 
    1690     Projection    UTM
    1691     Zone          56
    1692     Datum         WGS84
    1693     Zunits        NO
    1694     Units         METERS
    1695     Spheroid      WGS84
    1696     Xshift        0.0000000000
    1697     Yshift        10000000.0000000000
    1698     Parameters
    1699     """
    1700 
    1701     kwargs = {'basename_out': basename_out, 'verbose': verbose}
    1702 
    1703     if use_cache is True:
    1704         from caching import cache
    1705         result = cache(_convert_dem_from_ascii2netcdf, basename_in, kwargs,
    1706                        dependencies=[basename_in + '.asc',
    1707                                      basename_in + '.prj'],
    1708                        verbose=verbose)
    1709 
    1710     else:
    1711         result = apply(_convert_dem_from_ascii2netcdf, [basename_in], kwargs)
    1712 
    1713     return result
    1714 
    1715 
    1716 ##
    1717 # @brief Convert an ASC file to a DEM file.
    1718 # @param basename_in Stem of input filename.
    1719 # @param basename_out Stem of output filename.
    1720 # @param verbose True if this function is to be verbose.
    1721 def _convert_dem_from_ascii2netcdf(basename_in, basename_out = None,
    1722                                    verbose = False):
    1723     """Read Digital Elevation model from the following ASCII format (.asc)
    1724 
    1725     Internal function. See public function convert_dem_from_ascii2netcdf
    1726     for details.
    1727     """
    1728 
    1729     import os
    1730     from Scientific.IO.NetCDF import NetCDFFile
    1731 
    1732     root = basename_in
    1733 
    1734     # Read Meta data
    1735     if verbose: log.critical('Reading METADATA from %s' % (root + '.prj'))
    1736 
    1737     metadatafile = open(root + '.prj')
    1738     metalines = metadatafile.readlines()
    1739     metadatafile.close()
    1740 
    1741     L = metalines[0].strip().split()
    1742     assert L[0].strip().lower() == 'projection'
    1743     projection = L[1].strip()                   #TEXT
    1744 
    1745     L = metalines[1].strip().split()
    1746     assert L[0].strip().lower() == 'zone'
    1747     zone = int(L[1].strip())
    1748 
    1749     L = metalines[2].strip().split()
    1750     assert L[0].strip().lower() == 'datum'
    1751     datum = L[1].strip()                        #TEXT
    1752 
    1753     L = metalines[3].strip().split()
    1754     assert L[0].strip().lower() == 'zunits'     #IGNORE
    1755     zunits = L[1].strip()                       #TEXT
    1756 
    1757     L = metalines[4].strip().split()
    1758     assert L[0].strip().lower() == 'units'
    1759     units = L[1].strip()                        #TEXT
    1760 
    1761     L = metalines[5].strip().split()
    1762     assert L[0].strip().lower() == 'spheroid'   #IGNORE
    1763     spheroid = L[1].strip()                     #TEXT
    1764 
    1765     L = metalines[6].strip().split()
    1766     assert L[0].strip().lower() == 'xshift'
    1767     false_easting = float(L[1].strip())
    1768 
    1769     L = metalines[7].strip().split()
    1770     assert L[0].strip().lower() == 'yshift'
    1771     false_northing = float(L[1].strip())
    1772 
    1773     #Read DEM data
    1774     datafile = open(basename_in + '.asc')
    1775 
    1776     if verbose: log.critical('Reading DEM from %s' % (basename_in + '.asc'))
    1777 
    1778     lines = datafile.readlines()
    1779     datafile.close()
    1780 
    1781     if verbose: log.critical('Got %d lines' % len(lines))
    1782 
    1783     ncols = int(lines[0].split()[1].strip())
    1784     nrows = int(lines[1].split()[1].strip())
    1785 
    1786     # Do cellsize (line 4) before line 2 and 3
    1787     cellsize = float(lines[4].split()[1].strip())
    1788 
    1789     # Checks suggested by Joaquim Luis
    1790     # Our internal representation of xllcorner
    1791     # and yllcorner is non-standard.
    1792     xref = lines[2].split()
    1793     if xref[0].strip() == 'xllcorner':
    1794         xllcorner = float(xref[1].strip()) # + 0.5*cellsize # Correct offset
    1795     elif xref[0].strip() == 'xllcenter':
    1796         xllcorner = float(xref[1].strip())
    1797     else:
    1798         msg = 'Unknown keyword: %s' % xref[0].strip()
    1799         raise Exception, msg
    1800 
    1801     yref = lines[3].split()
    1802     if yref[0].strip() == 'yllcorner':
    1803         yllcorner = float(yref[1].strip()) # + 0.5*cellsize # Correct offset
    1804     elif yref[0].strip() == 'yllcenter':
    1805         yllcorner = float(yref[1].strip())
    1806     else:
    1807         msg = 'Unknown keyword: %s' % yref[0].strip()
    1808         raise Exception, msg
    1809 
    1810     NODATA_value = int(lines[5].split()[1].strip())
    1811 
    1812     assert len(lines) == nrows + 6
    1813 
    1814     if basename_out == None:
    1815         netcdfname = root + '.dem'
    1816     else:
    1817         netcdfname = basename_out + '.dem'
    1818 
    1819     if verbose: log.critical('Store to NetCDF file %s' % netcdfname)
    1820 
    1821     # NetCDF file definition
    1822     fid = NetCDFFile(netcdfname, netcdf_mode_w)
    1823 
    1824     #Create new file
    1825     fid.institution = 'Geoscience Australia'
    1826     fid.description = 'NetCDF DEM format for compact and portable storage ' \
    1827                       'of spatial point data'
    1828 
    1829     fid.ncols = ncols
    1830     fid.nrows = nrows
    1831     fid.xllcorner = xllcorner
    1832     fid.yllcorner = yllcorner
    1833     fid.cellsize = cellsize
    1834     fid.NODATA_value = NODATA_value
    1835 
    1836     fid.zone = zone
    1837     fid.false_easting = false_easting
    1838     fid.false_northing = false_northing
    1839     fid.projection = projection
    1840     fid.datum = datum
    1841     fid.units = units
    1842 
    1843     # dimension definitions
    1844     fid.createDimension('number_of_rows', nrows)
    1845     fid.createDimension('number_of_columns', ncols)
    1846 
    1847     # variable definitions
    1848     fid.createVariable('elevation', netcdf_float, ('number_of_rows',
    1849                                                    'number_of_columns'))
    1850 
    1851     # Get handles to the variables
    1852     elevation = fid.variables['elevation']
    1853 
    1854     #Store data
    1855     n = len(lines[6:])
    1856     for i, line in enumerate(lines[6:]):
    1857         fields = line.split()
    1858         if verbose and i % ((n+10)/10) == 0:
    1859             log.critical('Processing row %d of %d' % (i, nrows))
    1860            
    1861         if len(fields) != ncols:
    1862             msg = 'Wrong number of columns in file "%s" line %d\n' % (basename_in + '.asc', i)
    1863             msg += 'I got %d elements, but there should have been %d\n' % (len(fields), ncols)
    1864             raise Exception, msg
    1865 
    1866         elevation[i, :] = num.array([float(x) for x in fields])
    1867 
    1868     fid.close()
    18691658
    18701659

  • anuga_core/source/anuga/shallow_water/sww_file.py

    r7736 r7742  
    77from anuga.config import max_float
    88from anuga.utilities.numerical_tools import ensure_numeric
     9import anuga.utilities.log as log
    910from Scientific.IO.NetCDF import NetCDFFile
    1011

  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r7736 r7742  
    1818from anuga.shallow_water.data_manager import *
    1919from anuga.shallow_water.sww_file import SWW_file
    20 from anuga.shallow_water.file_conversion import tsh2sww, \
    21                         asc_csiro2sww, pmesh_to_domain_instance, \
    22                         dem2pts, _get_min_max_indexes
     20from anuga.file_conversion.file_conversion import tsh2sww, \
     21                        asc_csiro2sww, pmesh_to_domain_instance
     22                       
    2323from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
    2424from anuga.abstract_2d_finite_volumes.util import file_function
    … …  
    735735        #Cleanup
    736736        os.remove(sww.filename)
    737 
    738 
    739 
    740     def test_dem2pts_bounding_box_v2(self):
    741         """Test conversion from dem in ascii format to native NetCDF format
    742         """
    743 
    744         import time, os
    745         from Scientific.IO.NetCDF import NetCDFFile
    746 
    747         #Write test asc file
    748         root = 'demtest'
    749 
    750         filename = root+'.asc'
    751         fid = open(filename, 'w')
    752         fid.write("""ncols         10
    753 nrows         10
    754 xllcorner     2000
    755 yllcorner     3000
    756 cellsize      1
    757 NODATA_value  -9999
    758 """)
    759         #Create linear function
    760         ref_points = []
    761         ref_elevation = []
    762         x0 = 2000
    763         y = 3010
    764         yvec = range(10)
    765         xvec = range(10)
    766         z = -1
    767         for i in range(10):
    768             y = y - 1
    769             for j in range(10):
    770                 x = x0 + xvec[j]
    771                 z += 1
    772                 ref_points.append ([x,y])
    773                 ref_elevation.append(z)
    774                 fid.write('%f ' %z)
    775             fid.write('\n')
    776 
    777         fid.close()
    778 
    779         #print 'sending pts', ref_points
    780         #print 'sending elev', ref_elevation
    781 
    782         #Write prj file with metadata
    783         metafilename = root+'.prj'
    784         fid = open(metafilename, 'w')
    785 
    786 
    787         fid.write("""Projection UTM
    788 Zone 56
    789 Datum WGS84
    790 Zunits NO
    791 Units METERS
    792 Spheroid WGS84
    793 Xshift 0.0000000000
    794 Yshift 10000000.0000000000
    795 Parameters
    796 """)
    797         fid.close()
    798 
    799         #Convert to NetCDF pts
    800         convert_dem_from_ascii2netcdf(root)
    801         dem2pts(root, easting_min=2002.0, easting_max=2007.0,
    802                 northing_min=3003.0, northing_max=3006.0,
    803                 verbose=self.verbose)
    804 
    805         #Check contents
    806         #Get NetCDF
    807         fid = NetCDFFile(root+'.pts', netcdf_mode_r)
    808 
    809         # Get the variables
    810         #print fid.variables.keys()
    811         points = fid.variables['points']
    812         elevation = fid.variables['elevation']
    813 
    814         #Check values
    815         assert fid.xllcorner[0] == 2002.0
    816         assert fid.yllcorner[0] == 3003.0
    817 
    818         #create new reference points
    819         newz = []
    820         newz[0:5] = ref_elevation[32:38]
    821         newz[6:11] = ref_elevation[42:48]
    822         newz[12:17] = ref_elevation[52:58]
    823         newz[18:23] = ref_elevation[62:68]
    824         ref_elevation = []
    825         ref_elevation = newz
    826         ref_points = []
    827         x0 = 2002
    828         y = 3007
    829         yvec = range(4)
    830         xvec = range(6)
    831         for i in range(4):
    832             y = y - 1
    833             ynew = y - 3003.0
    834             for j in range(6):
    835                 x = x0 + xvec[j]
    836                 xnew = x - 2002.0
    837                 ref_points.append ([xnew,ynew]) #Relative point values
    838 
    839         assert num.allclose(points, ref_points)
    840 
    841         assert num.allclose(elevation, ref_elevation)
    842 
    843         #Cleanup
    844         fid.close()
    845 
    846 
    847         os.remove(root + '.pts')
    848         os.remove(root + '.dem')
    849         os.remove(root + '.asc')
    850         os.remove(root + '.prj')
    851 
    852 
    853     def test_dem2pts_bounding_box_removeNullvalues_v2(self):
    854         """Test conversion from dem in ascii format to native NetCDF format
    855         """
    856 
    857         import time, os
    858         from Scientific.IO.NetCDF import NetCDFFile
    859 
    860         #Write test asc file
    861         root = 'demtest'
    862 
    863         filename = root+'.asc'
    864         fid = open(filename, 'w')
    865         fid.write("""ncols         10
    866 nrows         10
    867 xllcorner     2000
    868 yllcorner     3000
    869 cellsize      1
    870 NODATA_value  -9999
    871 """)
    872         #Create linear function
    873         ref_points = []
    874         ref_elevation = []
    875         x0 = 2000
    876         y = 3010
    877         yvec = range(10)
    878         xvec = range(10)
    879         #z = range(100)
    880         z = num.zeros(100, num.int)     #array default#
    881         NODATA_value = -9999
    882         count = -1
    883         for i in range(10):
    884             y = y - 1
    885             for j in range(10):
    886                 x = x0 + xvec[j]
    887                 ref_points.append ([x,y])
    888                 count += 1
    889                 z[count] = (4*i - 3*j)%13
    890                 if j == 4: z[count] = NODATA_value #column inside clipping region
    891                 if j == 8: z[count] = NODATA_value #column outside clipping region
    892                 if i == 9: z[count] = NODATA_value #row outside clipping region
    893                 if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region
    894                 ref_elevation.append( z[count] )
    895                 fid.write('%f ' %z[count])
    896             fid.write('\n')
    897 
    898         fid.close()
    899 
    900         #print 'sending elev', ref_elevation
    901 
    902         #Write prj file with metadata
    903         metafilename = root+'.prj'
    904         fid = open(metafilename, 'w')
    905 
    906 
    907         fid.write("""Projection UTM
    908 Zone 56
    909 Datum WGS84
    910 Zunits NO
    911 Units METERS
    912 Spheroid WGS84
    913 Xshift 0.0000000000
    914 Yshift 10000000.0000000000
    915 Parameters
    916 """)
    917         fid.close()
    918 
    919         #Convert to NetCDF pts
    920         convert_dem_from_ascii2netcdf(root)
    921         dem2pts(root, easting_min=2002.0, easting_max=2007.0,
    922                 northing_min=3003.0, northing_max=3006.0,
    923                 verbose=self.verbose)
    924 
    925         #Check contents
    926         #Get NetCDF
    927         fid = NetCDFFile(root+'.pts', netcdf_mode_r)
    928 
    929         # Get the variables
    930         #print fid.variables.keys()
    931         points = fid.variables['points']
    932         elevation = fid.variables['elevation']
    933 
    934         #Check values
    935         assert fid.xllcorner[0] == 2002.0
    936         assert fid.yllcorner[0] == 3003.0
    937 
    938         #create new reference points
    939         newz = num.zeros(19, num.int)       #array default#
    940         newz[0:2] = ref_elevation[32:34]
    941         newz[2:5] = ref_elevation[35:38]
    942         newz[5:7] = ref_elevation[42:44]
    943         newz[7] = ref_elevation[45]
    944         newz[8] = ref_elevation[47]
    945         newz[9:11] = ref_elevation[52:54]
    946         newz[11:14] = ref_elevation[55:58]
    947         newz[14:16] = ref_elevation[62:64]
    948         newz[16:19] = ref_elevation[65:68]
    949 
    950 
    951         ref_elevation = newz
    952         ref_points = []
    953         new_ref_points = []
    954         x0 = 2002
    955         y = 3007
    956         yvec = range(4)
    957         xvec = range(6)
    958         for i in range(4):
    959             y = y - 1
    960             ynew = y - 3003.0
    961             for j in range(6):
    962                 x = x0 + xvec[j]
    963                 xnew = x - 2002.0
    964                 if j <> 2 and (i<>1 or j<>4):
    965                     ref_points.append([x,y])
    966                     new_ref_points.append ([xnew,ynew])
    967 
    968 
    969         assert num.allclose(points, new_ref_points)
    970         assert num.allclose(elevation, ref_elevation)
    971 
    972         #Cleanup
    973         fid.close()
    974 
    975 
    976         os.remove(root + '.pts')
    977         os.remove(root + '.dem')
    978         os.remove(root + '.asc')
    979         os.remove(root + '.prj')
    980 
    981 
    982     def test_dem2pts_bounding_box_removeNullvalues_v3(self):
    983         """Test conversion from dem in ascii format to native NetCDF format
    984         Check missing values on clipping boundary
    985         """
    986 
    987         import time, os
    988         from Scientific.IO.NetCDF import NetCDFFile
    989 
    990         #Write test asc file
    991         root = 'demtest'
    992 
    993         filename = root+'.asc'
    994         fid = open(filename, 'w')
    995         fid.write("""ncols         10
    996 nrows         10
    997 xllcorner     2000
    998 yllcorner     3000
    999 cellsize      1
    1000 NODATA_value  -9999
    1001 """)
    1002         #Create linear function
    1003         ref_points = []
    1004         ref_elevation = []
    1005         x0 = 2000
    1006         y = 3010
    1007         yvec = range(10)
    1008         xvec = range(10)
    1009         #z = range(100)
    1010         z = num.zeros(100, num.int)     #array default#
    1011         NODATA_value = -9999
    1012         count = -1
    1013         for i in range(10):
    1014             y = y - 1
    1015             for j in range(10):
    1016                 x = x0 + xvec[j]
    1017                 ref_points.append ([x,y])
    1018                 count += 1
    1019                 z[count] = (4*i - 3*j)%13
    1020                 if j == 4: z[count] = NODATA_value #column inside clipping region
    1021                 if j == 8: z[count] = NODATA_value #column outside clipping region
    1022                 if i == 6: z[count] = NODATA_value #row on clipping boundary
    1023                 if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region
    1024                 ref_elevation.append( z[count] )
    1025                 fid.write('%f ' %z[count])
    1026             fid.write('\n')
    1027 
    1028         fid.close()
    1029 
    1030         #print 'sending elev', ref_elevation
    1031 
    1032         #Write prj file with metadata
    1033         metafilename = root+'.prj'
    1034         fid = open(metafilename, 'w')
    1035 
    1036 
    1037         fid.write("""Projection UTM
    1038 Zone 56
    1039 Datum WGS84
    1040 Zunits NO
    1041 Units METERS
    1042 Spheroid WGS84
    1043 Xshift 0.0000000000
    1044 Yshift 10000000.0000000000
    1045 Parameters
    1046 """)
    1047         fid.close()
    1048 
    1049         #Convert to NetCDF pts
    1050         convert_dem_from_ascii2netcdf(root)
    1051         dem2pts(root, easting_min=2002.0, easting_max=2007.0,
    1052                 northing_min=3003.0, northing_max=3006.0,
    1053                 verbose=self.verbose)
    1054 
    1055         #Check contents
    1056         #Get NetCDF
    1057         fid = NetCDFFile(root+'.pts', netcdf_mode_r)
    1058 
    1059         # Get the variables
    1060         #print fid.variables.keys()
    1061         points = fid.variables['points']
    1062         elevation = fid.variables['elevation']
    1063 
    1064         #Check values
    1065         assert fid.xllcorner[0] == 2002.0
    1066         assert fid.yllcorner[0] == 3003.0
    1067 
    1068         #create new reference points
    1069         newz = num.zeros(14, num.int)       #array default#
    1070         newz[0:2] = ref_elevation[32:34]
    1071         newz[2:5] = ref_elevation[35:38]
    1072         newz[5:7] = ref_elevation[42:44]
    1073         newz[7] = ref_elevation[45]
    1074         newz[8] = ref_elevation[47]
    1075         newz[9:11] = ref_elevation[52:54]
    1076         newz[11:14] = ref_elevation[55:58]
    1077 
    1078 
    1079 
    1080         ref_elevation = newz
    1081         ref_points = []
    1082         new_ref_points = []
    1083         x0 = 2002
    1084         y = 3007
    1085         yvec = range(4)
    1086         xvec = range(6)
    1087         for i in range(4):
    1088             y = y - 1
    1089             ynew = y - 3003.0
    1090             for j in range(6):
    1091                 x = x0 + xvec[j]
    1092                 xnew = x - 2002.0
    1093                 if j <> 2 and (i<>1 or j<>4) and i<>3:
    1094                     ref_points.append([x,y])
    1095                     new_ref_points.append ([xnew,ynew])
    1096 
    1097 
    1098         #print points[:],points[:].shape
    1099         #print new_ref_points, len(new_ref_points)
    1100 
    1101         assert num.allclose(elevation, ref_elevation)
    1102         assert num.allclose(points, new_ref_points)
    1103 
    1104 
    1105         #Cleanup
    1106         fid.close()
    1107 
    1108 
    1109         os.remove(root + '.pts')
    1110         os.remove(root + '.dem')
    1111         os.remove(root + '.asc')
    1112         os.remove(root + '.prj')
    1113737
    1114738
Note: See TracChangeset for help on using the changeset viewer.