[gdal-dev] Projecting points on a mercator projection

Brice Lambi lambi at ucar.edu
Thu Apr 29 12:31:24 EDT 2010


This actually works just as expected but I was being dense.  The odd
points with the negative coords are in Alaska and Hawaii which is
outside of this domain.  After removing these and accounting for the
image size:

gx = gx % 1073
gy = gy % 689

everything works as expected.



lambi at ucar.edu wrote:
> 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
> 
> 
> 
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev



More information about the gdal-dev mailing list