I'm trying to implement a script to project an image having the two matrices of Lat and Lon points. I've followed the steps proposed in: Unable to warp HDF5 files and everything is fine. I create the three .vrt files and than by using:
gdalwarp -geoloc -t_srs EPSG:4326 img.vrt output.tif
the final image is correctly projected. Also, just for visualization purposes the following code does is job:
plt.pcolormesh(lon,lat,img, cmap='rainbow') plt.show()
Now, I want to create a script without having to create the .vrt files every time. So I started following Writing numpy array to raster file but I'm not able to obtain the desired projection (Actually, the output image is similar to the input one, just horizontally translate of some pixels). Here is the key part of code:
ny = img.shape nx = img.shape xmin, ymin, xmax, ymax = [lon.min(), lat.min(), lon.max(), lat.max()] xres = (xmax - xmin) / float(nx) yres = (ymax - ymin) / float(ny) geotransform = (xmin, xres, 0, ymin, 0, -yres) # create the 1-band raster file dst_ds = gdal.GetDriverByName('GTiff').Create('myGeoTIFF.tif', ny, nx, 1, gdal.GDT_Byte) dst_ds.SetGeoTransform(geotransform) # specify coords srs = osr.SpatialReference() # establish encoding srs.ImportFromEPSG(4326) # WGS84 lat/long dst_ds.SetProjection(srs.ExportToWkt()) # export coords to file dst_ds.GetRasterBand(1).WriteArray(img) # write r-band to the raster dst_ds.FlushCache() # write to disk dst_ds = None
I think I'm missing something when I create the geotransform matrix because in this way I introduce only one corner of the image without giving any info about the transformation.
How do I solve the problem?