texts: Importing Shapefile to PostGIS using QGIS DB Manager

texts:
I am trying to import a shapefile from my QGIS project to my PostGIS schema using DB manager.

I want to create a new table in my public schema and I have an ERROR: relation "spatial_ref_sys" does not exist.

Checking options 'primary key', 'geom column' & 'create spatial index' do not alter the error. Full error message below.

Creation of data source "public"."fwamanc" failed:

Erroneous query: SELECT AddGeometryColumn('public','fwamanc','geom',27799,'MULTIPOLYGON',2) returned 7

[ERROR:  relation "spatial_ref_sys" does not exist

LINE 1: SELECT SRID         FROM spatial_ref_sys WHERE SRID = new_sr...

^

QUERY:  SELECT SRID         FROM spatial_ref_sys WHERE SRID = new_srid

CONTEXT:  PL/pgSQL function addgeometrycolumn(character varying,character varying,character varying,character varying,integer,character varying,integer,boolean) line 48 at SQL statement SQL statement "SELECT public.AddGeometryColumn('',$1,$2,$3,$4,$5,$6,$7)" PL/pgSQL function addgeometrycolumn(character varying,character varying,character varying,integer,character varying,integer,boolean) line 5 at SQL statement ]

I have already created the postgis extension & checked validity of output.
This is the first time I have tried to import without restoring via pgdump files.

Anybody help with what I may be missing?


texts: Is it possible update a shapefile instead create a new one when using Processing?

texts: Using Processing to add or calculate a field it's always necessary create a new shapefile? Is there a way to just update and keep using same layer?


texts:
I'm using ArcMap 10.2, I'd like to restrict the type of records that makes their way into a feature class to only those with a geometry. Is there a way to create feature class in a geodatabase where the SHAPE field is set to non-nullable?

Ie. "Allow NULL Values" changed to "No"?


texts: Exporting GEE satellite imagery to Numpy array

texts:
I'm trying to extract LANDSAT-8 images around a lat/lon location, using the Google Earth Engine API for Python.
Specifically, I want to generate a bounding box of 1km by 1km around my point, and extract the pixels for the RGB bands, so my output should be 33 by 33 by 3, given that LANDSAT-8 resolution is 30m.
However, I can't get this to work. I'm generating my box using ee.Geometry.Point.buffer and ee.Geometry.Polygon.bounds, and then I'm trying to overlay my geometry over LANDSAT-8 images of 2020, of which I want to get the average value pixel-wise (so the annual average image should be 33 by 33 by 3). After that, I plan to extract the data as a numpy.Array using the ee.Image.sampleRectangle function, but it is not working. Here is a minimal working example:
#importing packages
import numpy as np
import ee
ee.Authenticate()
ee.Initialize()

#converting a lon/lat pair to an ee.point object
point=ee.Geometry.Point( [-122.2036486, 37.4237011] )
# converting the point to a patch: we define a circle with a 500m radius, and then we put a box around it
patch=point.buffer(500).bounds()

#checking that the patch has approximately the right area 1000^2 m^2
patch.area(1).getInfo()

#defining the image: Landsat 8 collection of images from 2020, with RGB bands. We take the mean of them
imagery=ee.ImageCollection("LANDSAT/LC08/C01/T1").filterDate('2020-01-01', '2020-12-31').select(['B3','B4','B2']).mean()

#getting a matrix from the image for our patch of interest
rect_image = imagery.sampleRectangle(region=patch)

#we should obtain a 33*33 matrix for each one of the 3 bands

#extracting one band to check
band_b4 = rect_image.get('B4')
#checking the shape of the resulting matrix
np.array(band_b4.getInfo()).shape

#the result is 1*1

I feel that one problem here is that the mean() function is reducing collapsing the dimensions of the images, so my "fix" to that is instead using the first() function to get a proper image, even if it is not the average image I want. Here is what happens:
imagery2=ee.ImageCollection("LANDSAT/LC08/C01/T1").filterDate('2020-01-01', '2020-12-31').select(['B3','B4','B2']).first()

#checking the dimensions
imagery2.getInfo()['bands'][1]['dimensions']

#getting a matrix from the image for our patch of interest
rect_image2 = imagery2.sampleRectangle(region=patch)

#in theory, since this is 1000m*1000m and the resolution of the satellite is 30m,
#we should obtain a 33*33 matrix for each one of the 3 bands

#extracting 1 band to check
band_b4_2 = rect_image2.get('B4')

#checking the shape of the resulting matrix
np.array(band_b4_2.getInfo()).shape

#now we get an error

How could I solve this?
I have tried a lot of things like using ee.ImageCollection.filterBound over the image collection, ee.Image.clip, and I have tried to export the images in TIFF or TFRecord formats with the ee.batch.Export.image.toDrive method, but that also hasn't worked.


texts: Edit map legend properties from pause drawing screen ArcMap

texts: Is there a way to edit a legend's properties from the Paused Drawing screen in ArcMap? The reason I ask is so I can update a legend's properties without having it redraw each time.


texts: Calculate ground nadir line from PV coordinates

texts:
I'm using an API which returns a satellite's PVCoordinates (Position/Velocity pairs).

{P(1036869.5330731269, 6546536.584960848, 9.892293545390519E-11), V(-9994.355195123433, 1582.9503539842497, -1242.4491245907245)}

Given an array of these PV coordinates, how can I calculate or derive latitude and longitude coordinates for the nadir line (ground sub-track)?

Does anyone know of a Java API which might help me or is the math easy enough that I can write a function myself? (Full Disclosure: I have a math disorder)

I've looked around the web and everything I see tends to make my eyes glaze over and humble me to the point of wishing my pathetic self should have taken more math in college.


texts: How to define precisely the radius size when using 'select features by radius' in QGIS

texts:
I'd like to select features with circles having precise radius in QGIS.
By default this tool is freehand so it is hard to get a precise value.


texts: ArcGIS JS API: How to get the lat/long of the viewport vertices?

texts: I would like to get the min lat, min long,  max lat, max long of my viewport in ArcGIS JS API.  Is this possible?


texts: Populate field based on value within a field in another table

texts:
I have two datasets: dat_1 and dat_2.
I would like to populate field_1 in dat_1 based on values within field_1 in dat_2. I don't know if i need it, but they have ogr_ID field in common, with the same values in both tables.

i'm trying with this script but the the tool is populating field_1 in dat_1 with the same value for all the rows, and it's also printing always the same value.

fc = r"path/dat_1"
tbl = r"path/dat_2"
field = "field_1"

s_cursor = arcpy.da.SearchCursor(tbl, field)
for row in s_cursor:
    load = row[0]

with arcpy.da.UpdateCursor(fc, field) as u_cursor:
    for row in u_cursor:
    row[0] = load
    u_cursor.updateRow(row)
    print (row)

del row
del u_cursor


texts: Unable to modify form in Access when linked to ArcGIS using database connection?

texts:
I created a split Access 2010 database with a back-end for just tables and the front end with queries and forms.  The front-end is linked using an OLE DB connection to ArcGIS, which works fine.  When the database connection is connected I am unable to modify my forms but I can update and create queries.

Could this be an issue with having the forms bound to my table query? or is it because I am using a table query for the form?  I don't know how to update my table records from the form without it being bound to the table. I am also using the Microsoft Office 12.0 Access Database Engine OLE DB Provider.

This is a multiple user database and other users maybe linking this to ArcGIS and might also need to make modifications.


texts: Georeferencing an old map in QGIS using Thin Plate Spline

texts:
![old map georeferenced badly][1]

[1]:

I georeferenced this old map using Thin Plate Spline, and Cubic Spline as transformation settings. However, the result was quite bad. Do I have to add points througout the whole rasterfile for this to be even remotely readable, or is the map just too bad to georeference in a good way?


texts: Is it possible to rotate annotation labels in Print Composer?

texts: Is it possible to rotate a labels (not feature labels but manual annotations) in Print Composer in QGIS 1.8?


texts: Performance regression: fiona.transform much slower since GDAL 3

texts:
I have run into some severe performance degradation when upgrading an environment to GDAL3. I could track the issue down to fiona.transform, which is a lot slower (about 15 (!) times) now than it was with GDAL 2.4.
The issue can be illustrated using this line, which only transform one point (the actual script transform a geometry):
python -m timeit -s "from fiona.transform import transform" "transform('EPSG:31287', 'EPSG:4236', [419908], [333400])"

These are my performance measurements with the images from perrygeo/gdal-base and using fiona 1.8.13:
latest  Python 3.8.5 | GDAL 3.1.3 | GEOS 3.8.1 | PROJ 7.1.1 | 20 loops, best of 5: 20 msec per loop
20181219-6f5f6a29   Python 3.6.7 | GDAL 2.4.0 | GEOS 3.7.1 | PROJ 5.2.0 | 1000 loops, best of 3: 675 usec per loop
20181219-f379ec62   Python 3.6.7 | GDAL 2.4.0 | GEOS 3.7.1 | PROJ 5.2.0 | 1000 loops, best of 3: 705 usec per loop
20181221-40f73e30   Python 3.6.7 | GDAL 2.4.0 | GEOS 3.7.1 | PROJ 5.2.0 | 1000 loops, best of 3: 698 usec per loop
20181221-bc2d4bbd   Python 3.6.7 | GDAL 2.4.0 | GEOS 3.7.1 | PROJ 5.2.0 | 1000 loops, best of 3: 688 usec per loop
20181221-f7a0a299   Python 3.6.7 | GDAL 2.4.0 | GEOS 3.7.1 | PROJ 5.2.0 | 1000 loops, best of 3: 703 usec per loop
20190312-f69f8699   Python 3.6.8 | GDAL 2.4.0 | GEOS 3.7.1 | PROJ 6.0.0 | 1000 loops, best of 3: 634 usec per loop
20190322-800eed8a   Python 3.6.8 | GDAL 2.4.1 | GEOS 3.7.1 | PROJ 6.0.0 | 1000 loops, best of 3: 589 usec per loop
20190509-da2e635a   Python 3.6.8 | GDAL 3.0.0 | GEOS 3.7.2 | PROJ 6.0.0 | 100 loops, best of 3: 10.4 msec per loop
20191110-6cc84c7e   Python 3.6.9 | GDAL 3.0.2 | GEOS 3.8.0 | PROJ 6.2.1 | 100 loops, best of 3: 11.4 msec per loop
20200301-8437abbb   Python 3.8.2 | GDAL 3.0.4 | GEOS 3.8.0 | PROJ 7.0.0 | 20 loops, best of 5: 10.5 msec per loop
20200509-50546ca8   Python 3.8.2 | GDAL 3.1.0 | GEOS 3.8.1 | PROJ 7.0.1 | 20 loops, best of 5: 10.2 msec per loop
20200907-c7ec91bc   Python 3.8.5 | GDAL 3.1.3 | GEOS 3.8.1 | PROJ 7.1.1 | 20 loops, best of 5: 10.8 msec per loop

Once can clearly see that the line performs at ~0.7 msec before GDAL3 and beginning with GDAL3 the line takes >10 msec to finish.
Does anyone have a hint, what could be the root of the issue and how it could be fixed?


texts: User-defined function for catalog_apply does not return las with newly computed additional las attributes

texts:
I am trying to compute the actual incidence angles.
The following is the code that I used based on the example in the lidR book
Code:
ang_calc = function(cluster, traj)
{
  las = readLAS(cluster)
  if (is.empty(las))return(NULL)
  las <- lastraj(las, traj, use.fun = TRUE)
  print(names(las@data))
  las <- lasangles(las, traj = c("xtr", "ytr", "ztr") )
  print(names(las@data))
  return(las)
}

output <- catalog_apply(lasc, ang_calc, traj)

The print commands confirm that the additional attributes are being added after each function. However, the final point cloud that is written on disk has none of the additional attributes. There was the following warning:
1: In regularize.values(x, y, ties, missing(ties), na.rm = na.rm) :
  collapsing to unique 'x' values

Just to test, I only read the cluster in the function and returned it without doing any thing. There was no warning and there was no difference in the files written on disk in both the cases.


texts: What's an efficient way of detecting road junctions on a map?

texts:
I may receive maps from any source for my project, and for the C++ algorithm I'm making, it's imperative for the algorithm to know where road junctions (nodes where more than two road segments meet) are. Besides, in many of the maps I get, there are roads overshooting and undershooting (ie: roads aren't properly connected). One way to correct the problem is by 'cleaning' the road network using Autocad.
The questions I have are:

Is there an efficient way (any
software which can do it?) to get my
map in a GML format, which contains
information about road junctions in
the map? (maybe have a GML tag which
says that a particular node is a
junction)
Is there any other way to 'clean'
the road network?

Colleagues have suggested FME, but that involves writing scripts and we're not sure about whether the script would be flexible enough to cater to all maps. The only other way to detect junctions is to use brute force and find which road segments have common nodes. Would ArcGIS help? (haven't used it, but have heard of it) I'm sure there should be a better way...


texts: I have X,Y data that I am trying to buffer to find clusters. My definition for a cluster is " 2 or more points with a maximum distance of 200 meters".  My problem is:  Point A, Point B, and Point C included into a a single cluster.  Point A to Point B is 50m, Point A to C is 190m, but point B to C is 225 meters.  Therefore I do not believe that all points in my cluster are meeting my definition of a cluster?


texts: Common custom toolbar for several different arcgis desktop addins using C#

texts:
Suppose, I have two  arcgis desktop add ins C1Addin and C2Addin developed in c# with the following components:

C1Addin:
 Button_1

C2Addin:
  Button_2

When the 1st addin which may be any of the two addins is installed, it will create a toolbar and add the UI components of the corresponding addin to the toolbar. When the other addin is installed, it will add its components to that toolbar.

I have several problems to implement the above behavior:

How can I override the default initialization behavior of addin so that I can check whether the toolbar exists and if not I can dynamically create the toolbar?
How can a toolbar is uniquely identified by two or more addins?


texts: Is it possible to run a cartesian join in QGIS?  I use Access to do this but I would like to know if I can do it directly in QGIS (or GRASS GIS) to avoid the back and forth between the GIS and Access.


texts: Understanding how horizontal factor in ArcGIS Path Distance tool works?

texts:
I work with ArcGIS 10.1 (advanced license), and I am finding difficult to grasp how the horizontal factor works in the Path Distance tool. I have read the help documentation several time. These are the sources:
(1) http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/understanding-path-distance-analysis.htm
(2) http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/how-the-horizonal-and-vertical-factors-affect-path-distance.htm
(3) http://desktop.arcgis.com/de/arcmap/10.3/analyze/arcpy-spatial-analyst/an-overview-of-horizontal-factor-classes.htm

Now, sticking with the ArcGIS help documentation's example of a car moving from location A to location B, I do understand that the use of a horizontal factor in the context of path distance analysis is actually aimed at increasing a given movement's cost on the basis of an additional horizontal friction. In other words, to stick with the car example, if we only consider the slope as a cost and the car is moving uphill, there will be a given cost (say 10); if we add an additional friction caused by the wind blowing at 180 degree offset (i.e., in a direction opposite the car movement's direction), the total cost will be increased by a factor f (in total: fx10).

So far so good.
YET, I am rather confused here:
To determine the HF for moving from one cell to the next, the prevailing horizontal direction at the processing cell must be established from the horizontal direction raster. The horizontal relative moving angle (HRMA), the relationship between the moving direction with respect to the horizontal direction, is calculated and plotted on the horizontal factor graph to determine the horizontal factor. The resultant value establishes the horizontal factor from the center of the processing cell (the FROM cell) to the edge of the TO cell (source: http://desktop.arcgis.com/de/arcmap/10.3/analyze/arcpy-spatial-analyst/an-overview-of-horizontal-factor-classes.htm).

That paragraph is totally obscure to me. In particular, I do not understand what is the horizontal direction raster to which reference is made. Is that the cost raster which is one of the input required by the tool?

In summary: let's assume I have a raster representing the cost of the car movement in a given area, and I want to use the Path Distance tool to calculate the accumulated cost from a origin location also accounting for the additional horizontal friction caused by the wind.

Questions:

(A) Am I correct in understanding that I need to only input the main cost raster AND the horizontal factor data?

(B) is the main cost raster equivalent to the aforementioned (rather obscure) horizontal direction raster?

(C) As for the horizontal factor data, assuming I want to use a customized table in which I assign a factor value to each wind direction, am I correct in understanding that I have to use 0 degrees to "represent" wind blowing from behind (i.e., blowing in the car movement's direction)? I ask this because this point is not clear to me from the help documentation:To define the HRMA, the angle between the horizontal direction defined on the horizontal direction raster and the moving direction to the cell to which the cost distance is being calculated is determined. Even though the values range from -180 to +180, the HRMAs are symmetrical with respect to the horizontal factor axis; therefore, they are represented from 0 to 180 on the x-axis of the horizontal factor graph. The HRMA is the number of degrees the travel direction is from the specified horizontal direction, regardless of positive or negative orientation. Therefore, the range of HRMAs is 0 to 180 degrees. (source: http://desktop.arcgis.com/de/arcmap/10.3/analyze/arcpy-spatial-analyst/an-overview-of-horizontal-factor-classes.htm).


texts:
I am looking to plot the ideal placement of police vehicles within a specific precinct (polygon). Rather than placing them myself, calculating a buffer of travel time (e.g. within 3 minutes), and using trial and error I'd like for the best placement to be determined pragmatically. By the ideal placement I mean the positions of the points (cars) that covers the highest percentage of the roads within the polygon.

Is there a way to accomplish this using Python/model builder?


texts: Can 2 GIS layers which representing the same area be compared directly while having the same datum but different coordinate system?

texts:
My question is about how coordinate systems affect comparison of different layers - I have tried to refresh my memory about coordinate systems in GIS, but am still unsure about this situation.

I am comparing the land cover classifications (estimates of the % each land cover type represents in the given area) from a raster layer versus from a vector shape used in photo-interpretation. This comparison will be done on many raster-classification/photo-interpretation combinations throughout the United States.

The area of the vector layer comes from delineating a watershed (generally 25 to 250 square kilometers in area), and that delineation and vector layer is then used to clip that area out of a larger land classification raster layer.

The vector layers used in photo-interpretation are all in lat/long (unprojected). The photo-interpretation is already done, so I need to keep these source shapes as they are. These vector layers have the following coordinate system:

Geographic Coordinate System:   GCS_North_American_1983
Datum: D_North_American_1983

All of the raster layers I'm working with have the following coordinate system:

Spatial Reference:  Albers_Conical_Equal_Area
Datum: D_North_American_1983

Are the areas of these 2 GIS layers the same, given they have the same Geographic Coordinate System, or am I in error to compare the vector layer to the raster layer, given that the vector layer is not projected and the raster layer is?
If the latter, then can I simply project the vector layer to the same coordinate system as the vector layer, and if so can that be done after clipping the raster to the vector shape or must the projection be done before clipping?


texts:
I was working on setting development workspace for GvSIG. I referred to https://confluence.prodevelop.es/display/WSS/GvSIG+Installation.

While compiling a simple program it says "extension directory is empty".
In that tutorial it does not explain how to add extensions in extension directory.

Is their any good tutorial on how to setup development environment for GvSIG
or does anyone have experience with this?


texts: PyQGIS Script failing to project layer to EPSG:54002

texts:
I'm trying to write a Python script to convert a vector layer using 4326 to EPSG:54002 World Equidistant Cylindrical for road distance calculations, but the script fails to do the conversion, I've tried numerous methods to get it right:

Using the re-project layer algorithm:

reprojectedLayer = processing.run("native:reprojectlayer", {'INPUT': roadLayer ,'TARGET_CRS': 'EPSG:54002', 'OUTPUT':'TEMPORARY_OUTPUT'})['OUTPUT']

reprojectedLayer = processing.run("native:reprojectlayer", {'INPUT': roadLayer ,'TARGET_CRS': QgsCoordinateReferenceSystem('EPSG:54002'), 'OUTPUT':'TEMPORARY_OUTPUT'})['OUTPUT']

Calling the setCrs method on the layer:

roadLayer.setCrs(QgsCoordinateReferenceSystem(54002, QgsCoordinateReferenceSystem.PostgisCrsId))

roadLayer.setCrs(QgsCoordinateReferenceSystem('EPSG:54002'))

Even just writing the layer to a file while specifying a different CRS:

QgsVectorFileWriter.writeAsVectorFormat(roadLayer, sys.argv[2] + "/roads", "utf-8",  QgsCoordinateReferenceSystem('EPSG:54002'), "GPKG")

Unfortunately, none of these have worked. If I save the file after any of the operations above and open it in QGIS it either just opens it fine as 4326 or opens it and displays a message that the CRS is invalid and that it will be defaulting to 4326.

One interesting thing to note is that when I specify the crs while writing the file, this pops up in the console:

proj_create_from_database: crs not found

It doesn't do that with the other methods, only when I write the file. This message doesn't make sense since if I try to reproject the layer to 54002 in QGIS, it works fine. Its only in the script that it fails.


texts: Redistribute raster pixel value in underlying polygons

texts:
Given a raster layer (e.g 10km resolution) and a vector layer with many features (with finer resolution), is it possible to redistribute the value of each raster pixel in an attribute among the underlying polygons, so that the sum of that attribute in the underlying features is equal to the overlaying pixel value?

If so, can this be done either homogeneously across underlying features, or based on a weight given by another attribute in the polygon layer?

QGIS-based solutions are favoured.


texts: Trying to insert geometry in existing table (geom field stay empty) using PostGIS?

texts:
I'm trying to insert a new line in a big table with the following SQL using PostgreSQL 9.3 and PostGis 2.1.3:

INSERT INTO "Edigeo"."GEOM"(
            cnig, geom_surf,
           etat, codcom)
    VALUES  ('H_11_4_6',
        ST_Difference((SELECT "GEOM_1".geom_surf FROM "Edigeo"."GEOM" AS "GEOM_1" WHERE "GEOM_1".codcom = '73101' AND "GEOM_1".cnig = 'H_1_6_0' LIMIT 1),(SELECT ST_Union("GEOM_2".geom_surf) FROM "Edigeo"."GEOM" AS "GEOM_2" WHERE "GEOM_2".codcom = '73101' AND "GEOM_2".cnig IN ('H_11_4_0','D_1_0_8')  AND ST_IsValid("GEOM_2".geom_surf) = True)),
        NULL,
        '73101')

my text field  are filled without problem but i have an issue to fill the geometry column (geom_surf) with the ST_Difference statement... It remains empty in the table after the execution.

If i do the same in a select statement, the geometry seems ok.

SELECT ST_Difference((SELECT "GEOM_1".geom_surf FROM "Edigeo"."GEOM" AS "GEOM_1" WHERE "GEOM_1".codcom = '73101' AND "GEOM_1".cnig = 'H_1_6_0' LIMIT 1),(SELECT ST_Union("GEOM_2".geom_surf) FROM "Edigeo"."GEOM" AS "GEOM_2" WHERE "GEOM_2".codcom = '73101' AND "GEOM_2".cnig IN ('H_11_4_0','D_1_0_8')  AND ST_IsValid("GEOM_2".geom_surf) = True))

returns

So it looks nice.

But in my table after the insert:

Is there an issue with my geometry? or something else? I don't understand why it's not working...


texts: Select only features within the layout view

texts:
I want to create a tool layer in ArcGIS10 that will select only features within the current layout view so that I can add the attribute table to the layout and only show attributes for the selected features.

Has anyone got any ideas on how to do this?  I want the selection to be done automatically each time the mxd is opened, and not have to do it manually.


texts:
I have few hazard maps (with color coding) of a particular state in a pdf report. I am a newbie in GIS and need to use these different maps from the pdf for combined interpretation. Can the colors representing different range of values in each map be imported in ArcGIS?

I have converted one of the maps into an image:


texts:
After this QGIS doesn't start

System config
When double click on QGIS Desktop 3.14.1 is shows the screen as attached but after that the program does not start or run. I have also tried "run as administrator" but still the problem persists. I also run the software with Internet connection ON but in vain. Once I was able to recover the crash file for v3.10 (attached herewith). With v3.10 and v3.14.1, same problem is occurring. I have tried to search the internet for various blogs, sites listed in your website support page, still no solution is there for my particular problem. I have installed software again with "Run as Administrator" but the issue is not solved.
.NET 4.8 already installed.
C++2005 to 2008 all are installed.

Error for 3.10 Crash ID: 5c75d2351852793e7098138ee0f44da41b6fa5f4
Stack Trace
CleanupThreadTlsData : DllMain :
_DllMainCRTStartup : LdrpCallInitRoutine : LdrShutdownThread : RtlExitUserThread : BaseThreadInitThunk : RtlUserThreadStart :
QGIS Info QGIS Version: 3.10.7-A Coru�a QGIS code revision: 7b4ca4c8d0
Compiled against Qt: 5.11.2 Running against Qt: 5.11.2 Compiled
against GDAL: 3.0.4 Running against GDAL: 3.0.4
System Info CPU Type: x86_64 Kernel Type: winnt Kernel Version:
10.0.18363

What to do?


texts: How can i create a link to a whole folder?

texts:
I know how to create a hyperlink to a single pdf-file: Create a hyperlink to a pdf file in the attribute table

But now I want to make a link to a folder with 3 or 4 pdf files in it. Is that possible?
I tried it the same way as linking a pdf but it didn't work.

I'm working ith qgis 2.0.1.


texts: Applying single-band color map from file using pyQGIS?

texts:
I am importing a number of rasters using the python console in QGIS 2.8. I would like to color them all using the same color map, which I have saved to a file. I would like to do this using python so that I don't have to navigate the gui for each file (i.e., right click on the layer -> properties -> style; Render type=Singleband psuedocolor; load color map from file).

My current approach is to follow the Using Raster Layers tutorial. However, I just read the contents of my color map file and convert it to the input of setColorRampItemList().

Here's my code so far

from qgis.core import QgsRasterLayer
from qgis.core import QgsColorRampShader
from PyQt4.QtCore import QFileInfo
from PyQt4.QtGui import QColor

def GetColorRampItemListFromText(file_contents):
    file_contents = file_contents.split('\n')[2:] # get rid of header
    file_contents = [ map(int,line.split(',')[0:4]) for line in file_contents ]
    lst = [ QgsColorRampShader.ColorRampItem(line[0], QColor(line[1],line[2],line[3])) for line in file_contents ]
    return lst

raster = "D:\\example.tif"

fileInfo = QFileInfo(raster)
path = fileInfo.filePath()
baseName = fileInfo.baseName()

layer = QgsRasterLayer(path, baseName)
QgsMapLayerRegistry.instance().addMapLayer(layer)

file_contents = """# QGIS Generated Color Map Export File
INTERPOLATION:INTERPOLATED
0,0,0,0,255,0.000000
1,0,255,0,255,1.000000
2,0,170,0,255,2.000000
3,170,85,0,255,3.000000
4,243,205,255,255,4.000000
5,0,255,255,255,5.000000
6,255,255,0,255,6.000000
7,255,255,255,255,7.000000
8,255,0,0,255,8.000000
9,0,0,255,255,9.000000
10,130,133,123,255,10.000000
11,11,86,29,255,11.000000
12,255,0,255,255,12.000000
13,0,0,0,255,13.000000
14,255,174,0,255,14.000000"""
lst = GetColorRampItemListFromText(file_contents)

fcn = QgsColorRampShader()
fcn.setColorRampType(QgsColorRampShader.INTERPOLATED)

fcn.setColorRampItemList(lst)
shader = QgsRasterShader()
shader.setRasterShaderFunction(fcn)
renderer = QgsSingleBandPseudoColorRenderer(layer.dataProvider(), 1, shader)
layer.setRenderer(renderer)

layer.dataProvider().reloadData()
layer.triggerRepaint()

My question is either:

Why is my code not rendering the color table?
Is there an easier way to apply a color table from a file?

EDIT:
The reason it wasn't rendering has something to do with layer.dataProvider().reloadData(). Removing it works as expected.

For what it's worth, Here's another similar question and answer: Assigning Color Ramp using PyQGIS


texts:
We are calculating Population and geocoverage in mapinfo. We have calculated population coverage as an area overlap. How do we calculate population coverage if we have population defined for each polygon?

1.  Are there any particular method to calculate population coverage for each boundary?

2. Any documentation and detail for analysis (not only steps but reason) would be helpful.


texts: Editing validation code for Python script tool in ArcGIS Pro using IDLE?

texts:
I created and successfully tested a Python script tool, for which I included some validation code, in ArcGIS Pro 1.3.1.

To write the validation code I right-clicked on the Python script tool in its standard toolbox to open its Tool Properties, went to the Validation tab and clicked Open in Script Editor.

This opened IDLE, allowed me to edit the code, and after I saved and closed IDLE, clicked OK to close the Tool Properties, and ran the tool the validation that I had just set clearly worked.

Five minutes later I went to enhance the validation code using the same editing workflow but found that clicking Open in Script Editor no longer opened IDLE to allow me to edit it.

I have rebooted my machine and the Open in Script Editor button is still no longer opening IDLE.

Via comments at the time I determined that I could get a very suboptimal script editor in Notepad to work by using the Ribbon to access ANALYSIS | Geoprocessing Options to set the Script Editor to C:\Windows\System32\notepad.exe.

However, doing the same thing to try and set the Script Editor to C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\lib\idlelib\idle.pyw returns the Open in Script Editor button to doing nothing.

I few weeks later I moved to another machine which had been upgraded to ArcGIS Pro 1.4.1 and I was hoping that would fix whatever was astray but it has not.

Is there any other method that can be used to be able to always edit the validation code for a Python script tool of ArcGIS Pro by using a Python IDE better than Notepad?

I have considered re-writing the tool as a Python Toolbox (*.pyt)but the site I am working with are more comfortable using Standard Toolboxes (*.tbx) so I would like to avoid doing that.


texts: ST_Polygonize/Shapely cannot polygonize when overlapping with shared node/point


texts: Populate field with information from a different field based on matching record

texts:
I have a feature class created by merging multiple feature classes together, after the merge there are a lot records that have the same spatial attributes but different tabular attributes, except one field.

I want an arcpy script that will look for matching records in 'field 1'(there are two or more matches) and for the first match it should be the full value of 'field 2', for the second match it should be only the last three characters.

Example:
'field 1' ('subcatchme') has matching records = '68103' and '68103' the values in 'field 2' ('NAME') are 'W18_0_02' and 'W18_0_03' respectively. I want to populate both the matching records in 'field 3' ('CatchmentName') with 'W18_0_02_03'.

from arcpy import env
env.workspace = r"D:\Users\Test2.gdb"

origFC = "W18_0_Merge"
dissolvedFC = "W18_0__Dissolve"

##Get a list of NAME Field
idList2 = []
with arcpy.da.SearchCursor(dissolvedFC, ["subcatchme"]) as cursor:
    for row in cursor:
        idList2.append(row[0])
del cursor
##Remove 'u' from front of list
idList = [i.encode('utf-8') for i in idList2]

##iterate through each subcatchme
for id in str(idList):
    print "Hi"
    ## Use a search cursor to query only the individual subcatchme
    with arcpy.da.SearchCursor(origFC, ["NAME"], '"subcatchme = {0}".format(id)') as cursor:
        print "Hi2"
        for row in cursor:
            print "Hi3"
            ## Use an update cursor to update the dissolved feature class fields
            with arcpy.da.UpdateCursor(dissovledFC, ["CatchmentName"], '"subcatchme = {0}".format(id)') as cursor2:
                for row2 in cursor2:
                    if row2[0] == None or row2[0] == '':
                        row2[0] = row[0]
                    else:
                        row2[0] = row[0] + row2[0] ##I want it to only show last three characters here
                    cursor2.updateRow(row2)
    del cursor, cursor2

I have re-purposed code previously pointed to me (Here) but I keep getting an error after the second cursor. the "Hi2" will print but it's not going through it.

  RuntimeError: An invalid SQL statement was used. [SELECT NAME,
  OBJECTID FROM W18_0_Merge WHERE "subcatchme = {0}".format(id)]

I have also been trying to change the last cursor so it only grabs the last three characters but the only format I can find is slicing strings but I don't know how to put it together when 'row[0]' already has the square brackets.

s = "Pirates of the Caribbean"
print(s[21:24])

would result in 'bean'

I am using ArcMap 10.4 with all the extensions. the program I am writing python in is pythonwin.


texts:
I am trying to perform a raster calculation of a DEM and I would like it to have any desired resolution of my choice. How can I write the code so that I adjust the pixel height and width in feet using any number from the interactive dropdown as shown?

Also, how can I avoid the message error below? I would truly be grateful for any practical solution provided here.
Message Prompt: AttributeError: 'NoneType' object has no attribute 'isValid' .

Here is my code below:
# pixel calculator layer (if displayed and checked)
if layer.type() == QgsMapLayer.RasterLayer and self.dlg.checkPixel.isChecked():

    # Get layer object from an a just created layer
    fn = 'output1'
    # rlayer = iface.addRasterLayer(fn)
    rlayer = QgsRasterLayer(fn, '')
    lyr = rlayer

    # get extension about the raster
    filename, file_extension = os.path.splitext(layer.source())
    output2 = self.folderName + "/ThinDEM/rasCal_" + layer.name() + file_extension

    # check file isn't openned and is writable
    version = 0
    while self.isFileOpened(output2):
        output2 = self.folderName + "/ThinDEM/rasCal_" + layer.name() + "("+ str(version) + ")" + file_extension
        version +=1

    entries = []
    ras = QgsRasterCalculatorEntry()
    ras.ref = 'ras@1'
    ras.raster = lyr
    ras.bandNumber = 1
    entries.append( ras )

    calc = QgsRasterCalculator( '(ras@1 + ras@1) / 2.0', output2, 'GTiff', lyr.extent(), lyr.width(), lyr.height(), entries )
    calc.processCalculation()

    # load layer
    if self.dlg.checkBox_6.isChecked():
        out = iface.addRasterLayer(output2, "")
        if not out.isValid():
            iface.messageBar().pushMessage(self.tr("Error"), self.tr("Could not load ") +  output2, level=Qgis.Warning)


texts: Looking for USDA plant hardiness zone shapefile

texts:
Given that this is a USDA resource I really thought that it would be available online somewhere, but for the life of me I cannot find a shapefile for the USDA plant hardiness zones anywhere.  Does anyone know where there might be one?

This  link is the closest I found but it doesn't seem to offer any sort of shapefile download:

http://www.usna.usda.gov/Hardzone/ushzmap.html


texts: PROJ.4 custom projection that is Transverse Mercator with Affine post-process

texts:
I am trying to define some of my local TMAF grids in PROJ.4 so I can reproject points between the local grid and MGA95 and/or WGS84.

For these grids I have the MapInfo CoordSys clause which defines a lot of parameters.
eg: “CoordSys Earth Projection”, 8, 33, 7, 117, 0, 0.9996, 500000, 10000000, 7, 0.890953, -0.455062, 2903977.24, 0.455063, 0.890954, -6919253.68, -100000, -100000, 100000, 100000

From the MapInfo doco I can see these parameters align to Datum, Origin Latitude, Scale at Origin, False Easting/Northing, and the Coefficients for the Affine Process.

What I'm having trouble with is trying to find out how to use these coefficients when creating a new custom projection, and the parameters in the PROJ.4 string that they could correctly map to.


texts: How to add date and time to a geotiff to enable time dimension in geoserver?

texts:
I have a collection of geotiff files, each one corresponding to a date YYYYMMDD.

How can I had the information date to those geotiff, so that geoserver allows the enabling the "Time" dimension for this layer?


texts: Data is not correctly appending in cloumn

texts:
for pt in pointLayer.getFeatures():
    for pg in polygonLayer.getFeatures():
        if pg.geometry().intersects(pt.geometry()):
            intersection = pg.geometry().intersection(pt.geometry())
            intersectionLength = intersection.length()
            intersectionArea = intersection.area()

            fieldIndex = pointLayer.fieldNameIndex(pointLayerColumName)
            print  "LOT Value:: " + str(pt[pointLayerColumName])
            if bool(pt[pointLayerColumName] == None):
                updateValue = pg[polygonLayerColumToUpdate]
                print  " 1 updateValue : " + updateValue
            else:
                updateValue = str(pt[pointLayerColumName]) + ","+ pg[polygonLayerColumToUpdate]
                print  "2 updateValue : " + updateValue

            pointLayer.startEditing()
            pointLayer.changeAttributeValue(pt.id(), fieldIndex, updateValue)
            pointLayer.commitChanges()

I want to update a column but first want to check that if data present in that column then add comma (,) after it and append my new data in their. But the problem is only last (looped) data is appending .

And this line

 if bool(pt[pointLayerColumName] == None):

always return true although I previously update data.


texts: Architecture for serving OSM-Mapnik Tiles

texts:
I am currently trying to set up a Mapping-Service for serving OSM-Mapnik Tiles. In my project we want to use the Mapnik-Layer as a BaseLayer for our GI-Application. I am quite new to Mapnik (only used a public Mapnik-Service as a BaseLayer in OL so far) and i am now struggling with the decision about a good service architecture for my use case. Maybe somebody could give me a hint.

Currently from my point of view i have to decide between two possible service architectures:

1.) In the first solution i serve mapnik tiles via mod_tile. The OSM data lies in a PostGIS database and is updated regularly. Mod_tile in combination with rendered does than the rendering and tiling of the OSM data and serves the tiles via tileserver interface. This data could than be integrated in an OL-Application via OpenLayer.Layer.OSM.

2.) In the second possible service architecture the OSM data again lies in a PostGIS database and is updated regularly. But the data is now served via a wms interface (are there any good solutions for this?). The wms interface would then be cached by a separate mapping cache. After that the maps could be used in the OL-Application via wms or wmts interface.

From my initial point of view, the second solution would better fit in my existent architecture, where i already use a mapserver which is cached by a mapcache. Also if using wms and osm layers together in Openlayers, i always get some troublesome behavior of the map extent, where the mapnik tileset define the map extent. The problem is here, that the mapnik map extent is always quite bigger than the actual visible map extent. This also speaks for a wms / cache solution, where the behavior is from my point of view a bit more like i expect it.

On the other hand i am new to the mapnik world and so i am not sure if this solution does better fit to the osm data and idea. So what is your opinion and experience? What is are benefits of the first architecture compared to the second one? What are good software tools (in a debian environment) for publishing osm / mapnik data via wms interface?

I hope i make my point clear and some people could share there experience with me.


texts: Carto isochrones funtion returns syntax error

texts:
Using this snippet of code:

INSERT INTO my_new_table
    (the_geom) SELECT (cdb_isodistance(the_geom, 'walk',
    string_to_array(distance, ',')::integer[])).the_geom
FROM my_original_pts_table

I am trying to create isochrones for a table of points.  I get an error after applying an INSERT.

Why is the first geom in parenthesis? Is the distance parameter something that will vary from table to table? After this does work, I still have to wrap it in a ST_Transform correct?


texts: Why is geography faster than geometry for distance calculations?

texts:
Intuitively, I thought geometry would be faster than geography distance calculations, but this says otherwise: https://stackoverflow.com/questions/3547699/need-performance-on-postgis-with-geodjango#answer-3551142

The docs here also seem to suggest that distance calculations are more efficient using geography:
https://docs.djangoproject.com/en/1.9/ref/contrib/gis/model-api/#srid

Why is this? Is it some sort of optimization or indexing with geography that allows faster distance calculations?


texts:
I'm working on a geocoding project using the plugin in MMQGIS QGis, and I want to use OpenStreetMap data in this project.

I have already downloaded the data of my country and I converted to .CSV files, the problem is that the file contain roads and street names (no province, city or country) so how to configure this data to use them?

knowing that the fields of the province, city or country are obligatory to geocode my data.


texts: Overview over the current activities and results to represent GeoData semantically to allow formal reasoning?

texts:
I am currently interested in representing some geo-data in a semantic way via RDF/OWL to allow formal reasoning on a geometry level. However I am lacking a good overview  on the current standards to do so and on the benefits and limitations of each.

If you could recommend a good overview webpage or paper I would be very thankful!


texts:
I am trying to georeference a building using 4 coordinates, which produces the following result:

After georeferencing, shouldn't the image be distorted to match the red lines? The loaded result is as follows:

I am not sure what I might be missing. Both layer and project's CRS are set to WGS 84.

Here is the report table:

Any ideas?

Answer:
Just changed the transformation type to Polynomial 1, which produced the following result:


texts: Creating a selection set from selected features using select by location

texts:
I am trying to convert my VBA code to vb.net. I use “Select by location” to select all the points that intersect lines. Then I want to click a button to put the points into a selection set. Then I can use two other buttons to cycle through the selection set.
When I try to convert my code I keep on getting the following error in vb.net.:
In vb I keep on getting error for ‘pFeatureCursor’ saying “Variable ‘pFeatureCursor’ is passed by reference before it has been assigned a value. A null reference exception could result at runtime.
Here is the VBA code:
Public Sub GetSelSet()
  Dim pMxDoc As IMxDocument
  Dim pMap As IMap
  Dim pActiveView As IActiveView
  Dim pFeatureLayer As IFeatureLayer
  Dim pFeatureSelection As IFeatureSelection
  Dim pSelectionSet As ISelectionSet
  Dim pFeatureCursor As IFeatureCursor
  Dim pCurFeat As IFeature
  Dim intLayer As Integer

  Set pMxDoc = Application.Document
  Set pMap = pMxDoc.FocusMap
  Set pActiveView = pMap
  Set pSelSet = New esriSystem.Set
  CurrentIndex = 0

  For intLayer = 0 To pMap.LayerCount - 1
    If TypeOf pMap.Layer(intLayer) Is IFeatureLayer Then
      Set pFeatureLayer = pMap.Layer(intLayer)
      If Not pFeatureLayer.FeatureClass Is Nothing Then
        If pFeatureLayer.FeatureClass.featureCount(Nothing) > 0 Then
          Set pFeatureSelection = pFeatureLayer
          Set pSelectionSet = pFeatureSelection.SelectionSet
          pSelectionSet.Search Nothing, False, pFeatureCursor
          Set pCurFeat = pFeatureCursor.NextFeature

          Do While Not pCurFeat Is Nothing
            pSelSet.Add pCurFeat
            Set pCurFeat = pFeatureCursor.NextFeature
          Loop

        End If
      End If
    End If

  Next intLayer

End Sub


texts:
I want to add more controls to myOl3ViewTemplate application. If you can help me add:

draw feature
delete feature
zoom[center] at a certain location

I looked at app.js and tried to add more controls but I couldn't get them to work. As for to zoom at a certain location I tried Bethlehem[South Africa] but the center goes to Congo I don't know why because these coordinates[-28.23333,28.3] says Bethlehem[SA] .

I would like you to help me on application launch to view only South-Africa and zoom at Bethlehem.


texts: Finding Square Polygons using ArcGIS Desktop?

texts:
Anyone know of an ArcGIS tool (or trick!) to determine square polygons?

This is sort of a - find the geometric shape of a polygon, if all lines of a polygon are equidistant then YES (select)

Seems easy, yet .. .

Thoughts?


texts: Customizing Distance Units in ArcGIS Desktop?

texts:
Is it possible to add custom distance units (ex. varas) to ArcGIS Desktop 10 so that they can be used when inputting distances in Editor and COGO tools?

I often work with deeds that reference distances in varas. It would be super handy if I could setup the conversion factor and have ArcMap do the math rather than needing to convert all my distances to feet before entering them.


texts:
Style: Rule-Based
 Please be informed that the Add Categories to Rule under the "Refine Current Rules" does not filter in accordance to the selected layer that you have chosen. Instead, all subcategories are being displayed. It does not automatically trim down to specific content as filtered/classified and that makes it harder to select only what is really necessary. I hope that they will make an update to this so that using QGIS will make it more fun to use.


texts: Missing mapsettings.ui while loading QGIS Cloud 3.1.8

texts:
I have an issue with QGIS Cloud 3.1.8. The plugin installs nicely but doesn't load due to the following error message.

FileNotFoundError: [Errno 2] No such file or directory:
'C:/Users/ahmed/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\qgiscloud\mapsettings.ui'

I was able to find the file mapsettings.ui on GitHub.
Should I save it and add it manually? Should I save it in a .txt-format? Is there perhaps a more elegant solution?
I'm not used to working with GitHub.
QGIS cloud was working fine before on the previous version 3.1.7 but when I updated it I had this issue.
QGIS v3.10.7
QGIS Cloud v3.1.8


texts: Downloading NED Data: GridFloat versus IMG versus ESRI Grid (ArcGrid)

texts: What are the advantages/disadvantages of each of these three file types?  They are all available as formats for downloading NED data from The National Map.  Which should generally be downloaded for what purposes?


texts: QGIS and Mapinfo reporting different projections for some files

texts:
When I open a tab file that's been generated in QGIS using EPSG 28355 (MGA 94, z55s) using QGIS 3.10 the projection is transformed to 5551

But when I open it in Mapinfo v17 it's correctly reported as 28355

ogrinfo thinks that it's EPSG:9122 and 9001

This is not just one dataset, it's been happening to any data that starts in 28355 or 28354.

It works fine for shp/gpkg etc...the issue is just with the .tab format.

As per my comment to a reply -this has been happening in all QGIS 3x versions and the test is simple.

Create a temp table in 28355

Add some data

Save it as a tab file and a shp file using save-as rather than "Make Permanent"
In a new project try and bring the data in. The SHP is fine but the TAB asks to transform.

TAB

SHP


texts:
I'm new to Mapbox and I want to add a GeoJSON made of polygons stored locally on my computer and I don't know how to do.
I've follow the example given by Mapbox (https://docs.mapbox.com/mapbox-gl-js/example/geojson-polygon/) but it's not doable to add 100 polygons that way. My file is structured like this:
{ "type": "FeatureCollection",
    "features": [
      { "type": "Feature",
         "geometry": {
           "type": "Polygon",
           "coordinates": [
             [ [100.0, 0.0], [101.0, 0.0], [101.0, 1.0],
               [100.0, 1.0], [100.0, 0.0] ]
             ]
         },
         "properties": {
           "prop0": "value0",
           "prop1": {"this": "that"}
           }
         },
    { "type": "Feature",
         "geometry": {
           "type": "Polygon",
           "coordinates": [
             [ [100.0, 0.0], [105.0, 0.0], [108.0, 1.0],
               [100.0, 1.0], [100.0, 0.0] ]
             ]
         },
         "properties": {
           "prop0": "value0",
           "prop1": {"this": "that"}
           }
         },
       ]
     }

I've tried this but my polygons don't appear on my map. I think i'm doing something wrong on map.addSource.
<script>
    mapboxgl.accessToken = 'my_token';
    var map = new mapboxgl.Map({
        container: 'map', // container id
        style: 'mapbox://styles/mapbox/streets-v11', // style URL
        center: [0, 0], // starting position [lng, lat]
        zoom: 9 // starting zoom
    });
    map.on('load', function () {

        map.addSource('scmpd-precinct-polygons', {
          type: 'geojson',
          data: 'my_file.geojson'
        });
     map.addLayer({
            'id': 'precinct-layer',
            'type': 'fill',
            'source': 'scmpd-precinct-polygons',
            'layout': {},
            'paint': {
                'fill-color': '#088',
                'fill-opacity': 0.8
            }
        });
    });

</script>


texts: QField - Get attribute from selected feature

texts:
My objective is to get a automatic layout with QFIELD based on a selected object (in QField) and with some attributes of the selected object.
The project works fine on QGIS but it seems that the layout of QField doesn't recognize the selected object.
Exemple of code using to get an attribute from the selected object on the layout :
concatenate(CASE WHEN is_selected () is true THEN "name_of_column" ELSE '' END)

An idea?


texts: Add multiple rasters using QgsRasterLayer with for loop

texts:
I am trying to create a for loop for coloring up multiple rasters simultaneously. I am able to do it for one layer, but trying to implement the for loop to do it for more than one layer is proving difficult.

import glob,os

rast_path="/home/user/Desktop/data/gsi/output"
rasters=glob.glob(os.path.join(rast_path,"*le.tif"))
print(rasters)

for raster in rasters:
#    (rastdir,rastfile)=os.path.split(raster)
    lyr=QgsRasterLayer(rasters,rastfile)
    s=QgsRasterShader()
    c=QgsColorRampShader()
    c.setColorRampType(QgsColorRampShader.INTERPOLATED)
    i=[]
    i.append(QgsColorRampShader.ColorRampItem(-20,QtGui.QColor('#0c0101'),'-20'))
    i.append(QgsColorRampShader.ColorRampItem(-15.1,QtGui.QColor('#ff0000'),'-15'))
    i.append(QgsColorRampShader.ColorRampItem(-10,QtGui.QColor('#ff0000'),'-10'))
    i.append(QgsColorRampShader.ColorRampItem(-5,QtGui.QColor('#ff0000'),'-5'))
    i.append(QgsColorRampShader.ColorRampItem(0,QtGui.QColor('#ff0000'),'0'))
    i.append(QgsColorRampShader.ColorRampItem(5,QtGui.QColor('#ff0000'),'5'))
    i.append(QgsColorRampShader.ColorRampItem(10,QtGui.QColor('#ff0000'),'10'))
    i.append(QgsColorRampShader.ColorRampItem(15,QtGui.QColor('#ff0000'),'15'))
    i.append(QgsColorRampShader.ColorRampItem(15,QtGui.QColor('#0c0101'),'20'))
    c.setColorRampItemList(i)
    s.setRasterShaderFunction(c)

    ps=QgsSingleBandPseudoColorRenderer(lyr.dataProvider(),1,s)
    raster.setRenderer(ps)
    QgsMapLayerRegistry.instance().addMapLayer(lyr)
#    qgis.utils.iface.addRasterLayer(raster,rastfile)


texts:
I have a vector layer that i would like to join with a list. The problem is that my vactor layer does not hold any overlapping or dublicate objects, but my list does. The list does not have any geomatrical data but it has an id field that coresponds with that of the vector layer. When i'm using the join feature in the property menu i only get the first item in the list which has a corresponding id in the vector layer. Is there a way to get a new vector layer from the join that includes every item from the list with the corresponding information and polygon from the vector layer?

What i'm looking for is the following:

The list contains no geographical information, but the result does.

Hope this is clear enough, please send me a comment otherwise.


texts:
I am using QGIS to georeference B&W images. I have now done this and would like to set the black surround of the image (as a result of the georeferencing) to transparent but this does not work as bands 1,2 and 3 have the same values. I have noticed there is now a band 4 after georeferencing and would use this to drop out the black surround but I can only do this to bands 1,2 and 3.

I am not super experienced with rasters in QGIS so I am having difficulty determining if this question has been addressed.

This is a screen grab of the image I have, I am trying to remove the black triangles:

These are the values of the black surround:

These are the values of the black text and lines I would like to keep:

For anyone who reads this who might find this usefull I just worked out how to do it. I can set a band to be transparent. So I just set Band 4 to be transparent and it worked perfect.


texts: Finding point with highest elevation in area using ArcObjects?

texts:
Is there a way to find a point with highest elevation in specified area by using ArcObjects? Let's say that I have a raster with DEM and I would like to establish some areas of interest in this raster. In each AOI, I would need to find a point where terrain has the highest elevation and return coordinates of this point.

I looked into ArcObjects API and I found only a manual way to discover such point, that is:

Establish AOIs
For each AOI, find corresponding pixels in the raster
Having sets of pixels defined for AOIs, iterate over each set and return a point with highest elevation

This process is doable but absolutely impractical due to performance concerns.


texts: What does 'Course Made Good' in NMEA string mean?

texts:
I was looking at some GPS data that a client sent me. It has lines like this:

$GPRMC,183123,A,1725.4170,N,07825.4779,E,37.80,297,161214,,

I was looking at http://aprs.gids.nl/nmea/#rmc, which says that the 8th piece of text stands for Course Made Good, True.

What exactly does this mean?


texts: How can I match Geolife traces using TrackMatching or similar tools?

texts:
I am working on route prediction and I decide to use realistic trace of moving object, I have choose geolife dataset
it is composed of 278 users and every user has its trips presented in GPS or GPX I am not sure format in ".plt" file

Geolife trajectory
WGS 84
Altitude is in Feet
Reserved 3
0,2,255,My Track,0,0,2,8421376
0
39.991364,116.326605,0,223,39750.3962731481,2008-10-29,09:30:38
39.991551,116.326653,0,163,39750.3963310185,2008-10-29,09:30:43
39.991821,116.326695,0,75,39750.3963888889,2008-10-29,09:30:48
39.991771,116.326676,0,77,39750.3964467593,2008-10-29,09:30:53
39.991825,116.326657,0,77,39750.3965046296,2008-10-29,09:30:58

I have tested currently trackMatching tool but I get a problem that my format isn't supported

?xml version="1.0" encoding="UTF-8" standalone="yes"?>
<restResponse>
  <error>Invalid GPX format</error>
</restResponse>

I have also tested Project-OSRM but it requires an OSM map which I don't have


texts: Iterating over each pixel of an image in the time-series

texts:
I am using on GEE for a small analysis of Sentinel-1 data collection. I would like to go through each pixel of my AOI and check whether there was a drop in the time-series and classify that pixel accordingly.
Any suggestions?


texts: Reading shp using networkx, parallel edges are missing attributes

texts:
I am trying to use networkx.read_shp to create a MultiGraph from a shape file. I understand that the function needs to be modified since it usually creates a DiGraph I attempt by trying this:
def read_multi_shp(path, simplify=True, geom_attrs=True, strict=True):
    try:
        from osgeo import ogr
    except ImportError:
        raise ImportError("read_shp requires OGR: http://www.gdal.org/")

    if not isinstance(path, str):
        return

    net = nx.MultiGraph() #My first change
    shp = ogr.Open(path)
    if shp is None:
        raise RuntimeError("Unable to open {}".format(path))
    for lyr in shp:
        fields = [x.GetName() for x in lyr.schema]
        for f in lyr:
            g = f.geometry()
            if g is None:
                if strict:
                    raise nx.NetworkXError("Bad data: feature missing geometry")
                else:
                    continue
            flddata = [f.GetField(f.GetFieldIndex(x)) for x in fields]
            attributes = dict(zip(fields, flddata))
            attributes["ShpName"] = lyr.GetName()
            # Note:  Using layer level geometry type
            if g.GetGeometryType() == ogr.wkbPoint:
                net.add_node((g.GetPoint_2D(0)), **attributes)
            elif g.GetGeometryType() in (ogr.wkbLineString,
                                         ogr.wkbMultiLineString):
                for edge in edges_from_line(g, attributes, simplify,
                                            geom_attrs):
                    e1, e2, attr = edge
                    net.add_edge(e1, e2)
                    net[e1][e2][0].update(attr) #Second change and perhaps the one that is not working
            else:
                if strict:
                    raise nx.NetworkXError("GeometryType {} not supported".
                                           format(g.GetGeometryType()))

    return net

The rest of the function stays the same. When checking eges I can see that when there are parallel edges one of them seems to have an empty dictionary instead of the attributes (despite these being present in the shapefile). I have seen a similar question before but it seems to be using an outdated version of networkx. Is there any way to make a multigraph that has parallel edges with their attributes?


texts: How do I use a Simple Feature Collection (sfc) from the sf package with the elevatr package to load elevation data?

texts:
I'm trying to load elevation data from AWS using the elevatr package, with a region defined as a Simple Features Collection from the sf package. Here's a simple example:

library(elevatr)
library(sf)
library(USAboundaries)

counties <- us_counties(map_date = "1930-01-01", resolution = "high", states = c("NY"))
elevation <- get_elev_raster(counties$geometry, z = 9, src = "aws")

This gives an error:

Error in if (attributes(class(locations)) %in% c("raster", "sp")) { :
  argument is of length zero
In addition: Warning messages:
1: In if (class(locations) == "data.frame") { :
  the condition has length > 1 and only the first element will be used
2: In if (class(locations) == "SpatialPoints") { :
  the condition has length > 1 and only the first element will be used
3: In if (class(locations) == "SpatialPointsDataFrame") { :
  the condition has length > 1 and only the first element will be used

Presumably because the sfc is a collection of polygons and MULTIPOINT objects, which aren't supported by the get_elev_raster function of the elevatr package.

Is there any way to use these packages together? With the end result being an sfc object that also has various forms of elevation data for each point?

I'm getting the data this way because it doesn't seem like the st_read function from the sf package can read ESRI ASCII files.


texts: ArcGIS Service that returns various map images given coordinates?

texts:
I just started learning about ArcGIS.

I am building a webapp (django) and we need map images (jpg) of a given location (lat,long), using various basemaps (topology, satellite, etc). Later I need to layer shapefiles on top.

The Google static image API does exactly what I need. Unfortunately, they only provide 3 different types of maps and we need other basemaps.

ArcGIS has all the basemaps we need. But, how do I use it?

I can see they have an ArcGIS API for Python and an ArcGIS REST API. But, browsing those docs I cannot find a way to generate map images.

Basically, I need a URL that implements the following function: "Generate a .jpg of map range lat,long to lat1,long1 using Basemap XYZ (and layer shapefile A.shp on top)", where I can tell it what lat,long, XYZ and A.shp I want. Does such a thing exist? in arcgis, or anywhere? Paid services are OK.


texts: Why is arcpy overwriting datetime in __main__? How do I fix it without rearranging my imports?

texts:
Consider the following scripts:

import_test1.py:

from datetime import datetime
print(u'{}: {}'.format(__name__, datetime))
import import_test2
print(u'{}: {}'.format(__name__, datetime))
from datetime import datetime
print(u'{}: {}'.format(__name__, datetime))

import_test2.py:

from datetime import datetime
print(u'{}: {}'.format(__name__, datetime))
import arcpy
print(u'{}: {}'.format(__name__, datetime))
from datetime import datetime
print(u'{}: {}'.format(__name__, datetime))

Now run this script:

C:\temp>py -2 .\import_test1.py
__main__: <type 'datetime.datetime'>
import_test2: <type 'datetime.datetime'>
import_test2: <type 'datetime.datetime'>
import_test2: <type 'datetime.datetime'>
__main__: <module 'datetime' (built-in)>
__main__: <type 'datetime.datetime'>

This is an MCVE. (I encountered the issue in a more complex script, but this suffices to reproduce the behavior.)

As you can see, importing arcpy (not doing a from arcpy import *, just a plain import arcpy) is overwriting my module's import when it's the __main__ script. It replaces the type I imported with the module of the same name, which is not what I need.

I know it's arcpy because I can reproduce it with an even simpler script:

from datetime import datetime
print(u'{}: {}'.format(__name__, datetime))
import arcpy
print(u'{}: {}'.format(__name__, datetime))
from datetime import datetime
print(u'{}: {}'.format(__name__, datetime))

Gives:

C:\temp>py -2 .\import_test.py
__main__: <type 'datetime.datetime'>
__main__: <module 'datetime' (built-in)>
__main__: <type 'datetime.datetime'>

I don't want to reorder my modules because I want to follow PEP8's guidelines on order of imports: standard library, third party, local application/library specific imports. This is also hard to manage if some other module performs the import of arcpy, as in my two-module example.

How can I get arcpy to stop overwriting variables that don't belong to it? Why is it even doing this in the first place?


texts: Openlayers plugin not displaying any source

texts: I was previously using QGIS 2.2 and using the openlayers plugin to display google streetmaps on an area. Reopening the project several months later I found that this layer no longer displays. In fact it displays as blank and white. If the google streetmaps is the top layer it blanks out all the layers below it. None of the other sources on openlayers plugin displays either. Usually, when selected, the map appears with a grid with little question marks in the middle of each square. I uninstalled 2.2 and installed 2.4, 2.6 and 2.8 in sequence, at each stage removing traces of the previous installation. However I can never get anything on the openlayers plugin to appear on any map. I am running windows 7


texts:
I know this kind of problem as been posted previously, but I can't solve my problem with the answers.

I'm running a simple gdal_merge command, but I get shut down every time by permission problems. I'm administrator of the directories, running CMD as Administrator. My account has writing permissions over all the folders I use for this. I'm running on Windows 7, 64bits.

Here is the command:

gdal_merge -o "C:\Users\Catherine Plante\Desktop\SuiviHydroGranby" --optfile "\\Server\DONNEES\MNE\liste.txt"

All I can think of is that GDAL (through OSGeo4W) runs the commands as another user, and that this one is blocked.

Does this make sense?

EDIT:
In the end, I had forget to clearly identify the output format. While doing so, I realized that I also had problems with the native format (.img). So this was mostly my fault! Just sad that GDAL identified the problem as a permission problem, rather then an incomplete query.


texts: QGIS - Counting vertexes in radius from each feature in line layer

texts:
I am looking to measure the impact of urban green on pedestrian route choice, and in that connection I want to create a new field within a line layer that counts how many trees exist in a radius of 10 meters from the center of each feature in the line layer. (could also be within the radius of 10 meters from the two extremes of the line, doesn't matter so much).
I am working in QGIS 3.4.13-Madeira.
Am at a total loss as to how to do this, and can't seem to find any articles that might point me in the right direction.

(the points are trees)


texts: Is there a way to according to the data's spatial reference to get the unit of data, for example: wgs1984 → Decimal degrees, UTM54 → Meters.


texts: Create normal vector and calculate distance between start point and intersection

texts:
I have one line-layer which pictures a street net and a point-layer which defines certain locations. Please see in the picture below:

What I am trying to do is to create a normal vector through a point and its nearest line for each point. Then I would like to calculate the distance between the line's starting point and the intersection of the normalvector. This should look like this:

As I have lots of points, I need to find a way to solve this generically.

I know that the plugin NNJoin gives me the nearest edge to a point, but how do I create normal vectors and how can I then calculate the lengths?

Do I have to write a python script to realize this or is there an other way?


texts: Keeping unmatched records using Join Field in ArcMap

texts: I would like to keep all records when joining fields in ArcMap.  I know that I could do this by adding a join, but I would like to use the Join Field tool so that I can dictate which fields are joined and so that it does not change the field name when joined.  Is there a way to keep unmatched records when using the Join Field tool?


texts: How to install GDAL 1.10 with OGR SQLite driver support on Ubuntu?

texts: How can I install GDAL 1.10 with OGR SQLite driver support, and JPEG and PNG drivers per http://www.gdal.org/frmt_mbtiles.html


texts: Shape length distance wrong in Robinson projection

texts:
I have a world map in the WGS 1984 Geographical Coordinate System. I want to measure the lengths of borders between countries in kilometers. To get the units in meters, I have first used the Define Projection tool and then the Project tool to change to the Robinson Projected Coordinate System. After that, I am using the Intersection tool and setting the output type to "line" to get a Shape_length column where I will have all the borders and their lengths. But, the value in the Shape_length column appears to be wrong as it measures the borders to be only a few meters. The Measure tool does the same thing.

I can't get my head around what I am doing wrong. I've attached a print screen of the data frame. What could the problem be?


texts: Loading Geodataframe to PostgreSQL with psycopg2 connection

texts:
I have created a Geodataframe from a PostgreSQL table with the following instruction (where 'con' is a psycopg2 connection):

gdf = gpd.GeoDataFrame.from_postgis("SELECT id, geom FROM {0}".format('pg_table'), con=con, geom_col='geom')

I have then processed the Geodataframe and I would like to load the result into a new table in PostgreSQL. However, I don't seem to find an equivalent instruction to 'from_postgis' (something like 'to_postgis' or similar). I have read a post about some SQLAlchemy, which seems over-complicated and I must admit I have not understood. Just wondering whether there's any simple way to do what I'm trying to do?


texts: Deleting feature classes in file geodatabase when context menu is disabled?

texts:
I have a file geodatabase with 20 feature classes and I need to delete a few of them, but the option to DELETE in ArcCatalog, if I select this feature class, and click right mouse, is disabled for the most of the feature classes in the file geodatabase.

Why is the reason for that? Can I fix it? This file geodatabase is created by xmlscheme of gdb.


texts:
"Natural Resources Canada (NRCan) has released the Canadian Geodetic Vertical Datum of 2013 (CGVD2013), which is now the new reference standard for heights across Canada. This new height reference system is replacing the Canadian Geodetic Vertical Datum of 1928 (CGVD28)"
-- http://www.nrcan.gc.ca/earth-sciences/geomatics/geodetic-reference-systems/9054

Should we GIS people who mainly work in feet or meters care, or is this something only really of interest to high precision (sub-centimeter) professions like land surveyors and geodesists*? If we do care, how do we use make use of it? I haven't found any .PRJ/.WKT/... reference files such as EPSG:5713, the CGVD for 1928.

* don't know the term for someone who studies and practices geodesy


texts: I have a PDF file that my partners gave me, this file has the exact vertex of the study area and heights lines, so I need to overlap it with some satellite image, so I can compare with the map in paper and make polygons and paint a geological map, the thing is I can't add the data in form of vectors because it appears all flipped, west in the east and north in the south, how can I load a pdf file as a vector layer without this issue?


texts: Calculating average monthly rainfall from raster using shapefile

texts:
I've been working with some raster data (19 years of precipitation data) and a shapefile. The raster extent represents a square-shaped area and the shapefile is within the extent of the raster.

I want to obtain the average monthly rainfall, I mean the average rainfall for each month for all the time series.

Each month has 19 relevant rasters (from the time series period) and I can get the average rainfall which results to be a single raster for each month but don't know how to get the areal average rainfall (a single value that represents the rainfall within the shapefile).

I loaded the rasters as single objects in R. The raster function mean(r1, r2, r3) gives a raster as output, I want a single value that represents only the area of study.

So far I've been working with R and was wondering if I can get that calculation using the raster package or any other package.


texts: I have a polyline feature with all the railway tracks of a country. However, there are many discontinuities in the polyline (maybe 15 feet), that is, at several locations the polyline has a break of 15 feet. How do I make this polyline continuous such that I can run network analyst tool with accurate results.


texts: Why does a z factor of 1 get turned into 0.999998 in my output

texts:
Just built a model that uses 3D Analyst - Surface Volume.  I am trying to use a z factor of 1.  Whenever I run the model my output table shows that a Z factor of 0.999998 was used instead of 1.  I know this is a very small difference but even so, it is skewing my results somewhat.  Really what I would like to know, is there a good explanation of why it gets changed to 0.999998 or is there a way to make it stay at 1?

Thanks again for your input.  When I've tried changing the type of raster, the tool grays out the z factor all together and locks it in at 0.999998. So... I'm guessing there is really no way around this using integer or floating point data.  I even tried converting to integer and using the raster calculator to multiply the data by 1000 to preserve the decimals in integer form, but was still stuck with the 0.999998 z factor. No dice.


texts: How to create a printable data sheet when using edit widget "value realtion" in QGIS 2.0?

texts:
I created a spatiaLite Database with a point layer (every point is a tree). Additional I created several tables to fill the columns of the spatial layer with attributes using the editing widget “value relations”. I need this because I want to speed up data entry! For some other attribute fields the data entry occurs via “line edit”.
I also designed a custom form like Nathan Woodrow describes it in his blog. This works pretty good.
Now my problem is not the data entry but the data output. I want to print a custom data sheet with all my attributes that occurs to one point (tree) that I’ve selected before. Can I do this with the print composer (it should also contain a small map with the point an the google satellite layer)?
When I tried to realize this I inserted an attribute table by choosing the relevant attributes. For those attributes where I’ve used the edit widget “value relation” the entry in my table is the entry of the key column in my external table. But I need the entries in the value column.
Maybe there is a better way to get this work done?
Thanks for any help,
Patrick


texts: Looking for fastest solution for Point in Polygon analysis of 200 million points

texts:
I have a CSV containing 200 million observations with the following format:

id,x1,y1,x2,y2,day,color
1,"-105.4652334","39.2586939","-105.4321296","39.2236632","Monday","Black"
2,"-105.3224523","39.1323299","-105.4439944","39.3352235","Tuesday","Green"
3,"-104.4233452","39.0234355","-105.4643990","39.1223435","Wednesday","Blue"

For each set of coordinates (x1/y1 and x2/y2), I want to assign the US Census Tract or Census Block that it falls within (I downloaded the Census tract TIGER shapefile here: ftp://ftp2.census.gov/geo/tiger/TIGER2011/TRACT/tl_2011_08_tract.zip).  So, I need to do a point-in-polygon operation twice for each observation.  It's important that the matches be very accurate.

What is the fastest way to do this, including time to learn the software?  I have access to a computer with 48GB of Memory -- in case that might be a relevant constraint.

Several threads recommend using PostGIS or Spatialite (Spatialite looks easier to use -- but is it as efficient as PostGIS?).  If those are the best options, is it imperative to populate an Spatial Index (RTree?)?  If so, how does one do so (e.g. using the Census Tract Shapefile)?  I would be extremely grateful for any recommendations that include example code (or a pointer to example code).

My first attempt (before finding this site) consisted of using ArcGIS to do a spatial join (x1/y1 only) of subsample of the data (100,000 points) on US Census Block.  That took over 5 hours before I killed the process.  I'm hoping for a solution that can be implemented on the entire dataset in under 40 hours of computing time.

Apologies for asking a question that has been asked before -- I've read through answers, and I'm left wondering how to implement the recommendations.  I've never used SQL, Python, C, and have only used ArcGIS once before -- I'm a complete beginner.


texts:
I am trying to get QGIS to work with the Anaconda distribution of Python but I do not know how to get this to work.
In an answer to Running QGIS through Python it was mentioned that QGIS does not work with Anaconda and I wanted to double check this, just in case there have been any updates in the last few months.
If it is possible to do so, could you direct me a resource that explains how one can go about doing this?
I am using Windows 7 (64 bit) with conda 3.7.1 and python 2.7.6.


texts: How to get square metric pixels for a GeoTiff, whose pixel spacing is given in degree?

texts:
I have received a GeoTiff file, whose pixel spacing is provided in degree (e.g. +1° to the right, -1° down). Of course, this means that the actual, metric pixel spacing changes in relation to the latitude, which is bad, if you want to relate pixel distances to actual, real-world distances. What is the straight-forward procedure to transform the image in such a way that the result is squared, metric pixels (e.g. +1m to the right, -1m down)?

For sake of completeness: The original image is provided in WGS84-latlon coordinates (EPSG:4326), the final image should be projected to the proper UTM zone (e.g. EPSG:32632).

I figure I'll probably first have to reproject the dataset. But then still the question remains about how to resample the resulting, reprojected image so that its pixel size can be set to an arbitrarily chosen, square size?

Using any distance measuring tool within QGIS is not an option, as I need to export the image again as a GeoTiff in order to process it externally with a different software.


texts: QGIS not binding value from field to widget in my custom form (Qt Designer)

texts:
I have created a custom form in QT Designer with a QSpinBox and a QLineEdit.
I named the objects the same as the column names so that qgis can bind the values.

Added the .ui to a layer.

The widgets are displayed with the value defined in the value field in QT Designer:

and I can edit the value. If I accept and reopen the form I always see the value in QT Designer but not the value I previously edited.

There must be something I am missing.

I am using QGIS 2.12. I have to use this version for now.


texts: Creating point with geometry to insert in database

texts:
I have a PostGIS table of points. With my php code and pdo, I want to add a point in it. I have the coordinates, I want to get the geometry, with ST_GeomFromText. Here is the code :

$x = $_POST['x'];
$y = $_POST['y'];
$geom="ST_GeomFromText('POINT($x $y)',2154)";
$r = $bdd->prepare("INSERT INTO ouvragesp (id, yl93, xl93,annee,geom) VALUES(:identifiant , :y , :x , 2017, :geom)");
$r->execute(array(
  'identifiant' => $id,
  'x' => $x,
  'y' => $y,
  'geom' => $geom));

However, when I execute it, it tells me that :

PDOStatement::execute(): SQLSTATE[XX000]: Internal error: 7 ERROR: parse error - invalid geometry HINT: "ST" <-- parse error at position 2 within geometry

I don't see why it takes only "ST" part.

Does someone know where the problem comes from and how to solve it ?


texts: Installing qgis-ltr or qgis-ltr-full with OSGeo4W or both?

texts:
I need some help with the installation of QGIS using OSGeo4W. I am using windows 7 and I want to install a 64 bit ltr release with all relevant dependencies and plugins.

I choose to install qgis with OSGeo4w to be able to update more conveniently, but I run into problems with the installation. So I do not want to use the standalone installer any more.

I choose the advanced install option in the first step to be able to install the ltr release. In the "select packeges" step under "desktop" there are a qgis-ltr option and a qgis-ltr-full option. I understood that I would need to select the qgis-ltr-full option to get all the relevant dependencies.However, my problem is that the installer shows that qgis-ltr-full is on version 2.8.2-1 while qgis-ltr is on version 2.8.9-1.

What do I need to choose here?

I do want the latest update of the ltr with all the bugfixes, but I also do not want to pick all the dependencies and plugins by hand.

Do I need to install both?
Do I get the latest ltr version that way with all dependencies?
Or do I end up with two separate qgis installations when I choose both?

I did google for an answer but couldn't find any.


texts: Missing GRASS image processing modules in QGIS 2.14 Processing

texts:
I am running QGIS 2.14 on debian. I have activated GRASS tools from Processing menu. Some functionality of the GRASS such as image processing modules are missing (i.evapo.pm and ...). how do I get/activate them?

The independent grass installation have the modules. But I am not sure if they work without any modification or how to integrate them into QGIS.


texts: $area is NULL for over-wrapping geometry

texts:
This question follows from $area gives NULL for some polygons in QGIS's Field Calculator. New question asks the same question as to why some area variables appear NULL. But the problem was resolved for one layer but not the other.
All $areas are found except for the polygon in the screenshot. The larger polygon and smaller polygon are over-wrapped but are separate entries. Area for smaller entry is found but not the bigger one. I tried fixing geometries but it does not help. Please let me know how I can proceed.
I am using the layer, all economic concessions (English) from the link: https://data.opendevelopmentcambodia.net/dataset/economiclandconcessions


texts: Carto SQL API: generate a point via query

texts:
Is it possible in Carto to create a point via a simple query?

If I run: https://UserName.carto.com/api/v2/sql?format=GeoJSON&q=SELECT ST_MakePoint(-50.10, 42.31)

Encoded:

https://UserName.carto.com/api/v2/sqlformat=GeoJSON&q=SELECT%20ST_MakePoint(-50.10%2C%2042.31)

It returns the following error: {"error":["column \"the_geom\" does not exist"]}.

Is there a way to create a point via a URL call?


texts: Creating sequential numbers in sorted field using ArcGIS Field Calculator?

texts:
I have a sorted field that I would like to populate with sequential numbers.
They are not in order of OID.

Is there any way to do this?


texts: Converting 3d polygon feature class into MultiPatch feature class using ArcGIS for Desktop?

texts:
I've got a polygon feature class with z values enabled. The polygon entities represent billboards. I'm looking for a way to convert those polygons into MultiPatch feature class entities. I'm using ArcGIS Desktop, version 10.3.

I've tried to use 'Layer 3D To Feature Class' Geoprocessing tool. I've inserted the required parameters which are the 3d layer (not feature class!) as input, an output multipatch feature class and a grouping definition by the OBJECTID field. It didn't work: my output had the same amount of entities only they all were having the same geometry.

How do I to convert polygon entities with z values (3d) into multipatch feature class?


texts:
I am working with MapInfo files.
I have to make a lot of geometrical changes and each geometry change needs a corresponding change in multiple layers of tabular data.
I am looking at automating the same in QGIS.
for e.g. I want to drag and drop point geometry data from one point to another, edit line geometry data by adding a new vertex, or move line from one place to another. And during these user actions, I want to update the MapInfo tabular data automatically.
How can I do this in QGIS customizations?


texts: Seeking open source point cloud software?

texts: Just acquired a Geoslam Zeb-Revo for scanning our stockpiles and wondered if there was another option to process point cloud data?


texts: Leaflet plot map.on(zoomend) function into addEventListener(function)

texts:
Continuing the previous query:

Leaflet set zoom for clickable layers on the map

I would like to fix my zoom feature, which doesn't work correctly as per in the picture below:

for this purpose I used the following solution:

      document.querySelector("input[name=cf]").addEventListener('change', function() {
            if(this.checked) map.addLayer(job)
              else map.removeLayer(job)
              map.on('zoomend', function() {
              if (map.getZoom() <7){
            map.removeLayer(job);
            }
              else
            map.addLayer(job);
            })
            })

Where I tried to bound the addEventListener('change', function() { with map.on('zoomend', function() {.
Console says nothing, but map stands as usual with no reaction.

The code:

    document.querySelector("input[name=cf]").addEventListener('change', function() {
            if(this.checked) map.addLayer(job)
              else map.removeLayer(job)
            })
             map.on('zoomend', function() {
              if (map.getZoom() <7){
            map.removeLayer(job);
            }
              else
            map.addLayer(job);
            })

also is not correct.

I tried also another implementation:

   document.querySelector("input[name=cf]").addEventListener('change', function() {
            if(this.checked) map.addLayer(job)
              else map.removeLayer(job)
            })

    document.querySelector("input[name=vm]").addEventListener('change', function() {
            if(this.checked) map.addLayer(job2)
              else map.removeLayer(job2)
            })

    document.querySelector("input[name=bt]").addEventListener('change', function() {
            if(this.checked) map.addLayer(job3)
              else map.removeLayer(job3)
            })

and the zoom code underneath:

  const zoomJobsMap = new Map([
   [6, job],
   [7, job2],
   [8, job3],
   ]);

  map.on('zoomend', function() {
let zoom = map.getZoom();   //# get zoom value at once
for (var [zoomLimit, jobObj] of zoomJobsMap) {
    (zoom < zoomLimit)? map.removeLayer(jobObj) : map.addLayer(jobObj);

}
 })

as per the example below:

https://codereview.stackexchange.com/questions/231557/setting-zoom-for-layer-disappearance-in-leaflet/231564#231564

however here the result is even worse (picture below:), because as zoom changes, the placemarks are shifting across the map.

Moreover  the console says:

  Layer.js:159 Uncaught Error: The provided object is not a Layer.
      at NewClass.addLayer (Layer.js:159)
      at NewClass. ((index):986)
      at NewClass.fire (Events.js:190)
      at NewClass._moveEnd (Map.js:1206)
      at NewClass. (Map.js:1654)

Is there any option to fix this zoom in temrs of the checked layer boxes?
Or alternatively shall I use try to uncheck each layer at the very beginning like here:

https://stackoverflow.com/questions/17732704/how-to-make-the-checkbox-unchecked-by-default-always


texts: Clipping TileLayer with georeferenced polygon (clipping mask)

texts:
I would like to add a mask to clip layers (Tile) below in a way that only the part bounded by a ol.Feature (polygon or multipolygon) is visible and the rest is hidden (left white).
The polygon is stored in a DB so I need a way to do that programatically.

Does anybody have some clues or examples to achieve this?

ol.Overlay?
ol.layer.Vector filled with white outside of the features

Ideally:
- add / remove vector features and interactions should update the mask
- swichting layer below the mask is possible

I'm not looking for a working code but more inspirations and ideas to start doing that in a "ol3" way.


texts: Generating seamless contours across tiles from TIN/LiDAR data?

texts:
I have both TINs and LiDAR data for an entire county, divided into a grid of multiple tiles. I want to generate contours from this data at 2 foot intervals, using ArcGIS 10.5.

I decided to do this by first generating an elevation raster using TIN to Raster, and then running the Contour tool to create contours from the raster.

The problem is that the resulting contours have a gap where the tiles meet. The width of the gap is equal to the cell size I specify when creating the raster.

How can I eliminate this gap and get contours that are seamless across tiles?

Would it be better to start with the LIDAR data instead, or would I still encounter a similar problem?


texts: Setting the map extent - spatial reference for wkid 27700

texts:
I want to set the extent of my basemap on application start but when I add the following code it doesn't seem to like it and renders a blank div?

map = new Map("map", {
    logo: false,
    extent: new esri.geometry.Extent("xmin="-377930, "ymin"=470972, "xmax"=849893, "ymax"=1269406,  new SpatialReference({ wkid: 27700 }))
});

initLayer();
baseMapServiceLayer = new my.BaseMap();

I wasn't sure if I should replace esri.geometry.Extent with Extent but that makes no difference - any help would be good!


texts: Is it possible to have a single foreign key column referenced to two or more different tables?

texts:
In Spatialite (3.0.1), having a normal (non-spatial) table, I would like to relate it to two different spatial tables, with different geometries types (multipolygons and polylines), using the same column as a foreign key, and constraint it to their primary keys.

The general idea would be to, while inserting a new record in the non-spatial table, one can choose to relate it to a line or a polygon from two distinct spatial tables.

Example:

CREATE TABLE sp1_polygons
  poly_id TEXT NOT NULL PRIMARY KEY;

SELECT AddGeometryColumn(
  'sp1_polygons', 'geometry',
  27493, 'MULTIPOLYGON', 'XY');

CREATE TABLE sp2_lines
  poly_id TEXT NOT NULL PRIMARY KEY;

SELECT AddGeometryColumn(
  'sp2_lines', 'geometry',
  27493, 'MULTILINESTRING', 'XY');

CREATE TABLE t1 (
  t1_id INTEGER NOT NULL PRIMARY KEY,
  geo_id TEXT NOT NULL,
  CONSTRAINT fk_sp1
    FOREIGN KEY (geo_id)
    REFERENCES sp1 (sp1_id)),
  CONSTRAINT fk_sp2
    FOREIGN KEY (geo_id) -- Equal to before
    REFERENCES sp2 (sp2_id));

In my particular case, the non-spatial table represent management actions (harvesting, planting, cutting, fixing, etc...) in a farm, while the spatial tables represent homogeneus plots (polygons) to be managed (farmland, forest, natural vegetation), or homogeneus lines  (water streams, natural fences, etc...).


texts:
I'm trying to merge several shapefile layers.

Each of the layers has three attributes:

CGRSName: Name of the grid cell (same for each shapefile)
CGRSType: explains the type of grid (not important for this case)
Taxon: with the name of the species for each of the shapefile.

In example one shapefile attribute table would look like this:

I would like to merge all the layers in one, while also merging all the different taxa that repeats in a the same grid, i.e. 35TMH3 in one cell (like Taxon and one cell saying Aesculus hippocastanum,Rosmarinus officinalis).

I've tried using join attributes by localization and merging vectorial layers, but it always gives me two row with the same CGRSName but different Taxon.


texts:
I want to create a thematic map in arcmap using the fields NAME and NOTES for 3 categories of polygons:

where the name cell is null
where the name cell is filled
where the name cell is null, but it has the note "NUME EBRAICA", which means "name in hebrew"; this note appears only to the polygons where the name cell is empty.

Any suggestions, please? Everything I've tried didn't work.


texts: Satellite looking angle calculation (non zero latitude for subsatellite point)

texts: When calculating looking angles for geostationary satellites, it is assumed that because the satellites have an equitorial orbit, the latitude of the subsatellite point is 0. In reality, the latitude is not 0 (though really close to 0). If i just assume that the latitude is equal to 0, would that make my calculation inaccurate? If yes, how can we calculate angle of azimuth and elevation with non zero latitude for the subsatellite point.


texts: Bulk reverse geocoding with Geopy using built-in rate_limiter

texts:
I have a dataset with latitude and longitude coordinates, and I would like to get actual addresses for those. Since I'm using Geopy and Nominatim, I am limited to 1 request per second.

Although the examples I have seen to find a way around this involve using time.sleep(1) method, geopy has a RateLimiter class specifically for those purposes.

The only issue is that it doesn't seem to work on reverse geocoding for me, it keeps saying that the RateLimiter doesn't have an attribute reverse.

Am I doing something wrong?

from geopy.geocoders import Nominatim
from geopy.extra.rate_limiter import RateLimiter

geolocator = Nominatim(user_agent="application")

geocode = RateLimiter(geolocator, min_delay_seconds=1)

location = geocode.reverse((50.6539239, -120.3385242), language='en', exactly_one=True)

print location.raw

AttributeError: 'RateLimiter' object has no attribute 'reverse'


texts: OGRInfo On Shapefile = zero length delimited identifier error

texts:
I have an ESRI Shapefile that I can open in ArcMap and QGIS and show the attributes fine. However when I try to load to PostGIS using ogr2ogr I get an error saying
ERROR 1: ERROR:  zero-length delimited identifier at or near """"
LINE 1: ALTER TABLE "attribute" ADD COLUMN "" VARCHAR(5)

I used ogrinfo -al File.shp and it lists all 20 records fine but none of the attributes seem to come through.
I have opened the dbf file in Excel and its fine and have even parsed the dbf header using Python and it is exactly the same header output as a previous shapefile which correctly loads.
Why is ogr2ogr having an issue parsing the file?


texts:
I have a road network data where roads(edges) are specified with 3 values to identify the category of road(one way).

Road with value 0 can be routed in both direction.
Road with value 1 can be routed in the directed way.
Road with value -1 can be routed in opposite direction.

As seen in above picture. Blue edges can be routed in both direction, green edges should directed in the same direction how they are digitised and edges with red colors should be directed in opposite direction of the arrow.

I am trying to get route from one node to another using pgRouting on the above data.

pgRouting works fine if it is given as undirected graph with bellow query.

select seq, node, edge, pt.cost, agg_cost, st_astext(the_geom) as geom from pgr_dijkstra('select id as id, source, target,cost from ways_noded', 60414, 68634, directed := false) as pt join ways_noded rd on pt.edge = rd.id;

But I am facing problem with routing with one ways in directed graph. Above query with directed:=true considers the edges which are directed towards same direction. I am not able to direct with the edges which are mentioned as two ways(0) also.

OSM data also has a one_way field and this link  says 0 for both ways, 1 for one way and -1 for reverse way and they also have an additional value 2 for specifying 2 way roads.

Can anyone help me in routing through one way and reverse ways?
I found this link but didn't get much info on directed graph.


texts: How can I switch line direction from the Linux commandline?

texts:
I want to reverse the line direction for all lines in a shapefile from the Linux commandline (that is to say, using a script).  Currently I'm using the SwapVectorDirection QGIS plugin which works well but has to be run manually.  Ideally I'd like a solution using ogr since I have it installed on my system.

This question is similar but I'm looking specifically for something I can include in a Linux script.


texts:
I have a shapefile that has been drawn to encompass certain points in another layer.  I can use the intersection feature to gather all the data from the points layer, for points that lie in this polygon I have drawn.  However, I need to be able to export this to CSV.

How can I write a Python script that will pull all the points data that lie within this polygon.

points layer = ("xyz")
polygon layer = ("abc")

I am looking to have all the data from the points layer that are within my polygon layer.


texts:
I have a classified raster with two classes and I'm using GDAL python API to convert the raster to polygon using this function
gdal.Fpolygonize(raster, raster, dataset, 0, [], callback = None)

The problem is that it takes a lot of time (about 30 minutes) to convert the raster. I then tried to convert it using ArcGIS raster to polygon tool which is faster (about 7 minutes).
Is there a way to speed up the GDAL polygonize function or is there an alternate open source tool to accomplish this in a faster way?


texts: Reduce query time when using both MySQL and PostgreSQL with PostGIS

texts:
I am doing a tracking query. The device is sending its location data, timestamp to store in MySQL Database. I want to show the query result of every point with its nearest POI. The POI database is stored in PostgreSQL database as PostGIS format. Its taking 8-9 seconds to return the query result. At first I am calling the longitude latitude values from MySQL Database then calling the nearest POI info for every point from PostgreSQL database using while loop. I want to know is there any way to reduce the query time? I am using the following query

SELECT  name, area, district, ST_Distance(geom::geography,ST_GeographyFromText('POINT(
$longitude $latitude
)')) as distance
FROM poi
ORDER BY distance
Limit 1;

I need to get a JSON result to visualize on map.


texts: Reclassifying RasterStack by function and overlaying in R

texts:
I've got a RasterStack composed of a four rasters. See below:

# Reproducible examples
kiwi <- raster(xmn = -180, xmx = 180, ymn = -56, ymx = 84, res = 1)
kiwi[] <- sample(2:35, ncell(kiwi), replace = T)

banana <- raster(xmn = -180, xmx = 179.9999, ymn = -56.00081, ymx = 83.99917,
res = 1)
banana[] <- sample(2:100, ncell(banana), replace = T)

apple <- raster(xmn = -180, xmx = 180, ymn = -56, ymx = 84, res = 1)
apple[] <- sample(4:100, ncell(apple), replace = T)

mango <- raster(xmn = -180, xmx = 180, ymn = -56, ymx = 84, res = 1)
mango[] <- sample(1:90, ncell(mango), replace = T)

# Raster stack
stacked<-stack(kiwi, banana, apple, mango)

I want to reclassify the RasterStack into four categories (tropical, mediterranean, boreal and savannah) based on some conditions (see below), and then combine into one raster and plot all four categories on a single map.

# Example conditions
 kiwi < 10 & banana > 20 & mango > 0.5 == "tropical"
 mango > 15 & apple < 10 == "boreal"

I think I need to create a function and apply it to overlay but I don't know how to write this function. Does anyone know how to do this?


texts: How to reduce the size of shapefiles when features are deleted?

texts:
This tiny little script deletes all small features in all layers. This works fine but the size of .shp files are not reduced. Please tell me why...

canvas = iface.mapCanvas()
expr = QgsExpression("$length<100")
for l in canvas.layers():
    l.startEditing()
    it = l.getFeatures(QgsFeatureRequest(expr))
    for ft in it:
        l.deleteFeature(ft.id())
    l.commitChanges()