Loading and reprojecting .nc file without projection

Question

I am trying to extract some bands from a .nc file. I don't know the file projection (if there is one), thus I created in GRASS GIS a newLocation with no projection (just a XY Cartesian plane) and tried to loaded the specific subset of the .nc I'm interested in. The file is states.nc, available here: https://luh.umd.edu/data.shtml specifically, LUH2 v2h Release 10/14/16. These are, for instance, the information of one subset of states.nc

gdalinfo NETCDF:"/myPath/states.nc":primf
Driver: netCDF/Network Common Data Format
Files: /myPath/states.nc
Size is 1440, 720
Origin = (-180.000000000000000,90.000000000000000)
Pixel Size = (0.250000000000000,-0.250000000000000)
Metadata:
  lat#axis=Y
  lat#long_name=latitude
  lat#standard_name=latitude
  lat#units=degrees_north
  lon#axis=X
  lon#long_name=longitude
  lon#standard_name=longitude
  lon#units=degrees_east

I am working at European level, thus at this point I switch to my location with EPSG:3035 projection and set my computational region (to a raster map of Europe, EU_3035) and resolution (1000 meters). I want to create a mask of the European raster I have and load (while reprojecting in EPSG:3035) each raster band I obtained from states.nc (using Europe as a mask).

for band in 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166;
do
echo "CREATE A MASK OF EU, CROP AND REPROJECT BAND $band TO LAEA3035"
r.mask --o raster=EU_3035
r.proj --o location=newLocation mapset=PERMANENT input='pastr.'$band output='pastr.'$band'_3035'
r.mask -r
done 

At this point, I get an error message:

WARNING: <PROJ_INFO> file not found for location <newLocation>
ERROR: Unable to get projection info of input map

I can see on my file browser some files I usually have in my other locations (e.g., PROJ_INFO, PROJ_UNITS) are missing. I cannot frame the exact problem and I am not sure if my workflow, so far, is correct.

I tried another approach with gdalwarp (that I used in the past to solve another projection problem). I first reproject the rasters I loaded from newLocation to WGS84.

gdalwarp -t_srs '+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0' -tr 0.01 0.01 /myPath/myRaster.tif /myPath/myRaster_WGS84.tif -overwrite 

I switch to my WGS84 location and load the files. Then if I check them with r.info and r.stats -c, I realize they are empty.

> r.info map=myRaster_WGS84
 +----------------------------------------------------------------------------+
 | Map:      myRaster_WGS84                 Date: Fri Aug 12 17:14:28 2022    |
 | Mapset:   PERMANENT                      Login of Creator: lisa            |
 | Location: WGS84                                                            |
 | DataBase: /myPath                                                          |
 | Title:                                                                     |
 | Timestamp: none                                                            |
 |----------------------------------------------------------------------------|
 |                                                                            |
 |   Type of Map:  raster               Number of Categories: 0               |
 |   Data Type:    FCELL                                                      |
 |   Rows:         100                                                        |
 |   Columns:      100                                                        |
 |   Total Cells:  10000                                                      |
 |        Projection: Latitude-Longitude                                      |
 |            N:         1N    S:          0   Res: 0:00:36                   |
 |            E:         1E    W:          0   Res: 0:00:36                   |
 |   Range of data:    min = NULL  max = NULL                                 |
 |                                                                            |
 |   Data Description:                                                        |
 |    generated by r.in.gdal                                                  |
 |                                                                            |
 |   Comments:                                                                |
 |    r.in.gdal --overwrite input="/myPath/myRaster_WGS84.tif" \              |
 |    output="myRaster_WGS84" memory=300 offset=0 num_digits=0                |
 |                                                                            |
 +----------------------------------------------------------------------------+

I guess it has to do with the XY location I created, but I don't know which is the appropriate projection and I cannot find this information. Searching on the internet I couldn't find a way to properly load and reproject bands of a .nc file I'm interested in. I am quite new to GRASS GIS and this is the first .nc file I manage. I am using GRASS 7.8.6 on a server running Ubuntu 18.04.6 from a Windows 10 machine.

Answer 1

Here's a way to work with that data. Following the hint in your question:

First with gdalinfo I found one of the variables, as you did. then I did the transform to WGS84 and clipped to the bounds of Europe along the way:

gdalwarp -t_srs EPSG:4326 -tr 0.5 0.5 -te 0 30 40 50 NETCDF:"states.nc":primf primf_WGS84.tif -overwrite

I chose a lower resolution (0.5 degrees) than what you want, to save some space. But the data itself is, I believe, only at 0.25 degrees anyway.

Then I started GRASS in a WGS84 defined location and did:

r.external primf_WGS84.tif output=primf

and I got 1000 Forest primary production maps.

Answer 2

Eventually, I tried two different approaches (example here with pastr subset and bands 1165 and 1166).

1. First approach

I assumed states.nc was in WGS84. I first loaded the pastr bands I was interested in into my WGS84 location (note the -o flag to override projection check):

for band in 1165 1166;
do
r.in.gdal --o -o input=NETCDF:"/myPath/states.nc":pastr band=$band
output='pastr_WGS84_'$band 
done

I switch to my EPSG:3035 location, set computational region (EU_3035) and resolution (1000 meters). I then reproject and cut the rasters to Europe:

for band in 1165 1166;
do
r.mask --o raster=EU_3035
r.proj --o location=WGS84 mapset=PERMANENT input='pastr_WGS84_'$band output='pastr_3035_'$band 
r.mask -r
done

I average (measures I'm interested in) the rasters:

r.series --o input=`g.list type=raster pattern=pastr_3035_* mapset=. separator=,` output=pastrAvg_3035 method=average

2. Second approach (based on Micha's answer)

In my WGS84 location, I reproject pastr from the unknown projection to WGS84 (but without directly clipping to EU):

gdalwarp -t_srs EPSG:4326 -tr 0.5 0.5 NETCDF:"/myPath/states.nc":pastr /myPath/pastr_WGS84.tif -overwrite

I load pastr_WGS84 into WGS84 location selecting only the bands I was interested in (with band):

for band in 1165 1166
do
r.in.gdal --o input=/myPath/pastr_WGS84.tif 
output='pastr_WGS84_'$band band=$band
done 

I switch to EPSG:3035 location, set computational region (EU_3035) and resolution (1000 meters). I import-reproject-clip the raster to Europe and average it, as in the first approach.

The two approaches led to slightly different pastrAvg_3035 rasters (if I look at the results of r.univar -e), in terms of non null cells and median values. I guess it depends from the approach but I cannot exactly understand why, and which one may be more correct (I guess the second - Micha's - approach).