Python crs.CRS类代码示例

def test_from_string():
    wgs84_crs = CRS.from_string('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
    assert wgs84_crs.to_dict() == {'init': 'epsg:4326'}

    # Make sure this doesn't get handled using the from_epsg() even though 'epsg' is in the string
    epsg_init_crs = CRS.from_string('+init=epsg:26911')
    assert epsg_init_crs.to_dict() == {'init': 'epsg:26911'}

def test_from_string():
    wgs84_crs = CRS.from_string('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
    assert wgs84_crs.to_dict() == {'no_defs': True, 'ellps': 'WGS84', 'datum': 'WGS84', 'proj': 'longlat'}

    # Make sure this doesn't get handled using the from_epsg() even though 'epsg' is in the string
    epsg_init_crs = CRS.from_string('+units=m +init=epsg:26911 +no_defs=True')
    assert epsg_init_crs.to_dict() == {'units': 'm', 'init': 'epsg:26911', 'no_defs': True}

def test_from_epsg_string():
    crs_dict = CRS.from_string('epsg:4326')
    assert crs_dict['init'].lower() == 'epsg:4326'

    # Test with invalid EPSG code
    with pytest.raises(ValueError):
        assert CRS.from_string('epsg:xyz')

def test_is_projected():
    assert CRS({'init': 'EPSG:3857'}).is_projected is True

    lcc_crs = CRS.from_string('+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0')
    assert CRS(lcc_crs).is_projected is True

    wgs84_crs = CRS.from_string('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
    assert CRS(wgs84_crs).is_projected is False

def test_from_proj4_json():
    json_str = '{"proj": "longlat", "ellps": "WGS84", "datum": "WGS84"}'
    crs_dict = CRS.from_string(json_str)
    assert crs_dict == json.loads(json_str)

    # Test with invalid JSON code
    with pytest.raises(ValueError):
        assert CRS.from_string('{foo: bar}')

def test_issue_1446():
    """Confirm resolution of #1446"""
    g = transform_geom(
        CRS.from_epsg(4326),
        CRS.from_epsg(32610),
        {"type": "Point", "coordinates": (-122.51403808499907, 38.06106733107932)},
    )
    assert round(g["coordinates"][0], 1) == 542630.9
    assert round(g["coordinates"][1], 1) == 4212702.1

def test_is_geographic_from_string():
    wgs84_crs = CRS.from_string('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
    assert wgs84_crs.is_geographic is True

    nad27_crs = CRS.from_string('+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs')
    assert nad27_crs.is_geographic is True

    lcc_crs = CRS.from_string('+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0')
    assert lcc_crs.is_geographic is False

def test_bare_parameters():
    """ Make sure that bare parameters (e.g., no_defs) are handled properly,
    even if they come in with key=True.  This covers interaction with pyproj,
    which makes presents bare parameters as key=<bool>."""

    # Example produced by pyproj
    crs_dict = CRS.from_string('+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0')
    assert crs_dict.get('no_defs', False) is True

    crs_dict = CRS.from_string('+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=False +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0')
    assert crs_dict.get('no_defs', True) is False

def test_symmetric_proj4(tmpdir):
    """ Test writing and reading proj4 string as attribute of variable """

    ds = Dataset(str(tmpdir.join('test.nc')), 'w')
    proj4 = '+proj=stere +units=m +datum=WGS84 +lat_ts=60 +lat_0=90 +lon_0=263 +lat_1=60 +x_0=3475000 +y_0=7475000'
    ds.createVariable('data', 'S1')
    set_crs(ds, 'data', Proj(proj4), set_proj4_att=True)
    out_proj4 = get_crs(ds, 'data')

    out_data = CRS.from_string(out_proj4).to_dict()

    assert len(out_data) == 9  # There should be 9 parameters
    assert CRS.from_string(proj4).to_dict() == out_data

def test_is_same_crs():
    crs1 = CRS({'init': 'EPSG:4326'})
    crs2 = CRS({'init': 'EPSG:3857'})

    assert crs1 == crs1
    assert crs1 != crs2

    wgs84_crs = CRS.from_string('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
    assert crs1 == wgs84_crs

    # Make sure that same projection with different parameter are not equal
    lcc_crs1 = CRS.from_string('+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0')
    lcc_crs2 = CRS.from_string('+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=45 +lat_0=0')
    assert lcc_crs1 != lcc_crs2

def test_get_crs(tmpdir):
    """ Test reading proj4 string from CF convention parameters """

    ds = Dataset(str(tmpdir.join('test.nc')), 'w')
    data_var = ds.createVariable('data', 'S1')
    data_var.setncattr('grid_mapping', 'crs_Lambert')
    crs_var = ds.createVariable('crs_Lambert', 'S1')

    in_proj4 = '+proj=lcc +units=m +lat_1=30 +lat_2=60 +lat_0=47.5 +lon_0=-97 +x_0=3825000 +y_0=3200000'

    # These parameters match the above proj4 string
    ncatts = dict()
    ncatts['grid_mapping_name'] = 'lambert_conformal_conic'
    ncatts['latitude_of_projection_origin'] = 47.5
    ncatts['longitude_of_central_meridian'] = -97
    ncatts['standard_parallel'] = [30, 60]
    ncatts['false_northing'] = 3200000
    ncatts['false_easting'] = 3825000
    set_ncattrs(crs_var, ncatts)

    out_proj4 = get_crs(ds, 'data')
    assert out_proj4 is not None

    out_data = CRS.from_string(out_proj4).to_dict()

    assert len(out_data) == 8  # There should be 8 parameters
    assert CRS.from_string(in_proj4).to_dict() == out_data

    # Test WGS84 lat/long
    data_var = ds.createVariable('data2', 'S1')
    data_var.setncattr('grid_mapping', 'crs_latlong')
    crs_var = ds.createVariable('crs_latlong', 'S1')

    in_proj4 = '+proj=latlong +a={0} +rf={1}'.format(pj_ellps['WGS84']['a'], pj_ellps['WGS84']['rf'])

    # These parameters match the above proj4 string
    ncatts = dict()
    ncatts['grid_mapping_name'] = 'latitude_longitude'
    ncatts['semi_major_axis'] = 6378137.0
    ncatts['inverse_flattening'] = 298.257223563
    set_ncattrs(crs_var, ncatts)

    out_proj4 = get_crs(ds, 'data2')
    assert out_proj4 is not None

    out_data = CRS.from_string(out_proj4).to_dict()

    assert len(out_data) == 4  # There should be 4 parameters
    # Note: pyproj adds units=m even for latlong, which is incorrect but not our problem
    assert CRS.from_string(in_proj4 + ' +units=m').to_dict() == out_data

def main():
    samplefile = r'bxk1-d-ck.idf'
    tiffile = samplefile.replace('.idf', '.geotiff')
    dtype = rasterio.float64
    driver = 'AAIGrid'
    crs = CRS.from_epsg(28992)

    # read data from idf file
    idffile = idfpy.IdfFile(filepath=samplefile, mode='rb')
    geotransform = idffile.geotransform
    height = idffile.header['nrow']
    width = idffile.header['ncol']
    nodata = idffile.header['nodata']
    transform = Affine.from_gdal(*geotransform)

    # write data from idf file to geotiff with rasterio
    profile = {
        'width': width,
        'height': height,
        'count': 1,
        'dtype': dtype,
        'driver': driver,
        'crs': crs,
        'transform': transform,
        'nodata': nodata,
    }

    # the default profile would be sufficient for the example, however the profile dict shows how to make the export
    # profile
    idffile.to_raster(tiffile, **profile)

    def index(self, *args, **kwargs):
        target_srs = self.crs
        # Annotations are always in WGS84
        source_srs = CRS.from_string("EPSG:4326")

        args = transform_coordinates(source_srs, target_srs,
                                     [args[0]], [args[1]])
        return self.dataset.index(*args, **kwargs)

    def __init__(self, ul, crs, res, size, desc=None):
        self.ul = ul
        if isinstance(crs, six.string_types):
            self.crs = CRS.from_string(crs)
        elif isinstance(crs, int):
            self.crs = CRS.from_epsg(crs)
        else:
            self.crs = crs

        if not self.crs.is_valid:
            raise ValueError('Could not parse coordinate reference system '
                             'string to a valid projection ({})'.format(crs))

        self.crs_str = self.crs.to_string()
        self.res = res
        self.size = size
        self.desc = desc or 'unnamed'
        self._tiles = {}

    def get_band_ix(self, indexes, x, y):
        # Reproject to native data coordinates
        target_srs = self.crs
        # Annotations are always in WGS84
        source_srs = CRS.from_string("EPSG:4326")

        transformed_x, transformed_y = transform_coordinates(source_srs,
                                                             target_srs,
                                                             [x], [y])
        return list(self.dataset.sample([(transformed_x, transformed_y)],
                                        indexes=indexes))[0]