nimrod_to_nc.py 文件源码

def warp(nimrod_dataset):
    """
    Warps the data array into one that has longitude/latitude as axes an fits
    the EPSG:4326 spatial reference system. The original array has the srs
    EPSG:27700 (OSGB 1936).
    :param nimrod_dataset: dictionary containing the data from the NIMROD file
    """

    # http://gis.stackexchange.com/questions/139906/replicating-result-of-gdalwarp-using-gdal-python-bindings
    # Create synthetic data
    gtiff_drv = gdal.GetDriverByName('MEM')
    cols, rows = nimrod_dataset['cols'], nimrod_dataset['rows']
    raster = np.reshape(nimrod_dataset['data'], (cols, rows))
    raster = np.int16(raster)
    top_left = (nimrod_dataset['start_easting'], nimrod_dataset['start_northing'])
    pixel_height = nimrod_dataset['column_interval']
    pixel_width = nimrod_dataset['row_interval']

    src_srs = osr.SpatialReference()
    src_srs.ImportFromEPSG(27700)

    src_geotran = [top_left[0], pixel_width,           0,
                   top_left[1],     0,      -pixel_height]

    rows, cols = raster.shape
    src_ds = gtiff_drv.Create(
        'test_epsg3413.tif',
        cols, rows, 1,
        gdal.GDT_Int16)
    src_ds.SetGeoTransform(src_geotran)
    src_ds.SetProjection(src_srs.ExportToWkt())
    src_ds.GetRasterBand(1).WriteArray(raster)

    # Transform to EPSG: 4326
    dest_srs = osr.SpatialReference()
    dest_srs.ImportFromEPSG(4326)

    int_ds = gdal.AutoCreateWarpedVRT(src_ds, src_srs.ExportToWkt(), dest_srs.ExportToWkt())
    nimrod_dataset['data_warped'] = int_ds.GetRasterBand(1).ReadAsArray()
    nimrod_dataset['GeoTransform'] = int_ds.GetGeoTransform()

    src_ds = None
    int_ds = None

    return nimrod_dataset