Loading GeoTiffs#

[1]:

import starepandas
import rasterio
import numpy
import geopandas
import pyproj
import pystare

Usage#

[2]:

file_path = '../tests/data/wind.tif'

[3]:

sdf = starepandas.read_geotiff(file_path, add_pts=True, add_latlon=True,
                               add_coordinates=True, add_xy=True, add_trixels=True)

[4]:

sdf

[4]:

	band_1	lon	lat	ix	iy	x	y	geometry
0	6.272615	-120.650055	34.972664	4.0	1.0	-59314.902344	-338019.9375	POINT (-120.65006 34.97266)
1	6.252211	-120.647316	34.972679	5.0	1.0	-59064.902344	-338019.9375	POINT (-120.64732 34.97268)
2	6.258501	-120.650040	34.970409	4.0	2.0	-59314.902344	-338269.9375	POINT (-120.65004 34.97041)
3	6.258501	-120.647301	34.970425	5.0	2.0	-59064.902344	-338269.9375	POINT (-120.64730 34.97043)
4	6.232330	-120.644562	34.970440	6.0	2.0	-58814.902344	-338269.9375	POINT (-120.64456 34.97044)
...	...	...	...	...	...	...	...	...
88	5.839790	-120.622513	34.954800	14.0	9.0	-56814.902344	-340019.9375	POINT (-120.62252 34.95480)
89	5.825665	-120.619774	34.954815	15.0	9.0	-56564.902344	-340019.9375	POINT (-120.61978 34.95482)
90	5.817735	-120.617035	34.954830	16.0	9.0	-56314.902344	-340019.9375	POINT (-120.61704 34.95483)
91	5.817735	-120.614296	34.954845	17.0	9.0	-56064.902344	-340019.9375	POINT (-120.61430 34.95484)
92	5.815444	-120.611557	34.954861	18.0	9.0	-55814.902344	-340019.9375	POINT (-120.61156 34.95486)

93 rows × 8 columns

Detail:#

This is a geotiff containing average wind speeds in Santa Barbara county. It is projected in California Albers (EPSG:3310). Note that the raster is clipped to the county outline, thus containing NaNs

[5]:

with rasterio.open(file_path) as src:
    values = {}
    for band in range(1, src.count+1):
        values[f'band_{band}'] = src.read(band)
    height = values['band_1'].shape[0]
    width = values['band_1'].shape[1]
    cols, rows = numpy.meshgrid(numpy.arange(width), numpy.arange(height))
    xs, ys = rasterio.transform.xy(src.transform, rows, cols)
    xs = numpy.array(xs)
    ys = numpy.array(ys)
    src_crs = src.crs

[6]:

dst_crs = 'EPSG:4326'
transformer = pyproj.Transformer.from_crs(src_crs, dst_crs)
lats, lons = transformer.transform(xs, ys)
sids = pystare.from_latlon_2d(lats, lons, adapt_level=True)

[7]:

for key in values:
    values[key] = values[key].flatten()

[8]:

pts = geopandas.points_from_xy(lons.flatten(), lats.flatten())
gdf = geopandas.GeoDataFrame(values, geometry=pts, crs=dst_crs)
gdf = starepandas.STAREDataFrame(gdf, sids=sids.flatten())
gdf.dropna(inplace=True)
gdf.reset_index(inplace=True, drop=True)

[9]:

#gdf.to_file('wind.gpkg', driver='GPKG')

[10]:

gdf

[10]:

	band_1	geometry	sids
0	6.272615	POINT (-120.65006 34.97266)	3331746293595969199
1	6.252211	POINT (-120.64732 34.97268)	3331745340069319375
2	6.258501	POINT (-120.65004 34.97041)	3331746259443112015
3	6.258501	POINT (-120.64730 34.97043)	3331745271579949615
4	6.232330	POINT (-120.64456 34.97044)	3331745341877499695
...	...	...	...
88	5.839790	POINT (-120.62252 34.95480)	3331745323129218159
89	5.825665	POINT (-120.61978 34.95482)	3331745322842829231
90	5.817735	POINT (-120.61704 34.95483)	3331745305821847439
91	5.817735	POINT (-120.61430 34.95484)	3331745306553027855
92	5.815444	POINT (-120.61156 34.95486)	3331745573111838703

93 rows × 3 columns

[ ]: