5

I am trying to create an image from a numpy array and georefenrece it so that I can import it into qgis. I have the latitude and longitude in degrees for the top left corner of the array as well as the spacing for the grids in metres. My understanding is that I can use the SetGeoTransform function to set these values and then set the projection which, in my case is lcc. The problem is that the image lands on a completely different part of the canvas. For the

SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))

is it correct to set the origin in degrees and the width in metres?

Can you see what elese I am doing incorrectly?

THis is an example of the code:

import osr
import numpy
import gdal
format = "GTiff"
driver = gdal.GetDriverByName( format )
dst_ds = driver.Create( 'dst_filename.tif', 134, 164, 1, gdal.GDT_Float64)
dst_ds.SetGeoTransform( [ -10, 8000, 0, 45, 0, -8000 ] )
srs = osr.SpatialReference()
srs.ImportFromProj4('+proj=lcc +lat_1=45 +lat_2=45 +lat_0=45 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs')
dst_ds.SetProjection( srs.ExportToWkt() )
raster = numpy.zeros( (164, 134), dtype=numpy.uint64)    
dst_ds.GetRasterBand(1).WriteArray(array)
# Once we're done, close properly the dataset
dst_ds = None
Scott
  • 81
  • 1
  • 3
  • 1
    your geotransform should have units that match your projection. also, your origin for the raster is not the same as the origin of our numpy array...bottom left is the raster origin which means [164,0] for your numpy...matching those could shift the display. – user1269942 Nov 03 '15 at 17:46
  • @user1269942 how do you set the bottom left to 164,0? thanks – Scott Nov 03 '15 at 17:53
  • it depends...I'm not saying that you're doing it wrong...it's hard to tell. But if 0,0 in your numpy array is the point that matches -10, 45 in your projection, then you need to assign your geotransform's orgin to be (-10, 45 - 164 * distperpixel). If I were you, I'd try using a projection using meters...it would help to get you up to speed with the various number a little easier than degrees. Since your pixel with is 8000m it would keep it nicely aligned. – user1269942 Nov 03 '15 at 18:04
  • @user1269942 thanks for this however the datasets that I was given have all the corners in degrees and the grid spacings were defined in meters. Would you know how I can convert these lat lons into metres? this is the projection that was used: +proj=lcc +lat_1=45 +lat_2=45 +lat_0=45 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs. Thanks – Scott Nov 03 '15 at 18:40
  • actually, that may be fine. excuse my oversight, I didn't see the units=m in your projection string. If it looks fine in qgis but is just shifted, then try a geotransform origin using the lower left corner coordinate. when in qgis, you can change your projection so that your mouse will read lat/lon and you can test the other corners to make sure they're good. you say you have the 4 corner coords...what is the lower left coordinate? – user1269942 Nov 03 '15 at 18:47
  • @user1269942 it is very strange because the top left corner of the image seems to be placed where the lat and lon are for the projection - ie the top left corner has lat = 45 and lon 0 and not what I set in the geotransform (-10, 45). The image shifts wehn I change the lat and lons in the CRS. I am confused as I thought that that the image was anchored at the lat lon defined by the geotranform. DOes this make sense to you? thanks – Scott Nov 03 '15 at 19:04
  • @user1269942 also if I try and plotting the image using matplotlib everything is fine so I cannot figure out where the discrepancies are – Scott Nov 03 '15 at 19:06
  • the geotransform sets the lower left. it should plot fine with matplotlib as it's just a numpy array. being spatially oriented is different. – user1269942 Nov 04 '15 at 01:05

0 Answers0