Apologies in advance as I am not a GIS specialist by any means. I have a set of 1 million points and I'm trying to find their values in a geotiff raster. I've tried various versions of this answer involving affine transformations: https://gis.stackexchange.com/a/221471/143163
def retrieve_pixel_value(geo_coord, data_source): """Return floating-point value that corresponds to given point.""" x, y = geo_coord, geo_coord forward_transform = \ affine.Affine.from_gdal(*data_source.GetGeoTransform()) reverse_transform = ~forward_transform px, py = reverse_transform * (x, y) px, py = int(px + 0.5), int(py + 0.5) pixel_coord = px, py data_array = np.array(data_source.GetRasterBand(1).ReadAsArray()) return data_array[pixel_coord][pixel_coord]
This isn't ideal since it's point by point querying but it's better than nothing. The problem I'm having, however, is that it's not returning the correct values.
As a sanity check, I used some WRF data that had lat/long layers and queried various points and the resultant lat/long that is supposed to be the nearest coordinates in the WRF data are very far off from where they should be. e.g. inputting 33.77864,-117.33142 returns 38.72556, -115.75209 (the range of this layer is 32.597065 to 39.3944,-121.413025 to -113.04607, so there should be much closer matches). Furthermore, switching lat long as inputs doesn't drastically change the return value (when switched, returns 34.820377,-120.55661). Like I said, not an expert, but to me this seems like I'm using the wrong coordinate system as inputs. Does anyone know a way to convert lat long to the appropriate coordinates to find values in a raster?
I realize this isn't the most efficient way to do a list query on a big db, given the original problem of finding raster values for 1 million points, is there a more GIS-ish way to do this?