last day (20 days later) » 

Feeds
Feeds
15:04
3
A: Python - How to write interpolated raster to netcdf file

ErikI changed your model solution a little bit; but it's like Ingvar Lukas wrote in his answer: you redefined xi and yi in the process, so when you later on define the netCDF values for lat and lon you try to fill small 1D arrays with a 2D array. That is the source of your error. import numpy as np f...

Creative crypter
Creative crypter
thanks for your answer Erik! I commented Ingvar´s answer as well, can you refer to it? Its probably just a transformation issue..
Erik
Erik
As you use scipy.interpolate.griddata, you get a linear interpolation of the data onto the grid; that's where the problem resides... not on the netCDF side. I don't have knowledge of whether netCDF and QGIS support unstructured (i.e., non-regular lat/long) grid data. but it sounds like that might be what you really want. Or you can, e.g., simply make a .CSV file with column1=lon, colum2=lat, column3=value, and load it as a set of points into QGIS? Wouldn't that be a better solution?
Creative crypter
Creative crypter
Hey, yea but i need to do it via a netcdf file generated in python.. is there no way to prepare the data in python so i can write it to netcdf? I thought i should first create a raster from the scattered data(latitude/longitude) and then write it to netcdf... When you look at the plot being generated, it actually looks correct..
Erik
Erik
Well, within the constraints of a raster file, then I'd suggest to set zi[:,:]=0 and loop over your data coordinates and fill zi only for those xi, yi coordinates that are minimal at the corresponding lon and lat coordinates. Assuming rectangular grids preserve Pythagorean distance, idx=np.argmin(np.sqrt( (xi-lon[i])**2 + (yi-lat[i])**2)) should give you the location of where to set zi[idx]=temp[i], where you should just loop over i...
I modified the answer to do what I described above.
Creative crypter
Creative crypter
Hi @Erik thank you so much, i start to understand how to work in this environment. I adjusted your source code with one little edit: I transformed the lat/lon coordinates into epsg:3857 format. The matplotlib output looks good, but when i place the file in QGIS there is no result. Can you check my updated source code?
I changed the spatial_ref string to:
crs.spatial_ref = """GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],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"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +un…
Erik
Erik
15:32
From
`from pyproj import CRS

CRS_3857 = CRS.from_epsg(3857)
CRS_wkt = CRS_3857.to_wkt(pretty=False)
print(CRS_wkt)`
I get
...something that is too big to print
run it on your own Python instance
So have this at the top:
`from pyproj import CRS

CRS_3857 = CRS.from_epsg(3857)
CRS_wkt = CRS_3857.to_wkt(pretty=False)`
And further down your code have
` crs.spatial_ref = CRS_wkt`
It shows up in QGIS at the correct location for me like that
I have modified the answer to show the other CRS
Creative crypter
Creative crypter
16:33
user image
Ah sorry i didnt research pyproj... ok i can follow your code! Do you know how i can display them as squares? Is this about the grid the data are projected in?
I am just wondering, when i visualize via matplotlib the coordinates appear as squares, but in QGIS more like vertical rectangles..

  last day (20 days later) »