3

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[0]
nx = img.shape[1]

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?

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
Thomas
  • 81
  • 3
  • Geotransform is just a linear transformation with scaling, while gdalwarp can do more sophisticated (/and accurate) reprojection. Keep in mind that the Earth is round ... – AndreJ Sep 08 '17 at 11:38
  • Thanks. That exactly what I was thinking. I'll try to explore a bit more gdalwarp and see how to extract the correct transformation from the lat/lon matrices – Thomas Sep 08 '17 at 12:39
  • Take a look at the files metadata or the download website. In most cases the projection is documented there. – AndreJ Sep 09 '17 at 06:38
  • I did a small step back just to well understand the gdalwarp function. I've created with the tool of the image provider the geotiff image and I tried to apply to it the gdalwarp function in this way: gdalwarp -t_srs EPSG:4326 original.tif warp_ref.tif. But the resulting image is not projected. Do you have expirience with Sentinel 3 data? the metadata file says: srsName="http://www.opengis.net/def/crs/EPSG/0/4326"> – Thomas Sep 10 '17 at 09:40
  • Can you point to the dataset you are having trouble with? – AndreJ Sep 13 '17 at 06:11
  • Hi, I'm working with a standard OLCI sentinel 3 tile available on eumesat weebsite https://www.eumetsat.int/website/home/Data/DataDelivery/CopernicusOnlineDataAccess/index.html. – Thomas Sep 13 '17 at 11:05
  • I tried the Data Conversion Tool From SNAP on the altitude band and it creates a misplaced tif file :-( . I guess the problem are the lat and lon tables being integer values that need a scaling factor of 1e06. – AndreJ Sep 13 '17 at 14:32
  • any advancement on the geolocation of the OLCI products ? I'm working right now on the same issue, and namely the WRR and WST products, among others, and I'd like to exchange with you on various options .. (so far, I've explored both the sapty library and a bash scripting based on .vrt but nothing really satisfactory ...) – Marco Mar 29 '18 at 14:51
  • This does not provide an answer to the question. Once you have sufficient reputation you will be able to comment on any post; instead, provide answers that don't require clarification from the asker. - From Review – Dan C Mar 29 '18 at 15:21

0 Answers0