在将2D Numpy数组转换为GeoTIFF-fi时如何处理nan

2024-10-03 11:23:11 发布

您现在位置:Python中文网/ 问答频道 /正文

我想把一个二维的Numpy数组转换成一个3波段(RGB)GeoTIFF文件。 您将在下面找到我的一段代码,这就是我要做的:

  1. 获取数组“topo”(由nan值组成)
  2. 将numpy数组转换为RGB(如here
  3. 列表项

我真的不知道如何处理南值。。。在步骤2和3中。 你能帮帮我吗? 谢谢

    import numpy as np
    import matplotlib.pyplot as plt
    from osgeo import gdal, osr
    # My image array
    topo = self.topo # topo is a numpy array
    lat = self.Y
    lon = self.X
    # normalize array values to [0,1]
    cmin = np.nanmin(topo)
    cmax = np.nanmax(topo)
    my_cm = plt.get_cmap('jet')
    normed_topo = (topo - cmin) / (cmax - cmin)
    mapped_topo = my_cm(normed_topo)                # => REMOVES MY NAN VALUES
    # which will give a MxNx4 array mapped between 0 and 1,
    mapped_topou8 = (255 * mapped_topo).astype('float32')
    r, g, b = mapped_topou8[:, :, 1], mapped_topou8[:, :, 2], mapped_topou8[:, :, 3]

    # HOW TO DEAL WITH NAN ?
    # Find indices that I need to replace
    inds = np.where(np.isnan(topo))
    # Place column means in the indices. Align the arrays using take
    r[inds] = np.nan
    g[inds] = np.nan
    b[inds] = np.nan
    # Convert to uint8 - IS IT USEFUL?
    r = r.astype('uint8')
    g = g.astype('uint8')
    b = b.astype('uint8')

    # geotransform needs coordinates of one corner and the resolution of the file
    xmin, ymin, xmax, ymax = [lon.min(), lat.min(), lon.max(), lat.max()]
    nrows, ncols = np.shape(topo)
    xres = (xmax - xmin) / float(ncols)
    yres = (ymax - ymin) / float(nrows)
    geotransform = (xmin, xres, 0, ymax, 0, -yres)
    # That's (top left x, w-e pixel resolution, rotation (0 if North is up),
    #         top left y, rotation (0 if North is up), n-s pixel resolution)

    # Create the 3-band raster file
    output_raster = gdal.GetDriverByName('GTiff').Create(sspath, ncols, nrows, 3, gdal.GDT_Float32)
    output_raster.SetGeoTransform(geotransform)  # Specify its coordinates
    srs = osr.SpatialReference()                 # Establish its coordinate encoding
    srs.ImportFromEPSG(3857)                     # This one specifies WGS84 lat long

    output_raster.SetProjection(srs.ExportToWkt())  # Exports the coordinate system to the file
   output_raster.GetRasterBand(1).WriteArray(r)    # Write r-band to the raster
    output_raster.GetRasterBand(2).WriteArray(g)    # Write g-band to the raster
    output_raster.GetRasterBand(3).WriteArray(b)    # Write b-band to the raster
    output_raster.FlushCache()                      # Write to disk.

Tags: thetooutputbandnpnanarraylat