<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40"><head><meta http-equiv=Content-Type content="text/html; charset=utf-8"><meta name=Generator content="Microsoft Word 15 (filtered medium)"><style><!--
/* Font Definitions */
@font-face
{font-family:"Cambria Math";
panose-1:2 4 5 3 5 4 6 3 2 4;}
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
@font-face
{font-family:Consolas;
panose-1:2 11 6 9 2 2 4 3 2 4;}
@font-face
{font-family:"Trebuchet MS";
panose-1:2 11 6 3 2 2 2 2 2 4;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0in;
font-size:11.0pt;
font-family:"Calibri",sans-serif;}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:blue;
text-decoration:underline;}
pre
{mso-style-priority:99;
mso-style-link:"HTML Preformatted Char";
margin:0in;
margin-bottom:.0001pt;
font-size:10.0pt;
font-family:"Courier New";}
span.HTMLPreformattedChar
{mso-style-name:"HTML Preformatted Char";
mso-style-priority:99;
mso-style-link:"HTML Preformatted";
font-family:Consolas;
mso-ligatures:none;}
span.EmailStyle22
{mso-style-type:personal-reply;
font-family:"Calibri",sans-serif;
color:windowtext;}
.MsoChpDefault
{mso-style-type:export-only;
font-size:10.0pt;
mso-ligatures:none;}
@page WordSection1
{size:8.5in 11.0in;
margin:1.0in 1.0in 1.0in 1.0in;}
div.WordSection1
{page:WordSection1;}
--></style><!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]--></head><body lang=EN-US link=blue vlink=purple style='word-wrap:break-word'><div class=WordSection1><p class=MsoNormal>I didn’t realize that netCDF also has a vector component.<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal><a href="https://gdal.org/drivers/vector/netcdf.html#vector-netcdf">https://gdal.org/drivers/vector/netcdf.html#vector-netcdf</a><o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>To read the vector component, you’d use the ogr_fdw extension to read that rather than postgis_raster extension.<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>Details here - <a href="https://github.com/pramsey/pgsql-ogr-fdw">https://github.com/pramsey/pgsql-ogr-fdw</a><o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>So perhaps that is what your python is doing reading the vector component. I don’t know if netcdf mixes those in one data set or not since I have no experience using netcdf.<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>I recall you said you are using the OSGeo Live 16 distribution. I just checked the OSGeoLive 16 does include ogr_fdw <o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>So do <o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>CREATE EXTENSION ogr_fdw;<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>The list of supported formats you can see with this query:<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>SELECT name FROM unnest(ogr_fdw_drivers()) AS f(name) ORDER BY name;<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>Which for osgeolive 16, I see netCDF listed<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal><o:p> </o:p></p><div style='border:none;border-left:solid blue 1.5pt;padding:0in 0in 0in 4.0pt'><div><div style='border:none;border-top:solid #E1E1E1 1.0pt;padding:3.0pt 0in 0in 0in'><p class=MsoNormal><b>From:</b> Regina Obe <lr@pcorp.us> <br><b>Sent:</b> Wednesday, November 15, 2023 6:19 PM<br><b>To:</b> 'PostGIS Users Discussion' <postgis-users@lists.osgeo.org><br><b>Cc:</b> 'Manaswini Ganjam' <manu.ganjam@gmail.com><br><b>Subject:</b> RE: [postgis-users] Extracting variable information from netcdf, imported as raster to a table<o:p></o:p></p></div></div><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>Just confirming some stuff, since I’m not completely following:<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>Raster_record.rast column is of type raster correct? IF so ST_X and ST_Y won’t work since those are for geometry types.<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>Also ST_Value(raster_record.rast, band_number), won’t work either since that expects as input a geometry or x,y on the raster you want the value you.<o:p></o:p></p><p class=MsoNormal>I would think you would have gotten an error with that, which makes me feel I’m missing something critical.<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>If you want to extract all the pixels in a raster, you’d do something like <a href="https://postgis.net/docs/RT_ST_PixelAsPoints.html">https://postgis.net/docs/RT_ST_PixelAsPoints.html</a><o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>SELECT pp.x, pp.y, pp.val, ST_X(pp.geom) AS lon, ST_Y(pp.geom) AS lat<o:p></o:p></p><p class=MsoNormal>FROM raster_record, <o:p></o:p></p><p class=MsoNormal>ST_PixelAsPoints(raster_record.rast, 1) AS pp<o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal><br> <o:p></o:p></p><p class=MsoNormal><o:p> </o:p></p><div style='border:none;border-left:solid blue 1.5pt;padding:0in 0in 0in 4.0pt'><div><div style='border:none;border-top:solid #E1E1E1 1.0pt;padding:3.0pt 0in 0in 0in'><p class=MsoNormal><b>From:</b> postgis-users <<a href="mailto:postgis-users-bounces@lists.osgeo.org">postgis-users-bounces@lists.osgeo.org</a>> <b>On Behalf Of </b>Manaswini Ganjam via postgis-users<br><b>Sent:</b> Wednesday, November 15, 2023 2:01 PM<br><b>To:</b> <a href="mailto:postgis-users@lists.osgeo.org">postgis-users@lists.osgeo.org</a><br><b>Cc:</b> Manaswini Ganjam <<a href="mailto:manu.ganjam@gmail.com">manu.ganjam@gmail.com</a>><br><b>Subject:</b> [postgis-users] Extracting variable information from netcdf, imported as raster to a table<o:p></o:p></p></div></div><p class=MsoNormal><o:p> </o:p></p><div><p class=MsoNormal>Hi, <o:p></o:p></p><div><div><p class=MsoNormal>I have been trying to download s3 cloud stored gridded climate data and generate tables with variables, lat, lon and timestamp (year, yearday). To achieve this I used raster2pgsql and imported multiple netcdf files into a database table. <o:p></o:p></p><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>Question: How to achieve the extraction of variables using postgis? I tried using ST_value, ST_pixelaspoints but I was getting errors, mainly due to the format in which netcdfs are stored in the database (the error says can't load some characters like 00E30100082...), I even tried changing the datatype to float but still did not work. I mean it is probably not simple like selecting a variable from the netcdf. I have enclosed my sql query below: <o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal> -- Iterate through all raster files in the table<br> FOR raster_record IN (SELECT * FROM gfdl_03_prcp) LOOP<br> -- Determine the year from the raster file name, assuming the format is '<a href="http://prcp_03_1950.nc">prcp_03_1950.nc</a>'<br> SELECT<br> regexp_replace(raster_record.filename, '.*_(\d{4})\.nc', '\1')::integer<br> INTO<br> year;<br> <br> -- Calculate the start date of the year<br> year_start := (year || '-01-01')::date;<br> <br> -- Determine if the year is a leap year<br> is_leap_year := EXTRACT(ISODOW FROM (year_start + interval '1 year')) = 7;<br> <br> -- Set the number of bands for the year (365 for non-leap years, 366 for leap years)<br> FOR band_number IN 1..(CASE WHEN is_leap_year THEN 366 ELSE 365 END) LOOP<br> -- Calculate the observation_time using the year and band number<br> observation_time := year_start + (band_number - 1) * interval '1 day';<br> <br> -- Extract X (lon) and Y (lat) coordinates from the raster<br> SELECT<br> ST_X(raster_record.rast) AS lon,<br> ST_Y(raster_record.rast) AS lat<br> INTO<br> lon,<br> lat;<br> <br> -- Insert the lat, lon, prcp, and observation_time into the extracted_values table<br> INSERT INTO extracted_values (lat, lon, prcp, observation_time)<br> VALUES<br> (lat, lon, ST_Value(raster_record.rast, band_number), observation_time);<br> <br> -- Increment the counter<br> counter := counter + 1;<br> <br> -- Commit the transaction periodically in batches<br> IF counter % batch_size = 0 THEN<br> COMMIT;<br> END IF;<br> END LOOP;<br> END LOOP;<br> <o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>The metadata for the two files is as follows: <o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>File from database:<o:p></o:p></p></div><div><div><div><div><div><div><pre>{'NC_GLOBAL#Conventions': 'CF-1.5',<o:p></o:p></pre><pre> 'NC_GLOBAL#GDAL': 'GDAL 3.6.4, released 2023/04/17',<o:p></o:p></pre><pre> 'NC_GLOBAL#history': 'Wed Nov 15 13:32:13 2023: GDAL CreateCopy( <a href="http://not_clipped_prcp.nc">not_clipped_prcp.nc</a>, ... )'}<o:p></o:p></pre></div></div></div></div></div><div><pre>File before loading into the database:<o:p></o:p></pre><pre>{'lat#units': 'degrees_north',<o:p></o:p></pre><pre> 'lon#units': 'degrees_east',<o:p></o:p></pre><pre> 'NC_GLOBAL#title': 'Daily statistically downscaled CMIP5 data for the United States and southern Canada east of the Rocky Mountains, version 1.0, realization 1, 0.1x0.1 degree spatial resolution.',<o:p></o:p></pre><pre> 'NETCDF_DIM_EXTRA': '{time}',<o:p></o:p></pre><pre> 'NETCDF_DIM_time_DEF': '{366,4}',<o:p></o:p></pre><pre> 'NETCDF_DIM_time_VALUES': '{0,1,2,3,4,5,6,7,8,9,10,11,12,13,14....362,363,364,365}',<o:p></o:p></pre><pre> 'prcp#add_offset': '819.17499',<o:p></o:p></pre><pre> 'prcp#long_name': 'daily precipitation accumulation',<o:p></o:p></pre><pre> 'prcp#missing_value': '-32768',<o:p></o:p></pre><pre> 'prcp#scale_factor': '0.025',<o:p></o:p></pre><pre> 'prcp#units': 'mm',<o:p></o:p></pre><pre> 'prcp#_FillValue': '-32768',<o:p></o:p></pre><pre> 'time#units': 'days since 1952-1-1 0:0:0.0'}<o:p></o:p></pre><pre><o:p> </o:p></pre><pre><span style='font-family:"Arial",sans-serif'>In case this information is useful: Previously I used python to extract variable information and generate a csv or table using this variable information, and the code is enclosed below. In the code I extracted variable values using lon = dataset.variables['lon'][:] and iterated for loops to write them all in csv. </span><o:p></o:p></pre></div><div><div><p class=MsoNormal>Python code:<o:p></o:p></p></div><div><p class=MsoNormal>import netCDF4 as nc<br><br># Step 1: Read the NetCDF file<br>filename = "/home/manaswini/<a href="http://prcp_03_1950.nc">prcp_03_1950.nc</a>"<br>dataset = nc.Dataset(filename)<br>dataset.set_auto_mask(False)<br>dataset.set_auto_scale(True)<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>lon = dataset.variables['lon'][:]<br>lat = dataset.variables['lat'][:]<br>time = dataset.variables['time'][:]<br>prcp = dataset.variables['prcp'][:]<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>import numpy as np<br>import csv<br><br># csv_buffer<br>csv_buffer = open('output.csv', 'w', newline='')<br>csv_writer = csv.writer(csv_buffer)<br><br># Iterate through grid points and write to CSV buffer<br>for i in enumerate(lon):<br> for j in enumerate(lat):<br> for k in enumerate(time):<br> csv_writer.writerow([lat[j], lon[i], prcp[i][j][k], year[k], yearday[k]])<br><br><br># Close the CSV buffer<br>csv_buffer.close()<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>Thank you, <o:p></o:p></p></div><div><p class=MsoNormal><span style='font-family:"Trebuchet MS",sans-serif'>Manaswini Ganjam</span><o:p></o:p></p></div></div></div></div></div></div></div></div></div></body></html>