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. [email protected] 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 > [email protected] > http://lists.osgeo.org/mailman/listinfo/gdal-dev _______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
