I have a few NumPy arrays: lat, lon, and 6 data bands. I'd like to convert these to GeoTiff and then apply the gdalwarp function, based on a specific projection, to each band.
Is there a standard way to do this in Python?
4
-
possible dublicate: gis.stackexchange.com/questions/37238/…– ThomasCommented Jun 11, 2016 at 22:27
-
Possible, yes. A major difference is that my lat/lon resolutions aren't always uniform. Defining the geotransform is the key, I think.– rakbarCommented Jun 11, 2016 at 23:04
-
So they don't form a complete grid? In this case you will either have to interpolate a surface from the observations, or be left with missing values. Any preference? Can you better characterize the presence of lat/lon observations for us?– Jeff GCommented Jun 12, 2016 at 14:34
-
1To give you more context: I have a series of HDF5 files, with swath grided data and lat/lon values. The coordinate description is " geodetic lat/lon of the center of each 1 km cell." I need to first apply some flagging and filtering operations to the data, then re-project to EPSG:6933. I was making virtual rasters (.vrt files) then using gdalwarp to do all this. But because of the flagging operations, things got a bit complicated, hence numpy arrays. The HDF5 don't have an specific geolocation meta data either. Hope this helps!– rakbarCommented Jun 12, 2016 at 22:11
Add a comment
|
1 Answer
Can be done fairly easy. First, create the dataset:
drv = gdal.GetDriverByName("GTiff")
ds = drv.Create("name.tif", width, height, 6, gdal.GDT_Float32)
You then need to set the geotransform using ds.SetGeoTransform where the argument is a six element tuple: (upper_left_x, x_resolution, x_skew, upper_left_y, y_skew, y_resolution)
Set the projection using ds.SetProjection, argument is the way string, then
Write the array for the first band:
ds.GetRasterBand(1).WriteArray(arr)
That should be enough to get you going. If the purpose is just to reproject the data, then you could look at using an in memory dataset
-
2Make sure to use ds.SetGeoTransform instead of ds.SetGeotransform– ArnJacCommented Feb 14, 2018 at 15:04