[gdal-dev] Projecting points on a mercator projection

[lambi at ucar.edu]

Hello,

I'm trying to figure out how to take a list of points and put them on a
mercator projected image using the gdal python bindings.  I'm trying to
match it with a NDFD grid that I've been able to create and put on google
maps from a GRIB2 file.  This is what I have:

geo = [-14483681.357462916523218,
          7192.681255611132656,
          0,
          6915896.954714655876160,
          0,
          -6738.891395318105424]
nx = 1073
ny = 689
driver = gdal.GetDriverByName("GTiff")
dst_ds = driver.Create(fname, nx, ny, 4, gdal.GDT_Byte)
dst_ds.SetGeoTransform(geo)
oproj = osr.SpatialReference()
oproj.ImportFromWkt('PROJCS["unnamed",GEOGCS["WGS
84",DATUM["WGS_1984",SPHEROID["WGS
84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]')

iproj = osr.SpatialReference()
iproj.ImportFromWkt('GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]')

for i in range(1,5):
    dst_ds.GetRasterBand(i).WriteArray(numpy.zeros((ny, nx)))

transform = osr.CoordinateTransformation(iproj,oproj)
dst_ds.SetProjection(self.oproj.ExportToWkt())
Rband = numpy.zeros([dst_ds.RasterYSize,dst_ds.RasterXSize])
Aband = numpy.zeros([dst_ds.RasterYSize,dst_ds.RasterXSize])
geom = gdal.InvGeoTransform(geo)[1]

for point in points:
    lat = point[0]
    lon = point[1]
    (gx, gy, gz) = transform.TransformPoint(lon,lat)
    (gx, gy) = gdal.ApplyGeoTransform(geom, gx, gy)
    Rband[gy][gx] = 255
    Aband[gy][gx] = 255
dst_ds.GetRasterBand(1).WriteArray(Rband)
dst_ds.GetRasterBand(4).WriteArray(Aband)

This results in a semi-correct image but it also has lots of points that
are negative that wrap around the image and don't look right.  I also know
this isn't taking into account that the image size is 1073x689.  What I
expect from ApplyGeoTransform with the inverse geotransform is (x,y) where
(x,y) are all positive and within the 1073x689 domain or can be scaled to
fit in that domain but I don't understand the negative points.  Any
suggestions?

Thanks,
Brice