Changeset 4180


Ignore:
Timestamp:
Jan 17, 2007, 2:17:22 PM (18 years ago)
Author:
duncan
Message:

.txt files can now be in lats and longs.

Location:
anuga_core/source/anuga
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/source/anuga/coordinate_transforms/geo_reference.py

    r3614 r4180  
    228228
    229229    def reconcile_zones(self, other):
    230 
     230        if other is None:
     231            other = Geo_reference()
    231232        if self.zone == other.zone or\
    232233               self.zone == DEFAULT_ZONE and other.zone == DEFAULT_ZONE:

  • anuga_core/source/anuga/geospatial_data/geospatial_data.py

    r4179 r4180  
    77from types import DictType
    88from warnings import warn
    9 
     9from string import lower
    1010from Numeric import concatenate, array, Float, shape, reshape, ravel, take, \
    1111                        size, shape
    … …  
    1818from anuga.coordinate_transforms.geo_reference import Geo_reference, TitleError
    1919from anuga.coordinate_transforms.redfearn import convert_from_latlon_to_utm
     20from anuga.utilities.anuga_exceptions import ANUGAError
    2021
    2122MAX_READ_LINES = 500       
    … …  
    7576        have dimensions "points" etc.
    7677        .txt file is a comma seperated file with x, y and attribute
    77         data.
    78         the first line has the titles of the columns.  The first two columns
    79         are assumed to be the x and y, and the title names acually used are
    80         ignored.
    81 
    82         The
     78        data.
     79       
     80        The first line has the titles of the columns.  The first two
     81        column titles are checked to see if they start with lat or
     82        long (not case sensitive).  If so the data is assumed to be
     83        latitude and longitude, in decimal format and converted to
     84        UTM.  Otherwise the first two columns are assumed to be the x
     85        and y, and the title names acually used are ignored.
     86 
    8387       
    8488        The format for a .txt file is:
    … …  
    141145                   points_are_lats_longs:
    142146                data_points, geo_reference =  \
    143                              self._set_using_lat_long(latitudes=latitudes,
     147                             _set_using_lat_long(latitudes=latitudes,
    144148                                  longitudes=longitudes,
    145149                                  geo_reference=geo_reference,
    … …  
    324328        return clipped_G
    325329
    326    
    327     def _set_using_lat_long(self,
    328                             latitudes,
    329                             longitudes,
    330                             geo_reference,
    331                             data_points,
    332                             points_are_lats_longs):
    333        
    334         if geo_reference is not None:
    335             msg = """A georeference is specified yet latitude and longitude
    336             are also specified!"""
    337             raise ValueError, msg
    338        
    339         if data_points is not None and not points_are_lats_longs:
    340             msg = """Data points are specified yet latitude and
    341             longitude are also specified!"""
    342             raise ValueError, msg
    343        
    344         if points_are_lats_longs:
    345             if data_points is None:
    346                 msg = """Data points are not specified !"""
    347                 raise ValueError, msg
    348             lats_longs = ensure_numeric(data_points)
    349             latitudes = ravel(lats_longs[:,0:1])
    350             longitudes = ravel(lats_longs[:,1:])
    351            
    352         if latitudes is None and longitudes is None:
    353             msg = """Latitudes and Longitudes are not."""
    354             raise ValueError, msg
    355        
    356         if latitudes is None:
    357             msg = """Longitudes are specified yet latitudes aren't!"""
    358             raise ValueError, msg
    359        
    360         if longitudes is None:
    361             msg = """Latitudes are specified yet longitudes aren't!"""
    362             raise ValueError, msg
    363        
    364         data_points, zone  = convert_from_latlon_to_utm(latitudes=latitudes,
    365                                                         longitudes=longitudes)
    366        
    367         return data_points, Geo_reference(zone=zone)
    368330   
    369331    def get_geo_reference(self):
    … …  
    738700            self.header, self.file_pointer = \
    739701                         _read_csv_file_header(file_pointer)
    740 
     702            self.blocking_georef = None # Used for reconciling zones
    741703        if self.max_read_lines is None:
    742704            self.max_read_lines = MAX_READ_LINES
    … …  
    783745        else:
    784746            try:
    785                 pointlist, att_dict, self.file_pointer = \
     747                pointlist, att_dict, geo_ref, self.file_pointer = \
    786748                   _read_csv_file_blocking( self.file_pointer,
    787749                                            self.header[:],
    788750                                            max_read_lines=self.max_read_lines,
    789751                                            verbose=self.verbose)
    790                 geo = Geospatial_data(pointlist, att_dict)
     752
     753                # Check that the zones haven't changed.
     754                if geo_ref is not None:
     755                    geo_ref.reconcile_zones(self.blocking_georef)
     756                    self.blocking_georef = geo_ref
     757                elif self.blocking_georef is not None:
     758                   
     759                    msg = 'Geo reference given, then not given.'
     760                    msg += ' This should not happen.'
     761                    raise ValueError, msg
     762                geo = Geospatial_data(pointlist, att_dict, geo_ref)
    791763            except StopIteration:
    792764                self.file_pointer.close()
    … …  
    794766                del self.file_pointer
    795767                raise StopIteration
     768            except ANUGAError:
     769                self.file_pointer.close()
     770                del self.header
     771                del self.file_pointer
     772                raise
    796773        return geo
    797774
     775def _set_using_lat_long(latitudes,
     776                        longitudes,
     777                        geo_reference,
     778                        data_points,
     779                        points_are_lats_longs):
     780    """
     781    if the points has lat long info, assume it is in (lat, long) order.
     782    """
     783   
     784    if geo_reference is not None:
     785        msg = """A georeference is specified yet latitude and longitude
     786        are also specified!"""
     787        raise ValueError, msg
     788   
     789    if data_points is not None and not points_are_lats_longs:
     790        msg = """Data points are specified yet latitude and
     791        longitude are also specified!"""
     792        raise ValueError, msg
     793   
     794    if points_are_lats_longs:
     795        if data_points is None:
     796            msg = """Data points are not specified !"""
     797            raise ValueError, msg
     798        lats_longs = ensure_numeric(data_points)
     799        latitudes = ravel(lats_longs[:,0:1])
     800        longitudes = ravel(lats_longs[:,1:])
     801       
     802    if latitudes is None and longitudes is None:
     803        msg = """Latitudes and Longitudes are not."""
     804        raise ValueError, msg
     805   
     806    if latitudes is None:
     807        msg = """Longitudes are specified yet latitudes aren't!"""
     808        raise ValueError, msg
     809   
     810    if longitudes is None:
     811        msg = """Latitudes are specified yet longitudes aren't!"""
     812        raise ValueError, msg
     813   
     814    data_points, zone  = convert_from_latlon_to_utm(latitudes=latitudes,
     815                                                    longitudes=longitudes)
     816   
     817    return data_points, Geo_reference(zone=zone)
     818   
    798819def _read_pts_file(file_name, verbose=False):
    799820    """Read .pts NetCDF file
    … …  
    858879    file_pointer = open(file_name)
    859880    header, file_pointer = _read_csv_file_header(file_pointer)
    860     pointlist, att_dict,file_pointer  = _read_csv_file_blocking( \
     881    try:
     882        pointlist, att_dict, geo_ref, file_pointer = \
     883                   _read_csv_file_blocking( \
    861884                file_pointer,
    862885                header,
    863886                max_read_lines=1e30) #If the file is bigger that this, block..
    864        
     887    except ANUGAError:
     888        file_pointer.close()
     889        raise
    865890    file_pointer.close()
    866     return pointlist, att_dict, None   
     891    return pointlist, att_dict, geo_ref   
    867892
    868893CSV_DELIMITER = ','
    … …  
    894919    #This is to remove the x and y headers.
    895920    header = header[:]
    896     header.pop(0)
    897     header.pop(0)
     921    x_header = header.pop(0)
     922    y_header = header.pop(0)
    898923   
    899924    read_lines = 0
    … …  
    936961    for key in att_dict.keys():
    937962        att_dict[key] = array(att_dict[key]).astype(Float)
    938        
    939     return pointlist, att_dict,file_pointer
     963
     964    #Do stuff here so the info is in lat's and longs
     965    geo_ref = None
     966    x_header = lower(x_header[:3])
     967    y_header = lower(y_header[:3])
     968    if (x_header == 'lon' or  x_header == 'lat') and \
     969       (y_header == 'lon' or  y_header == 'lat'):
     970        if x_header == 'lon':
     971            longitudes = ravel(pointlist[:,0:1])
     972            latitudes = ravel(pointlist[:,1:])
     973        else:
     974            latitudes = ravel(pointlist[:,0:1])
     975            longitudes = ravel(pointlist[:,1:])
     976       
     977        pointlist, geo_ref = _set_using_lat_long(latitudes,
     978                                                 longitudes,
     979                                                 geo_reference=None,
     980                                                 data_points=None,
     981                                                 points_are_lats_longs=False)
     982    return pointlist, att_dict, geo_ref, file_pointer
    940983
    941984def _read_pts_file_header(fid, verbose=False):

  • anuga_core/source/anuga/geospatial_data/test_geospatial_data.py

    r4179 r4180  
    1111from anuga.coordinate_transforms.geo_reference import Geo_reference, TitleError
    1212from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
    13 
     13from anuga.utilities.anuga_exceptions import ANUGAError
    1414# Ignore these warnings, since we still want to test .xya code.
    1515import warnings
    … …  
    679679        assert allclose(results.get_attributes(attribute_name='speed'),
    680680                        [0.0, 10.0, 40.0])
    681        
    682     def not_test_loadcsv(self):
    683         """ not_test_loadcsv(self):
    684         comma delimited
    685         """
     681
     682    def test_load_csv(self):
     683       
     684        import os
     685        import tempfile
     686       
    686687        fileName = tempfile.mktemp(".csv")
    687688        file = open(fileName,"w")
    688         file.write("longitude,latitude,z \n\
    689 -35.3149601,150.9198238,452.688000\n\
    690 -35.3149791,150.9209232,459.126000\n\
    691 -35.3149980,150.9220226,465.613000\n")
    692         file.close()
     689        file.write("lat,long,elevation speed \n\
     6901.0 0.0 10.0 0.0\n\
     6910.0 1.0 0.0 10.0\n\
     6921.0 0.0 10.4 40.0\n")
     693        file.close()
     694        #print fileName
    693695        results = Geospatial_data(fileName, delimiter=',')
    694696        os.remove(fileName)
    … …  
    700702        assert allclose(results.get_attributes(attribute_name='speed'),
    701703                        [0.0, 10.0, 40.0])
    702 
     704       
     705       
    703706    def test_load_xya(self):
    704707        """ test_load_xya(self):
    … …  
    947950  ###################### .CSV ##############################
    948951
     952    def test_load_csv_lat_long_bad_blocking(self):
     953       
     954        fileName = tempfile.mktemp(".csv")
     955        file = open(fileName,"w")
     956        file.write("Lati,LONG,z \n\
     957-25.0,180.0,452.688000\n\
     958-34,150.0,459.126000\n")
     959        file.close()
     960       
     961        results = Geospatial_data(fileName, max_read_lines=1,
     962                                  load_file_now=False)
     963       
     964        try:
     965            for i in results:
     966                pass
     967        except ANUGAError:
     968            pass
     969        else:
     970            msg = 'Different zones in Geo references not caught.'
     971            raise msg       
     972       
     973        os.remove(fileName)
     974       
    949975    def test_load_csv(self):
    950         """ test_load_csv(self):
    951         space delimited
    952         """
    953         import os
    954        
     976       
    955977        fileName = tempfile.mktemp(".txt")
    956978        file = open(fileName,"w")
    … …  
    17591781
    17601782
     1783    def test_load_csv_lat_long(self):
     1784        """
     1785        comma delimited
     1786
     1787        """
     1788        fileName = tempfile.mktemp(".csv")
     1789        file = open(fileName,"w")
     1790        file.write("long,lat,z \n\
     1791150.916666667,-34.50,452.688000\n\
     1792150.0,-34,459.126000\n")
     1793        file.close()
     1794        results = Geospatial_data(fileName, delimiter=',')
     1795        os.remove(fileName)
     1796        points = results.get_data_points()
     1797       
     1798        assert allclose(points[0][0], 308728.009)
     1799        assert allclose(points[0][1], 6180432.601)
     1800        assert allclose(points[1][0],  222908.705)
     1801        assert allclose(points[1][1], 6233785.284)
     1802       
     1803     
     1804    def test_load_csv_lat_longII(self):
     1805        """
     1806        comma delimited
     1807
     1808        """
     1809        fileName = tempfile.mktemp(".csv")
     1810        file = open(fileName,"w")
     1811        file.write("Lati,LONG,z \n\
     1812-34.50,150.916666667,452.688000\n\
     1813-34,150.0,459.126000\n")
     1814        file.close()
     1815        results = Geospatial_data(fileName, delimiter=',')
     1816        os.remove(fileName)
     1817        points = results.get_data_points()
     1818       
     1819        assert allclose(points[0][0], 308728.009)
     1820        assert allclose(points[0][1], 6180432.601)
     1821        assert allclose(points[1][0],  222908.705)
     1822        assert allclose(points[1][1], 6233785.284)
     1823
     1824         
     1825    def test_load_csv_lat_long_bad(self):
     1826        """
     1827        comma delimited
     1828
     1829        """
     1830        fileName = tempfile.mktemp(".csv")
     1831        file = open(fileName,"w")
     1832        file.write("Lati,LONG,z \n\
     1833-25.0,180.0,452.688000\n\
     1834-34,150.0,459.126000\n")
     1835        file.close()
     1836        try:
     1837            results = Geospatial_data(fileName, delimiter=',')
     1838        except ANUGAError:
     1839            pass
     1840        else:
     1841            msg = 'Different zones in Geo references not caught.'
     1842            raise msg       
     1843       
     1844        os.remove(fileName)
    17611845       
    17621846    def test_lat_long(self):
Note: See TracChangeset for help on using the changeset viewer.