texts: In QGIS Atlas how to get the index and the number of features grouped/filtered by an Attribute with expression?

texts:
Using a QGIS Expression in the Layout Editor, how can one return the count of features matching a particular attribute value for an atlas feature, for example, "AREAS_Area".
In this particular example, the valid range of values for "Areas_Area" is 1-25.
I would like to return the count of features matching the "Areas_Area" attributes for 1, 2, 3, etc.
Here is what I have currently, which provides me with an index (current Area number out of total Area numbers), as well as the page count (current feature number out of total number of features).
concat('Area: ',"Areas_AREA",' / ',maximum("Areas_AREA"),' | ','Segment: ',@atlas_featurenumber,' / ',@atlas_totalfeatures)

Image preview:

So given the above information and the result I have in the preview, my expectation would be that it would instead read something more like:
Area: 1 / 25 | Segment: 1 / 13


texts: Aligning multiple images with different spatial resolutions

texts: I am using 2 different satellite images that have different spatial resolutions. One of them has 1 degree and other has 0.25 degree resolution. I want to subtract one image from another. So I should resample 0.25 to 1 degree resolution. When I do that, it is resampled to 1 degree but not overlay with other image. It is shifting below. Am I doing something wrong ? The both data have same projections and datum.


texts: GeoServer with SQL Table having lat/long

texts:
I have a table in SQL Server with geometry column.

I save geometry like : geometry::Point(@Lat, @Long, 4326)and this results the geometry for the given lat/long.

Now when I publish this layer in GeoServer with Declared SRS :EPSG:4326 and download the KML. It does place the point at right place on Google Earth and the layer does not render at right place in Open Layers as well.

What is wrong?


texts: Advanced Digitizing Tools Crashing in QGIS 2.8.1 x64 on Ubuntu 14.04

texts:
I'm hoping someone here can tell me what I'm doing wrong... Whenever I use the Advanced Digitizing tools in QGIS, I crash after about 5 clicks. I'm even doing it without OTF transformations on.

Any idea what's going on or if this is a known bug?


texts:
I have imported some shape files into SQL Server 2008. Now i want to use OpenLayers with my MVC 4 application to render those shape files from SQL Server Database.

Any help in this regard?


texts: Which points in projected raster corresponds to certain latitude,longitude location?

texts:
What I have is a SRTM 3 arc-second raster file in WGS84 geographic coordinate system (here is that file). I then project it to the Azimuthal Equidistant projection with GDAL:

gdalwarp -s_srs EPSG:4326 -t_srs "+proj=aeqd +lat_ts=40.81266 +lon_0=14.414252" -r near -of GTiff C:/vesuvius2_wgs84.tif C:/vesuvius2_aeqd.tif

I would like to know, how can I determine which point in the newly projected raster (C:/vesuvius2_aeqd.tif) corresponds to certain latitude,longitude location?

For example I would like to check where does the location of 14.4, 40.8 lay in the "C:/vesuvius2_aeqd.tif" file:

And for example it lays at -1051, -450 projected coordinates.

How can I calculate this?

EDIT:

Is it possible to know the exact pixel which corresponds to the required geographic coordinate location? For example, the 14.4, 40.8 geographic coordinate, corresponds to the pixel: 200x210 in the projected raster.
Is that possible?


texts: QGIS 3.2 Grass processing algorithms not working

texts:
All the Grass processing algorithms I have tried on windows 10 QGIS 3.2.1 give me the error below. Anyone else experiencing this? Any solutions?

Processing algorithm…
Algorithm 'v.buffer' starting…
Input parameters:
{ '-c' : False, '-s' : False, '-t' : False, 'GRASS_MIN_AREA_PARAMETER' : 0.0001, 'GRASS_OUTPUT_TYPE_PARAMETER' : 0, 'GRASS_REGION_PARAMETER' : None, 'GRASS_SNAP_TOLERANCE_PARAMETER' : -1, 'GRASS_VECTOR_DSCO' : '', 'GRASS_VECTOR_LCO' : '', 'angle' : 0, 'cats' : '', 'column' : None, 'distance' : 100, 'input' : 'G:/My Drive/GIS_DATA/Project Data/MMTC/pond.shp', 'minordistance' : None, 'output' : 'C:/Users/BAS/AppData/Local/Temp/processing_b45361d56a76487ba46a6c3186e7fd8f/5b2a3d4064364147ad44386aabf002d9/output.shp', 'scale' : 1, 'tolerance' : 0.01, 'type' : [0,1,4], 'where' : '' }

g.proj -c proj4="+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
v.in.ogr min_area=0.0001 snap=-1.0 input="G:\My Drive\GIS_DATA\Project Data\MMT\pond.shp" output="vector_5b78e50bc752f3" --overwrite -o
g.region n=5622558.817551329 s=5622558.817551329 e=-7832392.192254409 w=-7832392.192254409
v.buffer input=vector_5b78e50bc752f3 type="point,line,area" distance=100 angle=0 scale=1 tolerance=0.01 output=output472320f11531410082718db00660058f --overwrite
v.out.ogr -c type="auto" input="output472320f11531410082718db00660058f" output="C:\Users\BAS\AppData\Local\Temp\processing_b45361d56a76487ba46a6c3186e7fd8f\5b2a3d4064364147ad44386aabf002d9\output.shp" format="ESRI_Shapefile" --overwrite

'C:\PROGRA~1\QGIS' is not recognized as an internal or external command, operable program or batch file.
'C:\PROGRA~1\QGIS' is not recognized as an internal or external command, operable program or batch file.

File "C:\PROGRA~1\QGIS 3.2\apps\Python36\lib\site.py", line 177 file=sys.stderr)
SyntaxError: invalid syntax

Press any key to continue . . .

'C:\PROGRA~1\QGIS' is not recognized as an internal or external command, operable program or batch file.
'C:\PROGRA~1\QGIS' is not recognized as an internal or external command, operable program or batch file.

File "C:\PROGRA~1\QGIS 3.2\apps\Python36\lib\site.py", line 177 file=sys.stderr)
SyntaxError: invalid syntax

Press any key to continue . . .

Execution completed in 0.69 seconds
Results:
{'output': <QgsProcessingOutputLayerDefinition {'sink':C:/Users/BAS/AppData/Local/Temp/processing_b45361d56a76487ba46a6c3186e7fd8f/5b2a3d4064364147ad44386aabf002d9/output.shp, 'createOptions': {'fileEncoding': 'System'}}>}

Loading resulting layers
The following layers were not correctly generated.

C:/Users/BAS/AppData/Local/Temp/processing_b45361d56a76487ba46a6c3186e7fd8f/5b2a3d4064364147ad44386aabf002d9/output.shp

You can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.


texts:
Question: why does a raster object plot over an object of class "SpatialPolygonsDataFrame" but incorrectly over an object of class "sf"?

I am working with this raster. I'd like to overlay it on a map of Russia.

I'd like to use the graticule function in plot_sf, so I need to plot an sf object:

library(maptools)
library(sf)
library(raster)
# define CRS
crs.laea <- CRS("+proj=laea +lat_0=56 +lon_0=100
               +x_0=0 +y_0=0 +ellps=WGS84
             +datum=WGS84 +units=m +no_defs")
# load in raster
    agb<-raster(".../GEOCARBON_Forest_AGB_Map/GEOCARBON_AGB_Map.tif")
    #crop to desired area
    boreal<-crop(agb, extent(0, 180, 30, 80))
    russia_la<-projectRaster(boreal,crs=crs.laea)
# load in background map
    data("wrld_simpl", package = "maptools")
    wm <- crop(wrld_simpl, extent(0, 180, 30, 80))
    wm <- spTransform(wm, CRSobj = crs.laea)
    sf_worldmap<-st_as_sf(wm)

I now make a nice background plot:

plot(sf_worldmap[1], axes=T,key.pos=NULL,
     border="darkgrey", col="lightgrey",
     xlim=c(-3660660,st_bbox(sf_worldmap[1])[c(3)]-3000000),
     ylim=c(-1000000,ymax(russia_la)-1000000),
     graticule=st_crs("+proj=longlat +datum=WGS84"),
     lon= seq(20, 140, by = 20),
     lat=seq(40,60,by=10))

Then I go to add the raster, and it does not plot properly:

plot(russia_la, axes=F,key.pos=NULL, add=TRUE,
 xlim=c(-3660660,st_bbox(inline[1])[c(3)]-3000000),
 ylim=c(-1000000,ymax(russia_la)-1000000))

This is the resulting plot. (A similar image results from plot(russia_la,add=TRUE), but is not shown here for brevity.)

The same problem does not occur if I plot the raster over the SpatialPolygonsDataFrame, but as I noted earlier I want to use the graticule feature and that is only available with sf objects.

plot(wm, border="darkgrey", col="lightgrey", axes=F,
          xlim=c(-3660660,xmax(russia_la)),
          ylim=c(-1000000,ymax(russia_la)-1000000))
plot(russia_la,add=T)


texts: psexec remote file execution with ArcGIS Server

texts:
This is perhaps going a bit out of bounds, so i'm not confident I'll get an answer. The short of this is that I've created a python script which leverages psexec to start a bat file on a remote server. This works superb until I publish it. Fortunately I do know where the scripts and models get uploaded to the GIS server after publishing which can help a great deal with diagnostics. Prior to the psexec line in the script, I'm copying a file from one location to another, and I can see that is successful because of the time stamp of the copied file. So, i'm fairly confident that the server is getting hung up on the psexec portion.

Here is what the published script looks like.

The ArcGIS account does have read/execute access to the psexec.exe folder.
When the GP service is run, it just hangs until it times out. So, since I need to define a question: has anyone used psexec in the arcgis server environment? This should be a lot easier than it is. Is there a way to see if the process is hanging on the remote end?


texts:
I am trying to change the appearance of the QGIS-Vertexmarker that is used when snapping is active and vertex/segment is found where you can snap to.

I think i found the code snippet where in the cpp files the appearance of this vertexmarker is set:

void QgsMapToolCapture::cadCanvasMoveEvent( QgsMapMouseEvent * e )
{
  QgsMapToolAdvancedDigitizing::cadCanvasMoveEvent( e );
  bool snapped = e->isSnapped();
  QgsPoint point = e->mapPoint();

  if ( !snapped )
  {
    delete mSnappingMarker;
    mSnappingMarker = nullptr;
  }
  else
  {
    if ( !mSnappingMarker )
    {
      mSnappingMarker = new QgsVertexMarker( mCanvas );
      mSnappingMarker->setIconType( QgsVertexMarker::ICON_CROSS );
      mSnappingMarker->setColor( Qt::magenta );
      mSnappingMarker->setPenWidth( 3 );
    }
    mSnappingMarker->setCenter( point );
  }

(https://github.com/qgis/QGIS/blob/bb9e276d591594821e2368f6ca077850baf1220b/src/gui/qgsmaptoolcapture.cpp)

In the PyQGIS-API I haven't found an entry for the QgsMapToolCapture so I have no idea how to change the Snapping Marker using PyQGIS.

Any idea how I can change the color, size and marker-style for the SnappingMarker?


texts:
I am having some issues creating a geoprocessing service within arcmap using reportlabs pdf toolkit. The script runs fine outside of arcmap, however when i try and run it in arcmap to publish as a service I get this error:

Traceback (most recent call last):

  File "Y:\Server\Reporting\test\test_report.py", line 19, in <module>
    createPDF(pdf_path, user_input)
  File "Y:\Server\Reporting\test\test_report.py", line 16, in createPDF
    c.save()
  File "C:\Python27\ArcGIS10.3\lib\site-packages\reportlab\pdfgen\canvas.py", line 1215, in save
    self._doc.SaveToFile(self._filename, self)
  File "C:\Python27\ArcGIS10.3\lib\site-packages\reportlab\pdfbase\pdfdoc.py", line 219, in SaveToFile
    data = self.GetPDFData(canvas)
  File "C:\Python27\ArcGIS10.3\lib\site-packages\reportlab\pdfbase\pdfdoc.py", line 244, in GetPDFData
    return self.format()
  File "C:\Python27\ArcGIS10.3\lib\site-packages\reportlab\pdfbase\pdfdoc.py", line 424, in format
    IOf = IO.format(self)
  File "C:\Python27\ArcGIS10.3\lib\site-packages\reportlab\pdfbase\pdfdoc.py", line 869, in format
    fcontent = format(self.content, document, toplevel=1)   # yes this is at top level
  File "C:\Python27\ArcGIS10.3\lib\site-packages\reportlab\pdfbase\pdfdoc.py", line 80, in format
    f = element.format(document)
  File "C:\Python27\ArcGIS10.3\lib\site-packages\reportlab\pdfbase\pdfdoc.py", line 1558, in format
    return PD.format(document)
  File "C:\Python27\ArcGIS10.3\lib\site-packages\reportlab\pdfbase\pdfdoc.py", line 676, in format
    L = [(format(PDFName(k),document)+b" "+format(dict[k],document)) for k in keys]
  File "C:\Python27\ArcGIS10.3\lib\site-packages\reportlab\pdfbase\pdfdoc.py", line 80, in format
    f = element.format(document)
  File "C:\Python27\ArcGIS10.3\lib\site-packages\reportlab\pdfbase\pdfdoc.py", line 607, in format
    s.decode('pdfdoc')
  File "C:\Python27\ArcGIS10.3\lib\site-packages\reportlab\pdfbase\rl_codecs.py", line 1047, in _rl_codecs
    if name.startswith(e): return RL_Codecs.__rl_codecs(e)
AttributeError: 'NoneType' object has no attribute '_RL_Codecs__rl_codecs'

My very simple test code:

from reportlab.pdfgen import canvas
import arcpy

base_path = r"\\hc-gis\GIS_Data\Server\Reporting\test"

pdf_path = base_path +  r"\output\Report.pdf"
arcpy.AddMessage(pdf_path)

user_input = 'Welcome to Reportlab!'#arcpy.GetParameterAsText(0)

arcpy.AddMessage("Generating Report...")

def createPDF(path, input):
    c = canvas.Canvas(path)
    c.drawString(100,750,input)
    c.save()

if __name__ == "__main__":
    createPDF(pdf_path, user_input)
    arcpy.SetParameterAsText(0, pdf_path)
    arcpy.AddMessage('Complete')

My goal here is to return a PDF as a result of the user input.


texts: Is it possible to label the start and end of a line in MapServer?

texts:
I would like to be able to have the upstream invert at the top of the main and the downstream at the bottom.

<---Down Label---------------------------------------------Up Label--->


texts: How to set QGIS action to open and search for features on Google Maps?

texts:
At the moment, I have an action that takes each term in an address field and conducts a search on EACH term in Google (ie. multiple tabs are opened for each term). The action code I have is the following:

"C:\Program Files (x86)\Mozilla Firefox\firefox.exe" https://www.google.co.uk/?gws_rd=ssl#q=[% "ADDRESS1" %]

I would like to ask: what is the correct code to include ALL terms in the address field and, if possible, conduct a search on Google Maps instead of a simple Google search?


texts: Seeking street/road shapefile with TF_TIME and FT_TIME attribute?

texts:
Doing a snow removal project within a county figuring out the fastest route for them to traverse on the roads to remove the snow from the sidewalk

I was able to find the roads/street shapefile for the county and using network analysis route analysis I was able to find the shortest distance for them to travel. The issue is the shapefile does not have a time or travel time attribute where I can know how long it takes for the travel the roads, because the shortest distance is not always the fastest time.

I looked on most government websites and even used the Tiger/Line shapefiles but they do not contain the time attributes.


texts: How to convert a single band to multi-band raster?

texts: I have paletted or unique values but they are not different bands. Now, I want to convert this band 1 (Pallette) to different bands based on the value (0-17). Can anyone help me in this regard?


texts: How to transpose coordinates in Mapbox (either manually or programmatically)?

texts:
So I have a route going from California to India, and it looks like this currently.

The arc is currently going to the left, but I want it to go to the right (over Africa). Can I add a value to the longitude coordinate of the point in India so that it will go to the right? Sorry if this question is worded weirdly, I'm not sure of the correct terminology here.


texts: Automatically output joined fields in input layer while using spatial join

texts:
For example, I have two layers:

A polygon layer city block with a ID field city block index.
A point layer parking lot with a field parking lot size.

I want to know the total parking lot size for each feature in the land use layer, so I use the spatial join function. Here's my workflow:

Use join attributes by location (summary) to create a new temporary joined layer who contains a field called parking lot size_sum.
Go to the property panel of the city block layer and join the parking lot size_sum field as a new field called joined layer_parking lot size_sum.
In the city block layer, create a new field called total parking lot size .
Copy values from joined layer_parking lot size_sum and paste them into total parking lot size.
Delete the temporary joined layer.

All I want to do is to add a new field into the original city block layer by summarizing the parking lot size field. But I have to go through all the steps above. Is there a simpler way to do that?


texts: Image Segmentation tiff files for Object Detection

texts:
I have tiff files generated fro 12 bands of sentinel-2a data. I want to label each pixel of the tiff files as belonging to copper mine or normal.
After labelling I want the image format to stay as it is, with only labels for each pixel.
Is there any free tool to do this?
I have been trying on AWS, but it doesn't seem to work on tiff files.
EDIT : It would be helpful if it can be done through scikit-learn or through Python libraries.


texts: Rasters lose quality when Georeferencing using ArcGIS for Desktop?

texts:
I've started to have a problem when georeferencing some rasters in ArcGIS 10.2.2 where the image appears in much lower quality with false colors and higher contrast. The images appearing in ArcGIS do not look like the original rasters at all.  These are usually JPG images, although I think I have also had the same problem with TIFFs in the past. The display and pyramid settings are both set to use bilinear interpolation.

Sometimes this problem fixes itself after I update the georeferencing.  If it does fix itself, it then sometimes goes back to the poor image quality after I zoom in or out.  I haven't been able to find any rhyme or reason as to which rasters have this problem and when.

Does anyone know what is happening with these rasters or how to fix this issue?

The raster is appearing like this:

Instead of this (not shown at the same angle or scale, but you get the idea):

EDIT:

Taking your comments into account, I attempted a few more things, but still hitting a wall with this.

Running raster calculator cleaned it up a bit but didn't fully resolve the problem.  This issue is also happening with rasters that are grayscale to begin with.  It's also more than a contrast issue because text becomes illegible--it almost seems like it's downsampling somehow. This sometimes happens as soon as I bring the raster into the MXD.

Looking at Layer Properties > Symbology > Show: RGB Composite > Stretch,  I set the Type to "minimum-maximum" and unchecked "apply gamma stretch," which seemed to help a fair amount, but it's still not looking entirely like the original image. These are also not setting I've ever turned on or off before, so I don't know how they might be changing automatically.

Did this for both an RBG and a grayscale image with the same results.  I'm going to be working with dozens of rasters, so I'm trying to avoid running a ton of processes on each one.

Is pyramid information stored anywhere else besides in the .ovr file?

So, in the following image, I altered those things above in the Layer properties, rebuilt the pyramids, and here is what it looks like.  Far more legible, but still not right:


texts: Trying to remove all data from intersections in FME

texts:
I have road data and my goal is to get rid of any data where intersections are.  That way when two road features meet they do not intersect.  I am using FME from safe to complete this.

I tried using a "areaoverarea" transformer and create a list from its parameters and then use a "listduplicateremover"to remove duplicate lists that contain features with the same name.  I thought this would get rid of the overlapped areas so intersections were clear.  So far it has not worked.  I was wondering if anyone had any other ideas.

sorry for the crude drawings.


texts: UniqueValueRenderer and ClassBreaksRenderer configuration

texts: So I can create one of these renderers based on parameters I let the user set. Some parameter combinations either create invalid conditions or no conditions at all - missing data, too many breaks or too few, and so on. This obviously unacceptable. I've considered doing what the desktop product does and display the list of breaks and let the user fiddle with the parameters until they do something useful. Has anyone done that? Is there a sample somewhere on how to do that? If not can someone offer advice on either doing that or suggest another way to ensure that a renderer "works?"


texts:
In ArcGIS the IProximityOperator.QueryNearestPoint method is available to get the nearest point.

Is there any similar function available in PyQGIS?


texts:
I'm having trouble running a loop to update the values in a particular field in attribute table by reading from a csv file. Any input is greatly appreciated. The following is the sample code:

import arcpy

featureclass = "C:\\xyz"
table1 = "C:\\pqr.csv"
alpha= list(csv.reader(open(table1,'r')))

crows=arcpy.UpdateCursor(featureclass) # guess here it should be like cursor( , )

for crow in crows:
    for row in alpha:
        if crow.STNAME==row[3] and crow.ZIP==row[6]:
            if int(crow.STNUM)>=int(row[13]) and int(crow.STNUM)<=int(row[14]):
                if row[4]==crow.STTYPE and row[15]=="B":
                   # fill a field in attribute table with a particular row[] value from csv
                else:
                    if int(crow.STNUM)%2==0 and row[4]==frow.STTYPE and row[15]=="E":
                  # fill the field in attribute table with a particular row[]value from csv


texts: Cannot draw two layers at once using Leaflet

texts:
I have two map layers that I want to show at once: One is the full data, the other is a filtered version of the data. The map will draw with each of these layers, but not with both at the same time. Any idea what is going on?

I brought in my data using:

var contour = $.getJSON("StateParks.json", function(data) {
    var myStyle = {
    "color": "#ff7800",
    "weight":2,
    "opacity": 0.65

    };

A. Full data here:

L.geoJson(data,{
     style: myStyle
     }).addTo(map);
     });

B. Filter of that data here:
L.geoJson(data,{
     filter: stateFilter
      }).addTo(map);
      });
      function stateFilter(feature){
         if (feature.properties.State === "Wyoming"){
            return true
         }
 }

In other words, A will draw, B will draw, but both A and B will not.

My StateParks.json has a format similar this:

{"type":"FeatureCollection","features":[
{"type":"Feature","geometry":{"type":"Polygon","coordinates":[[[x1,y1],[x2,y2],[x3,y3]]]]},"properties":{"OBJECTID":1,"ParkArea":855,"Name":"Park Name","Picnic":"Yes","State": "Wyoming", "OfficeAddress": "Address"}}


texts: Connect to certain timestamp for child version

texts:
In ArcMap you can connect either to a historical marker or to a specified timespan by changing the version within the subtab Historic. I am however not quite sure on which parent-version this new version is based on. I suppose it is the DEFAULT-version because we need historic data that is only available through the archive-tables of our SDE.

Is it possible to connect to a specific timestamp for a child version?


texts:
Is there a way to log data about changes provided by user to file? I want to log which user made a certain change, when this change occurred and the content of this change.

Use case:

User logged in my QGIS app.
User made some changes (e.g.: deleted some line)
Details about change was saved to file (e.g.: log.txt)


texts: Import horizontal and vertical measurements into GIS

texts:
I would like to import horizontal and vertical measurements into GIS. Basically, it is a recording of change in the stratigraphy of an archaeological site. I measured across the excavation from a starting point to get the horizontal measurement and then down from that point to the feature I wanted to record. Each point therefore has two measurements (horizontal and vertical). Since the recordings were recorded on a vertical wall at a site, I don't have a GPS co-ordinate for the starting point.

Does anyone know how to do this?


texts: SAGA Merge Vector Layers 'Add Source Information' not returning correctly

texts: Running "Processing -> Toolbox -> Saga -> Merge vector layers" in QGIS with 'Add source information' checked returns a merged file with an additional attribute "Source" but the values are a series of numbers (such as 1485259209.6737) rather than layer names. The number is the same for all the polygons in the same layer, but I don't know what this number refers to. I need the layer names or at least some way of identifying the source from that number. I have tried this with layers in .sqlite format and in .geojson format.


texts: Raster statistics change after re-projection -- which one is correct

texts:
I'm looking into some global DEMs (like ETOPO1 - http://www.ngdc.noaa.gov/mgg/global/global.html), and noticed that the statistics for the elevation values can change quite a bit after reprojecting the raster.

The original data comes in a LatLong pseudo-projection. Its stats are:

from GRASS 7.1 r.univar
min: -10803
max: 8333
mean: -1892.40422534
stddev: 2649.98339303

from gdalinfo -stats
min: -10803.000
max: 8333.000
mean: -1892.404
stddev: 2649.983

Now if I project the raster to a cylindrical equal-area projection:

GRASS
min: -10803
max: 8333
mean: -2382.28934158
stddev: 2508.93105538

gdalinfo
min: -10803.000
max: 8333.000
mean: -2382.289
stddev: 2508.931

So the projection process is causing this, although I'm using a nearest neighbor method, so the values shouldn't change.

I keep thinking on the elevation values as a big array of values, and that these values shouldn't change. Right?

I would expect some change with other interpolators, but not nearest neighbor. Or am I missing something really basic here?

In any case, this raises the question: which one is correct?
(ignore the fact that the mean value doesn't really mean much due the bimodal distribution of elevations on Earth)


texts: How do I remove username / password from an ArcGIS mxd?

texts:
When you create an mxd and use layers from a database connection where you are using "Database authentication" and clicked "Save username and password" and share the mxd with others, why are others then able to access the database with my account? Or, how do I removed the cached username / password / connection from my mxd? Or, why is my account getting locked out everyday?

This is a terrible problem, first and most importantly, for database security, and second, the account will continually be locked out until you track down everyone who has the mxd that was shared.

(Why will it be locked? A. Your database may restrict the number of connections that your account can make and with multiple people using the mxd with your account, your account will be logged in too many times. B. When your database password changes, others will be using the mxd with your account and your old password that tries to login, the login will fail (multiplied by however many layers in the mxd), and it will lock your account for failed attempts.)


texts: Smallworld > FME > PostGIS > QGIS, handling geometry types

texts:
I'm migrating data from Smallworld via FME/SpatialBiz Plugin to a PostgreSQL/PostGIS Environment and then trying to open the resulting tables with QGIS 2.12.1. Problem arises, because Smallworld supports multiple geometries per record and for areas it does not really differ between curvepolygons, multipolygons, donuts,... and simple polygons as well.

Basically my FME process (which deaggregates multiple geometries and only keeps 'area geometries' [synonym for any kind of polygon in Smallworld] of a special name) runs without errors, but when opening the PostGIS tables in QGIS I get two entries for each table (see below), one of which with a defined spatial type (CurvePolygon for example), the other lets me select the geometry type manually.

Added:

SELECT * FROM public.klrg_statistik_bezirk WHERE ST_isValid(geom) <> true

Gives me an error: Unknown geometry type: 10 - CurvePolygon

No matter if I do so or not, I'm getting the following error (warning rsp.):

("Unsupported geometry type 'unknown type' in layer '...'")
... which pops up again and again when adding other layers, reopen the QGIS-project and so on.

Strangely, the data seems to be displayed correctly in QGIS, as you can see here:

At first sight, no problem comes up, when editing these PostGIS data QGIS, but I'm not quit shure if no issues arise in the future.

In the FME process i tried to fix this by exposing the format attribute 'postgis_type' and setting its value to 'postgis_multisurface' for example as lined out in the FME postgis writers documentation (I think multisurface is the geometry type of choice when polygon types of any kind come mixed up all together ...).

My questions are: How do I have to interpret this behaviour of QGIS? Might there be a way to fix this in FME at all, or is some postprocessing in PostGIS required? What kind? And are these the right questions at all?

Edit: Added screenshots of FME workspace

FME Workspace left part with SpatialBizPlugin reader / transformers. Note what the Deaggregators / GeometryPropertyExtractors / Testers do!

FME Workspace right part (directly behind AttributeKeeprs). Note that all features pass the GeometryValidators!!!

... which are configured like this:


texts:
I want to download my data from QGIS to an existing table in Postgresql. I tried to use DB Manager, but the field add to existing table doesn't work. Maybe somebody knows what the problem is?

Or is there another way to solve the problem?


texts: Iterate Table produces single output instead of all files? (ModelBuilder, ArcGIS 10.1)

texts:
I am using ArcGIS Desktop 10.1 and I would like to export attribute tables of my features stored in a geodatabase. I was following the post here to create the model: http://www.burdgis.com/blog/view/exporting-multiple-tables-from-a-geodatabase-using/. Only, instead of Table to Excel Tool, I used the Table to Table, as I want my output tables as .dbf. I have roughly 10 feature classes in my .gdb so I expect to get ~10 .DBF files back.

I have used Inline model variable substitution  t_%Name%.dbf (http://resources.arcgis.com/en/help/main/10.1/index.html#//002w0000005w000000) to keep the names and all of my output tables.

However, as a result I've got back the single file, instead of all iterated features.

What am I doing wrong?

I don't see the difference between the posted model and mine:


texts:
I am updating a QGIS2 plugin to QGIS3. The following code in PyQGIS2 works:

        association_symbol = QgsSymbolLayerV2Registry.instance().symbolLayerMetadata("ArrowLine").createSymbolLayer({
        'arrow_width': '0.4',
        'head_length': '2.0',
        'head_thickness' : '0.75',
        'color': 'red',
        'line_color': 'red',
        'is_curved': '0',
        'arrow_start_width': '.2'
    })
    self.road_object_path_association_layer.rendererV2().symbols()[0].appendSymbolLayer(association_symbol)

I made an attempt at correcting for PyQGIS3, but it does not work:

        association_symbol = QgsApplication.symbolLayerRegistry().symbolLayerMetadata("ArrowLine").createSymbolLayer({
        'arrow_width': '0.4',
        'head_length': '2.0',
        'head_thickness' : '0.75',
        'color': 'red',
        'line_color': 'red',
        'is_curved': '0',
        'arrow_start_width': '.2'
    })
    self.road_object_path_association_layer.renderer().symbol().symbolLayers()[0].appendSymbolLayer(association_symbol)

The error returned is as follows:

  AttributeError: 'QgsSimpleLineSymbolLayer' object has no attribute 'appendSymbolLayer'

The symbols() object required an argument.  I was not sure what argument to provide.  Further research let me to believe I could get the info from the symbolLayers()[0] object, but that appears to have been a red herring. Any assistance would be helpful.


texts: I created centroids for zip codes and now want to geocode a table of addresses to the zip code centroid. I created an address locator and chose zip code as type. However I am unable to geocode to the zip code centroid. The table of addresses consist of city, state, and zip.


texts:
I saw the latest update of QField 1.18 wherein users are allowed to draw/digitize in freehand lines/polygons geometries.
I was wondering how can it be enabled in QField?
We badly need it in our fieldwork/geotagging activities of natural resources.


texts: Calculating TurfJS planepoint values for a set of points / Get containing feature with TurfJS

texts:
Turf.planepoint() calculates the z value for a point that falls on the plane of a triangle.

I am struggling to make this work with a tin from turf.tin() and a set of points.

Is there a turf function that returns the feature which contains a point?

Testing every triangle against every point isn't feasible.

Turf features doesn't seem to have IDs. I could create an ID for every triangle, then turf.tag() the IDs to the points and finally select the triangle from the tin feature set by this id before evaluating the turf.planepoint().

Ideally turf.planepoint() would accept feature collections but I could easily get there if I only could get the containing feature for each point.


texts: How calculate areas between certain contour lines in QGIS?

texts:
How do I calculate the area between certain contour lines (e.g. 0-800, 800-1400, 1400-2100...) from a given shapefile. What I need is the area of the certain elevation levels.

shapefile is from vogis-data download (via FTP vogis.cvn.at - Gelaendemodell - Hoehenschichten - Vlbg), vector polylines, altitude is ranging from 200 - 3300 m, but not equidistances
lines to polygone seems not to work (too complex/large data?)
fieldcalculator $area did not (only lines?)
there is no DEM in raster format


texts: How to approximate the elevation of a shapefile in any pixel of a grid using a digital terrain model?

texts:
Introduction
I am working on some flood maps of a given region. As the data is confidential, I will try to stay as general as possible without any loss of relevant information. The truth is, my knowledge on GIS is limited.

Available databases / files for my analysis
Some data on a river and its bank:

A precise digital terrain model (.adf / ArcInfo binary grid format). The pixel size is 1m x 1m.
a polygon (one feature in an ESRI shapefile). This polygon contains the 2D information on the flood boundary (it delimits where a river will overflow). No altitude or Z elevation is contained in this file.

Even though that information is not contained in the polygon, the elevation at each point on the boundary polygon is basically the same.

The problem I am facing

  I would like to deduce the flood depth (the water height) observed in
  each cell of a grid (of the same resolution as the DTM). And then map
  it.

I am assuming this would be simple if we had the information on the elevation of the data contained in the shapefile. In this case, I would convert that shapefile to a raster and do a little subtraction. But I do not have it and I must deduce it.

I am comfortable with any coherent approximation. I must know this information at any pixel inside the polygon, not just on its boundary.

The tools I would like to use
I exclusively work with QGIS (in its latest version - Las Palmas). I do not have an ArcGIS licence.


texts: How to get features "not getting" intersected in arcobjects?

texts:
I have two feature classes. I'm intersecting both of them. I have to return features which does not intersect with the other.
Lets take an example.
I have two feature classes: villages and colleges. I'm intersecting both of them.
My query is "Select villages where population greater then 500 and not having colleges."
So villages "having population greater then" is working finely but then it spatially search for colleges in those villages.
what all I need is to select villages which does not have colleges.

public ArrayList Attributespatialquery(String StrTableA, String StrTableB, String Vcode, String GISColumn, String SDE)
    {

        string TableA, TableB;

        ESRI.ArcGIS.RuntimeManager.Bind(ESRI.ArcGIS.ProductCode.Desktop);
        AoInitialize initA = new AoInitialize();
        IWorkspaceFactory pworkspaceFactory = null;
        IWorkspace pworkspace = null;
        ESRI.ArcGIS.Geodatabase.IFeatureWorkspace pfeatureWorkspace = null;
        ESRI.ArcGIS.Geodatabase.IFeatureClass TableBfeatureClass = null;
        ESRI.ArcGIS.Geodatabase.IFeatureClass TableAfeatureClass = null;
        IFeatureCursor pFeatCusr = null;
        IFeatureCursor pFeatCusr1 = null;
        IFeature pFeature = null;
        IFeature pFeature1 = null;

        //System.Data.SqlClient.SqlConnection con = new System.Data.SqlClient.SqlConnection();
        //System.Data.DataSet dataset = new System.Data.DataSet();
        //// MyConnection objCon = new MyConnection(System.Configuration.ConfigurationManager.ConnectionStrings["DBConnectionGIS"].ToString());
        //con.ConnectionString = System.Configuration.ConfigurationManager.ConnectionStrings["DBConnectionGIS"].ToString();
        try
        {
            // initialize an ArcGIS Server license if the license is available
            if (initA.IsProductCodeAvailable(esriLicenseProductCode.esriLicenseProductCodeStandard) == esriLicenseStatus.esriLicenseAvailable)
                initA.Initialize(esriLicenseProductCode.esriLicenseProductCodeStandard);
            // throw an exception if the license is not available
            else
                throw new Exception("ESRI ArcGIS Server License is unavailable or has failed");

            Type factoryType = Type.GetTypeFromProgID("esriDataSourcesGDB.SdeWorkspaceFactory");
            pworkspaceFactory = (IWorkspaceFactory)Activator.CreateInstance(factoryType);
            pworkspace = pworkspaceFactory.OpenFromFile(SDE, 0);

            pfeatureWorkspace = (IFeatureWorkspace)pworkspace;

            TableA = "dbo." + StrTableA; // Explict Cast
            TableAfeatureClass = pfeatureWorkspace.OpenFeatureClass(TableA);
            TableB = "dbo." + StrTableB; // Explict Cast
            TableBfeatureClass = pfeatureWorkspace.OpenFeatureClass(TableB);

            IQueryFilter pQueryFilter = new QueryFilterClass();
            pQueryFilter.SubFields = "*";
            pQueryFilter.WhereClause = GISColumn + " in (" + Vcode + ")";

            pFeatCusr = TableAfeatureClass.Search(pQueryFilter, false);

            while ((pFeature = pFeatCusr.NextFeature()) != null)

            {
                ISpatialFilter spatialFilter = new SpatialFilterClass();
                spatialFilter.Geometry = pFeature.Shape;
                spatialFilter.GeometryField = TableAfeatureClass.ShapeFieldName;
                spatialFilter.SpatialRel = esriSpatialRelEnum.esriSpatialRelIntersects;
                pFeatCusr1 = TableBfeatureClass.Search(spatialFilter, false);

                while ((pFeature1 = pFeatCusr1.NextFeature()) != null)

                    if (!arrSpatial.Contains(pFeature1))
                    {
                        arrSpatial.Add(pFeature1.OID);
                        // WriteLog("values added");
                    }

            }


texts: How to draw Geodesic lines in Openlayers 2?

texts:
I want to display multiple lines between two coordinates. I will pass number of lines and color of the line, it has to display like below image.

please see the below image.this is my previous post


texts: Correct network analysis settings needed for shortest route problem

texts:
Why does the route analysis take a longer path to reach the point when i click the solve button, even though I can see a shorter path with my naked eye?

I'm trying to digitize again and make a new network dataset, but it's still same like before.

Am I missing any steps during creating the network?


texts: I am trying to download pgrouting but the download links are not working. I have tried many links but no link is working. can you please help me in this regard. will be thankful


texts: Function takes too long to run (over 72 hours)

texts:
I have an errorless code, which works for 3 days of daily data. however I wish to use it for 6000 days. It has not stopped running (we are on hour 75) There are two if statements inside a function with which I change a value if it is below a maximum or minimum, which I thought would be the problem, however when I removed those, the code is still running for more than a day (and going).
The piece of the code I am showing is the piece from which I expect the long runtime comes. It can however also be earlier on in the code.

var unweighted = function(img){var uni = img.reduceRegions({
reducer: ee.Reducer.sum().unweighted(),
collection: grid2,
scale:500}).map(function(feat){
  var id = ee.String(feat.get('system:index')).cat('-').cat(ee.String(img.get('system:index')));
  var date = ee.String(img.get('system:index'));
  var bandNames = img.bandNames().map(function(name){return ee.Number.parse(name)});
  var values = ee.Feature(feat).toDictionary(img.bandNames()).values();
  var array = ee.Array.cat([values, bandNames], 1);
  var min = array.reduce(ee.Reducer.min(2), [0], 1);
  var key = ee.String(min.get([0,1]).toInt());
  var minh = feat.get('min')
  var maxh = feat.get('max')
  if (key>maxh){var key = maxh}
  else if (minh>key){var key = minh}
  return ee.Feature(feat.setMulti(ee.Dictionary.fromLists(['Minimum','id'], [key, id])));
});
return uni;
};

Preferably I would still have a minimum and maximum check. What would be even better is that every band under the minimum and over the maximum are removed for each image seperatly. However at the moment I mostly have to tackle the runtime. Can anyone see why it runs for so long, and what can I do to streamline it?

Link to code


texts: Route between two coordinate points via line

texts:
I tried to search for an answer but couldn't find any so I decided to ask.

I would like to ask which would be the best way to get the shortest (and the only one) route between two points from a part of multiple LineStrings?

I have road on the map and it consist of multiple LineStrings. I have user clicked in two places on the map and those are on that road or near it. Now I would like to get the distance between these two points via the road.
Does anyone have suggestions how to do it? Is it better to do it in Java backend beforehand somehow or can it be done in Javascript frontend?
LineStrings are actually Feature-objects and are stored in MongoDb, so pgRouting won't be an option.
I had an idea of using the info from Features and their order on the road, but the problem occurs when I don't know is the nearest point of the clicked point before or after the clicked point in the coordinate array.


texts: Creating a Grid on a polygon and find number of points in each grid

texts:
How can I create a grid on top of a polygon and then insert the points onto the grid?

As far as i know, i am using How to create a regular polygon grid in PostGIS?  guide to created the grid. However, I do not know how to move the grid around the map.

Must I set the SRID?

I am planning to create the grid on the australia and new zealand map.


texts: How to extract coherence and covariance matrices of TerraSAR-X data?

texts:
I have obtained quad-polarized level 1b TerraSAR-X data from DLR. When I tried to read the data using PolSARpro, I understood that it is only capable of reading SSC format

But my data is in MGD (Multilook Ground Detected format) which causes PolSARpro produce the following error when trying to read the product file.

ERROR IN THE TERRASAR-X DATA MODE and/or LEVEL

The data is not problematic and it can be read and displayed by NEST easily. What can I do to extract T3 and C3 matrices in PolSARpro (envi) format like the following?


texts: QGIS 2.4 Print Composer Will Not Export Image

texts:
I am running QGIS 2.4 on a Macbook Pro (Mavericks) and am having problems with the print composer.

Recently, my composer window goes totally white when trying to export images or pdfs and nothing else happens.

Yesterday, I was able to export images, but today nothing.

I am still learning QGIS and am not sure if this is something I did or an issue with the print composer.


texts: How to UPDATE with LATERAL Nearest-Neighbour query?

texts:
Say I have table gps like:

CREATE TABLE gps(
   gps_id serial primary key,
   measured_timestamp timestamp,
   geom Geometry(Point,4326),
   stop_id int --Null
);

And I have a table of stops like:

CREATE TABLE stops(
   stop_id serial primary key,
   geom Geometry(Point,4326)
);

If I want to do an UPDATE on gps to find the nearest stop to each point, is there a way to use a LATERAL query?

I tried something like

UPDATE gps
SET stop_id = nearest.stop_id
FROM LATERAL(SELECT stop_id FROM stops ORDER BY stops.geom <-> gps.geom LIMIT 1) nearest

but that told me

ERROR:  invalid reference to FROM-clause entry for table "gps"
                                               ^
HINT:  There is an entry for table "gps", but it cannot be referenced from this part of the query.

So is the only way to do?

UPDATE gps
SET stop_id = nearest.stop_id
FROM (SELECT gps.gps_id, stop_id
      FROM gps
      LATERAL JOIN (SELECT stop_id FROM stops ORDER BY stops.geom <-> gps.geom LIMIT 1) stops) nearest
WHERE nearest.gps_id = gps.gps_id

This feels like joining the same table to itself, which wouldn't need to happen with a SELECT INTO


texts: Spatial query without spatial database module?

texts:
I have the following scenario:

Oracle database without the spatial module.
Application in J2EE with the following architecture:
-- Entity
-- Persistence (DAO / JPA + JDBC)
-- Business (BC)
-- Vision (MB + JSF)

I have an entity called EQUIPMENT must have the attribute COORD.
Persist this entity in the database is easy, I'll use JPA.

When I enter the system, it should use geolocation to determine the position of the equipment and check if its position is a certain radius from what is recorded in the database.

Which APIs should I use to do this type of query considering I do not have the spatial module of the Oracle?


texts:
I tried to change the only the shapefile datasource encoding through properties and also How to encode shapefiles from LATIN1 to UTF-8? but it doesn't change it.

Without saving the QGIS project can't I directly change the shapefile datasource encoding?

I also tried this code, it shows change encoding but as I open the shapefile again it show the same old system encoding:

for layer in QgsMapLayerRegistry.instance().mapLayers().values():
    layer.setProviderEncoding(u'UTF-8')
    layer.dataProvider().setEncoding(u'UTF-8')
    print layer.name(), layer.dataProvider().encoding()

(I am not saving QGIS project, just accessing and saving shapefile directly. When I save the project and do this task then it work correctly)


texts:
I want to create the Voronoi diagram from a set of point (see below). When I use the "Voronoi Polygons" tool from QGIS, there is an error :

  File
  "C:/PROGRA~1/QGIS3~1.0/apps/qgis/./python/plugins\processing\algs\qgis\VoronoiPolygons.py",
  line 128, in processAlgorithm self.tr('There were no polygons
  created.'))
  _core.QgsProcessingException: There were no polygons created.

What do you think is the reason for that?

I use QGIS 3.0, and the tool works for other sets of points (random, for example).


texts: Matching colours in already mosaiced single image?

texts:
I have an already mosaiced image for a large area. I need to match colours in this single image. I do not have each image separately. I just have this large single image downloaded from GLS. How can I match the colour of this image if I do not have its image in parts? I did a Histogram match but it does not work. Please provide me steps for any of the following:

Erdas 2014
ArcGIS
ENVI
QGIS


texts: google maps V3 : in javascript api, render route obtained with web api

texts:
I have to obtain route directions with the Google maps v3 WEB (http request) api (computations are made server-side, based on this data, without client interaction, so no JS available). Then, the routes are sorted out with internal logic, and the server (apache/PHP) send a web html/js/jQuery page to the client, which displays a page with a Google map (displayed with JS (ajax) api)

On this map, i would like to render the routes obtained earlier by the server with the web api.

So I tried to give the DirectionsResult to the google.maps.DirectionsRenderer::setDirections() but it does nothing, and firebug doesn't trace any error (and google JS if obfuscated, of course)

Yes, I first set the map with setMap(), and the directionsRenderer is initialized correctly, all that is OK !

It seems that de DirectionsResult is not the same when obtained from Web API and from JS API. Is there a way to "translate" it, so that it could be rendered ?

The alternative would be to do the route computations again on the client side with the ajax api, only to render them, but I think this is stupid to recompute directions as I already have them; so I don't want to do it this way unless it is the only way.

Thank you for your help !


texts: Postgis can't find the optimum way, seems to be always using Btree index rather than spatial index

texts:
I am using the Postgis ST_DWithin function like this:

SELECT *,
    ST_X(geometry) as longitude,
    ST_Y(geometry) as latitude,
    ST_Distance_Sphere(
        ST_SetSRID(
            ST_MakePoint(-118.478255,34.303046999999999),
            4326),
        geometry) / 1609.34 as DISTANCE
FROM  poi
WHERE ( category =10980400 )
    AND ST_DWithin(geometry::geography,
        (ST_SetSRID(
            ST_MakePoint(-118.478255,34.303046999999999),
            4326)::geography
        ),
        80467.2,
        false)
ORDER BY DISTANCE
LIMIT 10

I am using it for finding nearest points (like ATMs, restaurants) within a radius.

I have indexes available for geometry, geography and a btree index for category.
When I "explain analyze" the query I can see that the query only uses category index. Not sure if it happens always or is it a runtime decision.

The problem is that there are cases where the query performs better if I remove category index and the query uses spatial index.

Similarly, there are cases when I have all indexes and the query is better by using category index. But then the cases above degrade.

Why can't postgis decide the optimum way ?? I have run VACUUM ANALYZE on the table after having all the indexes.

SELECT version() : "PostgreSQL 9.3.14"
SELECT postgis_version() : "2.1 USE_GEOS=1 USE_PROJ=1 USE_STATS=1"

Here are the indexes I have:

CREATE INDEX world_poi_idx4
  ON public.poi
  USING gist
  ((geometry::geography));

CREATE INDEX world_poi_new_idx3
  ON public.poi
  USING gist
  (geometry);

 CREATE INDEX world_poi_idx1
  ON public.poi
  USING btree
  (category);

Here are the results of explain analyze (example of just one case, results are different in other cases, this is just to provide an idea):

When category index is present:
https://explain.depesz.com/s/iZWj

When category index is not present:
https://explain.depesz.com/s/0WH


texts: Create user defined file output location for geoprocessing on ArcGIS for Server 10.3

texts:
I am working on a geoprocessing service  that allows a user to interactively select shapes and then download a .dbf file of the selected shapes.  The tool was created in modelbuilder and uses the copy rows tool to create the dbf file.

The part I cannot figure out is how to allow the user to select the output folder, as I need to fill in the output path parameter in ArcMap.  When I upload the service definition from the results window, the file path gets hardcoded into the tool and is not exposed to the user.

Any leads on where to look? My google-fu has been an abysmal failure.


texts: How to add a column to table with sdelayer command?

texts:
I am working with an oracle11g database and I would like to add a new column to a feature class stored in this database. I would like to use the sdelayer command but can any one of you give  me a tip whether I am doing this correct

sdelayer -o add - l TABLE,COLUMN -e p  -i sde:oracle11g -s test -D oracle11g - u gis_test
-p gis@test

One more thing, I would like to specify the data type for my new column. Where should I put it in my command?


texts:
So I am trying to rotate map icons based on an "orientation" attribute field (which provides the bearing in degrees). I've tested this in Mapbox studio and it is working fine. When I copy paste the code from the studio window into my html/js project however, there is no rotation and I can't work out why. Code looks like below. Thanks for your help!
  map.addSource('points', {
    type: 'vector',
    url: 'mapbox://tileset'
    });

  map.addLayer({
    'id': 'lights',
    'type': 'symbol',
    'source': 'points',
    'source-layer': 'layername',
    'layout': {
      'visibility': 'visible',
      'icon-size': 1,
      'icon-image': 'beer-15',
      'icon-rotate': [
        "interpolate",
        ["linear"],
        ["get", "orientation"],
        0, 0,
        350, 350
      ],
    },
  });


texts: How to split lines into equal segments with QGIS?

texts:
I'm attempting to split a curved (wiggle) line into equal length segments. Have used Grass's v.split which works to a degree but does not provide the output I'm after.

v.split will extract segmented lengths based on a "maximum segment length" or "maximum number of vertices in segment". So, if I want to split the line into 10 meter segments and use a maximum segment length of 10 then all my segments come up just short, e.g. 9.98, or 9.6. for individual lines, wiggle or straight.

What I want instead is given for example a 55m line:
Start at one end and place a node/point every 10m along the line until it reaches the end. In the case above it would get split into 6 segments with the last segment being 5m long and all the others being 10m.

Is there  a better module?


texts: How to display distance in meters (km) and not in Decimal Degrees in QGIS 2.4?

texts:
I am a GIS beginner, and I have a question about how to get the distance between two points in meters instead of decimal degree.

In QGIS 2.4 I have a project with two point layers (in WGS84). From both layers I used I wanted to have the distance in meters between the two closest point between the layers. So what I did: Vector -> analyse tools -> distance Matrix, ticked the box "only nearest (k) target points: 1". This created a CSV file with all the distances, but in the wrong units.

So I looked it up on the internet and I found out two things:

I have to set my CRV in a projected coordinate system (like a UTM zone).
And I have to tick or untick the box with "on the fly" --> which I cannot find in QGIS2.4 Project --> project setting --> CRS.

I tried the first option. I tried to set the layers in UTM 29 ZONE IRENET95, EPSG: 2158 by "save as" --> add new vector layer. It still gives the distance in the wrong units. My project only contains point in southwest Ireland, so that is why I picked this projection. Or should I project it in a different projection?

I also tried to find a converter to convert the decimal degree outcome to meters, but I did not find any. I guess there is none. Am I right?

For me it seems like there must me an easy solution, but I did not manage to find one. I think I am doing something wrong with setting the layers in a different projection.


texts:
I am having a problem in QGIS with some layers failing to display on the map canvas.
I am trying to learn about QGIS by working through the Training Manual found on docs.qgis.org. I am running version 3.16.1 under Linux Mint 19.3. I installed the Flatpak version of QGIS  using the Linux Mint software manager. Several months ago I downloaded the training data that goes with the manual and worked on some of the exercises, and it all seemed to work well. But recently when I looked at the project the layer that was imported from a SpatiaLite database, called landuse, no longer displays on the canvas even though the layer is check-marked in the layers list. When I look at the attribute table for this layer and select Zoom to Feature on a record it zooms in but does not display the feature. When I select Flash Feature it does flash. Another layer, called roads, which is from a Geopackage database, also will no longer display on the canvas. Other layers which were imported as shapefiles are displayed. I tried downloading the training data again and connected the SpatiaLite and Geopackage files under new names. I imported the layers into a new project but the landuse layer (SpatiaLite) still does not display. The roads layer (Geopackage) does display but if I zoom in it disappears and will re-appear if I go to full zoom. Does anyone have any ideas to fix this?


texts:
I've been wrangling around in GFS grid based data (GRIB - mostly 0.5 degree resolution). I'm interested in plotting these by converting to geojson with ogr2ogr (which is quite easy). However, is there a best practice for interpolating the data so that the geojson is has smooth arcs? As you know, it's quite blocky. I know there's a ton of GIS software that will do this and spit out a graphic, but I'd rather have the geojson smoothed version.

Any ideas?


texts: Altering text elements in multiple mxds?

texts:
I am trying to make changes to titles on multiple mxds instead of having to open 1 by 1 but keep getting an error.

import arcpy.mapping as mapping, os

for root, dirs, files in os.walk(r"F:\Final Maps\New folder"): # GetPArameterAsTExt(0)
    for name in files:
        filename = os.path.join(root, name)
        if ".mxd" in filename:
            mxd = mapping.MapDocument(filename)
            for map in mxd:
                for elm in arcpy.mapping.ListLayoutElements(map, "TEXT_ELEMENT"):
                    if elm.text == "TEST 90666":
                        elm.text = "9067"

This is the error that I am getting:


texts: Convert and plot python basemap and netCDF file in UTM zone 45N

texts:
I have a netCDF file with precipitation output for one year. I have a problem with aligning the basemap and netCDF file, and plot them in the same coordinate system. The netCDF file covers the Himalaya region.

Here is a sample of my longitudes:

array([-1595000., -1585000., -1575000., -1565000., -1555000., -1545000.,
       -1535000., -1525000., -1515000., -1505000., -1495000., -1485000.,
len(lons) = 320

and here my latitudes:

array([ 4195000.,  4185000.,  4175000.,  4165000.,  4155000.,  4145000.,
        4135000.,  4125000.,  4115000.,  4105000.,  4095000.,  4085000.,
len(lats) = 190

Next I define a projection with Proj4 notation, to transform from UTM_45N to latlon

latlon = pyproj.Proj("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

Define projection using EPSG codes

UTM_45N = pyproj.Proj("+init=EPSG:32645")

Transform from UTM to latlon

tr_lons, tr_lats = pyproj.transform(UTM_45N, latlon, lons, lats)

--> RuntimeError: x,y and z must be same size

When converting using transform, the lons and lats should be the same length, so this does not work.

Alternatively, I could plot the netCDF file and  basemap in UTM zone 45N. However, basemap does not support EPSG code 32645

--> ValueError: 32645 is not a supported EPSG code

Plot data

m = Basemap(projection='tmerc', resolution='l',  height=3000000, width=3000000, lat_0=27.5, lon_0=82.5)

cs = m.pcolormesh(lons, lats, np.squeeze(prec_mean[0]))

Any ideas what's going wrong, or suggestions for a different approach?


texts: Misalligned coordinate system of overlapping shapefile

texts:
I am using QGIS and have 2 shapefiles related to a particular city: The city blocks and the census tracks. When I bring them both into QGIS and overlay them, it's clear that the boundaries should perfectly align, but for some reason the boundaries are off by a few meters: Where the thin black lines are the census track boundaries and the grey lines underneath are the city blocks. I think that the issue might be the coordinate system, but am fairly new to GIS and don't know how to manually set the

coordinate system such that the maps overlap perfectly.

Could someone provide a easy explanation of how to fix the misalignment?

EDIT: The current Spatial Reference System for the 2-layers (taken from the metadata) are:

City Blocks: +proj=tmerc +lat_0=0 +lon_0=-60 +k=0.999996 +x_0=400000 +y_0=5000000 +ellps=aust_SA +towgs84=-57,1,-41,0,0,0,0 +units=m +no_defs

Census Tracts: +proj=longlat +ellps=GRS80 +no_defs


texts:
I have a problem with GetFeatureInfo (I'm working with openlayers).
In my HTML file I have this:

Sieraków MAPA

And I want to show some info in div #info.
In javascript I have layer from geoserver:

var theurl = 'http://localhost:8080/geoserver/sierakow/wms';
    var jeziora = new ol.layer.Image({
          source: new ol.source.ImageWMS({
            url: theurl,
            params: {'LAYERS': 'sierakow:jeziora'},
            serverType: 'geoserver'
          })
    });
    map.addLayer(jeziora);

    map.on('singleclick', function(evt) {
            document.getElementById('info').innerHTML = '';
            var viewResolution = view.getResolution();
            var url = jeziora.getGetFeatureInfoUrl(
              evt.coordinate, viewResolution, 'EPSG:4326',
              {'INFO_FORMAT': 'text/html', 'FEATURE_COUNT': 50});
            if (url) {
              document.getElementById('info').innerHTML =
                  '';
            }
    });

To do this I followed this examples: LINK
I also try THIS ONE with parsing but result is the same, nothing happend. My div #info is still empty and I don't know why.

Console show me this errors:

Uncaught TypeError: jeziora.getGetFeatureInfoUrl is not a function

Uncaught DOMException: Failed to execute 'getImageData' on 'CanvasRenderingContext2D': The canvas has been tainted by cross-origin data.

This is url:

http://localhost:8080/geoserver/sierakow/wms?SERVICE=WMS
&VERSION=1.3.0
&REQUEST=GetFeatureInfo
&FORMAT=image%2Fpng
&TRANSPARENT=true
&QUERY_LAYERS=sierakow%3Ajeziora
&LAYERS=sierakow%3Ajeziora
&INFO_FORMAT=text%2Fhtml
&FEATURE_COUNT=50
&I=50
&J=50
&CRS=EPSG%3A4326
&STYLES
&WIDTH=101
&HEIGHT=101
&BBOX=6919511.668971639%2C1779979.5230483634%2C6921441.703935839%2C1781909.558012564

And this error shows:

Refused to display 'http://localhost:8080/geoserver/sierakow/wms?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetFeatureInfo&FORMAT=image%2Fpng&TRANSPARENT=true&QUERY_LAYERS=sierakow%3Ajeziora&LAYERS=sierakow%3Ajeziora&INFO_FORMAT=text%2Fplain&FEATURE_COUNT=50&I=50&J=50&CRS=EPSG%3A4326&STYLES&WIDTH=101&HEIGHT=101&BBOX=6918701.436880034%2C1786965.867550262%2C6920631.4718442345%2C1788895.9025144626' in a frame because it set 'X-Frame-Options' to 'sameorigin'.

I will be greatful for every tips how to do that :)


texts: Error #1064 from MySQL searching database for results near lat/long

texts:
I'm trying to run a SQL query to search for a term and a users location. I need help to search my database for locations near a user with the description of 'X'. The users location (lat, long) is drawn when they click 'find' on a form for 'X'.

Individually each query works fine. I can find all with a specific %search_term% using:

 SELECT * FROM  `table_name` WHERE  `column_name` LIKE '%search_term%'

and also find all the locations from the distant of the user using:

  SELECT *, (
  3959 * acos (
  cos ( radians($lat) )
  * cos( radians( lat ) )
  * cos( radians( lng ) - radians($long) )
  + sin ( radians($lat) )
  * sin( radians( lat ) )
  )
  ) AS distance
  FROM locations
  HAVING distance < 50
  ORDER BY distance
  LIMIT 0 , 20;

but together I receive an error #1064.

This is the query:

  SELECT * FROM  `table_name` WHERE  `column_name` LIKE '%search_term%' AND
  *, (
  3959 * acos (
  cos ( radians($lat) )
  * cos( radians( lat ) )
  * cos( radians( lng ) - radians($long) )
  + sin ( radians($lat) )
  * sin( radians( lat ) )
  )
  ) AS distance
  FROM locations
  HAVING distance < 50
  ORDER BY distance
  LIMIT 0 , 20;


texts: Output the scale value for each data driven page (DDP)

texts:
I'm using data driven pages (DDP) in ArcMap 10.3. The index layer is a polygon feature class and I would like to know the scale that each data driven page exports at. Ideally there would be a way to append an attribute to this layer that indicates the scale for each page (200+ pages). (I know DDP can be set on previously existing scale attribute value, but when its based on the various sizes of each polygon, I want to know what the numeric value is.

While that is my primary question, the ultimate goal is to be able to apply a different size label to corresponding features in each map. If I know the scale of each page, I can then adjust label size accordingly by page.


texts: Make a query to get extent from a set of features

texts:
I need to get the extent of a set of filtered features on a feature layer and zoom to that.

I tried the queryExtent method on the layer object, but with no success. The promise got back with an object that had the right count set on it, but there was no extent present in the result:

    var query = new Query();
    query.where = "application_id= '" + applicationId + "'";
    query.returnGeometry = true;
    layer.queryExtent(query).then(function(result){
        console.log('result', result)
        // result: { count: 2, _ssl: undefined }
        if(result.count > 0) {
            zoomTo(result.extent);
        }
    });

How can I zoom to the extent of a filtered set of features on the layer?


texts: GeoServer WFS-Transaction authentication in web browsers

texts:
I am developing a webGIS application with editable layers which uses transactional WFS through GeoServer to do CRUD operations in a PostGIS table. I was able to set up the authentication for wfs-t transactions in GeoServer and it works well: when i want to save the edits/new features a popup window appears with authentication request:

if I type the correct username&password the transaction makes the geometry changes.

The security issue is: this authentication will be kept until the whole browser window is open, it's not enough to close the browser tab. It won't ask again for the authentication until you close the whole browser and reopen it.

So if you don't close the browser anybody can edit your layer.

Do you have any idea how to "sign out" to forget the authentication keys, without closing the whole browser window?


texts: Automated georeferencing in ArcGIS Desktop?

texts:
I am attempting to fully automate the process of georeferencing a load of .png files - i have xml file containing the co-ordinates - but how can i use them?

Does anyone have any advice?


texts: Creating input node file for V.net.salesman in GRASS?

texts:
I have imported line vector layer and point vector layer(12 points) to grass using v.in.ogr.
Now i want to use v.in.salesman to find shortest path between these 12 points, but i couldnt input anything in "Input nodes" field.


texts: Vertices not showing in Edit session (ArcMap 10.5)

texts:
No matter what I am editing (buildings, roads, etc.) I can not get the vertices to show when I click "Edit Vertices" inside of an edit session. It works fine on my work laptop (ArcMap 10.3) but they just do not exist on my personal license(10.5). I am editing features blindly, clicking in random spots not seeing where new vertices are placed and hoping that when I double-click to finish the sketch that it looks somewhat right. The same goes for when I try to create features. I have no way to see where the vertices are being placed when I click!

The red dots are where I clicked (I did this with snip tool). No vertices showed up to guide me along the way when creating the feature. When I double-click to finish the sketch it draws but I had to do this blind. When I go to edit the vertices of that new line there is nothing that pops up!


texts: Using the value of a field in an S57 file with Mapserver 7

texts:
I'm relatively new to Mapserver and I'm having an issue with the line style for depth contours in my map.  Basically I just want the depth to be written along the line.

I can get text to appear along the line, but I can't for the life of me get the value of the VALDCO field of the DEPCNT to show up.  I think my basic problem is that I am not understanding how to access this information in my map file.  Here's what I'm working with:

LAYER
    NAME "DEPCNT_CA276241"
    TYPE LINE
    DUMP true
    TEMPLATE fooOnlyForWMSGetFeatureInfo
    EXTENT -66.1660611 44.6667528 -63.6660696 46.0000945
    DATA "DEPCNT"
    METADATA
      'ows_title' 'DEPCNT_CA276241'
      "gml_include_items" "VALDCO"
      "gml_VALDCO_type" "Real"
      "gml_VALDCO_precision" "2"
    END #metadata"
    STATUS ON
    CONNECTIONTYPE OGR
    CONNECTION "./s57/CA276241.000"
    TRANSPARENCY 100
    PROJECTION
      'proj=longlat'
      'datum=WGS84'
      'no_defs'
    END #projection

    LABELITEM 'VALDCO'
    CLASS
      NAME 'DEPCNT'
      STYLE
        WIDTH 0.91
        COLOR 0 170 255
      END #style
      LABEL
        #TEXT "[VALDCO]"
        SIZE 10
        COLOR 255 0 0
        ANGLE FOLLOW
        POSITION AUTO
        FORCE TRUE
        ANTIALIAS TRUE
      END #label
    END #class
END #layer DEPCNT_CA276241

My understanding from reading the documentation I can find is that the depth value (which is the VALDCO field of the feature) should be used as the label. Unfortunately the number is not showing up.

I also tried commenting out the LABELITEM and replacing it with the TEXT "[VALDCO]" line you can see in the code.  Some of the documentation suggested that was how to use a field.... all that does is write the literal text on the contour line.

So, I am missing something here.... suggestions?


texts: Rescale 300 polygons within one shapefile, but keeping the original center position of each polygon

texts:
I have a shapefile type Fishnet in ArcMap, containing 300 or so polygons in it.
Each polygon represent an area of interest containing plants. I have a model extracting the zonal stats of the shapefile for each of the polygons.
I want to reduce the size of each polygon where they are, without changing their positions.
If I use the scale option after selecting several of my polygons, it will move the little cross midway between all these selected polygons. And when the size change, the polygons will increase or decrease in size but will become farer or closer respectively.
Globally, when I select different polygon, I want each of them having the cross in them signaling the center, and resize around each of these cross, where they were originally.

Any ideas?

Thank you guys. I have tried to use "Buffer". However, it was not working on using negative values.
and when I tried using positive values to see if it can work in increasing the size of the polygons, it was only drafting a circle above my polygon, and bigger than it.
so it is not the tool that I need.

I remember having used one in the past, but impossible to remember which one. I was having 1600 circles on my "fishnet" grid. They were supposed to be positions of plants of interest in field.
And a tool allowed me to decrease by a ratio all the circles at once while keeping each of those in place (just a smaller radius that they have had before)


texts: How to plot a Start and Stop point while drawing a line in CartoDB?

texts:
My visualization plots a bunch of lines on a map using CartoDB using a Start coordinate and a Stop coordinate.

My data contains geojson formatted coordinates  (for each row) I want to show a "Start" and a "Stop" illustrating the start point and the end point.

{"type":"MultiLineString","coordinates":[[[-73.978,40.762],[-73.991,40.75]]]}

where the first coordinate is a start and the next a stop.


texts: OpenStreetMap boundary relation lost when exporting to spatialite with ogr2ogr?

texts:
I am currently working in map for a province of Papua New Guinea. In OpenStreetMap (OSM) I have recently updated and added new relations for boundaries of admin_level=5 and admin_level=6, which were not present before.

The relations I am analysing are the following: 5,259,918 and 5,259,919. You can check them here:

http://ra.osmsurround.org/analyzeRelation?relationId=5259918

http://ra.osmsurround.org/analyzeRelation?relationId=5259919

I extract the OSM data with http://extract.bbbike.org in PBF format, and then I export it with ogr2ogr using a command similar to the following

ogr2ogr -f "SQLite" -dsco SPATIALITE=YES png-full-test.db planet_xxx.osm.pbf

The resulting .db is later used in QGIS for creating the map. I use a custom osmconf.ini to extract the tags I need for my map.

The relations were created using both non-existing lines (borders where there is no other feature, created just for this purpose) and existing lines (such as rivers, which are part of the relation). The latter do not have specific tags for boundary=administrative or admin_level=6 as these are included in the relation (and according to my best understanding, there is no need to put them twice).

However, when I open the resulting file, only the lines which have the tags explicitly added will be shown. The ones such as rivers or roads, which are part of the relations but that they don't have neither boundary or admin_level tag, will not be listed in none of the resulting tables (other_relations, lines, points, multilinestrings or multipolygons). They are lost in the process. A quick and dirty solution would be to add the missing tags to the rivers / roads, but I am wondering if there is any other cleaner way to solve the issue.


texts: Dissolving polygon shapefile in Python without ArcPy?

texts:
I need to figure out a way to dissolve a polygon in Python using other modules besides ArcPy. I've been using the Python GDAL/OGR for merging shapefiles but I don't see anything to dissolve a polygon. I've been looking into Shapley, Fiona and Pyshp.

I also don't want to use qgis, as importing the qgis.core module gives my script problems.

Any ideas?

**EDIT I'm trying the Fiona approach but am having a hard time configuring it. I installed a .whl of it and it appears to install but when I run something it crashes Python saying

'This application has requested the Runtime to terminate it in an unusual way'

When I try to update the module in Pycharm, I get an error 'gdal version is not defined'


texts: Customized Annotation display with Data Driven Pages

texts:
I want to only display annotation from an annotation dataset that match the first three characters of the field in the shapefile used to drive a Data Driven Page (DDP) map export. (I'm running the map export through a Python script.)

Let's call the shapefile field that is used in the DDP FIELDNAME_A. Values in FIELDNAME_A have the format 360-2, 386-4, 516-2, etc.

I have an annotation dataset with a field called FIELDNAME_B which contains values of 360, 386, 516, etc.

Would it be to display annotation that matches the first three values in the DDP shapefile? For example:

DDP Shapefile FIELDNAME_A value: 360-2
Annotation data displayed from FIELDNAME_B: 360

Below are two screen captures of the attribute tables.

Data Driven Page shapefile FIELDNAME_A attribute table below

Data Driven Page shapefile FIELDNAME_B attribute table below

The out of the box option for a Definition query does not seem to work, since the values in FIELDNAME_A and FIELDNAME_B do not match.

Is there a query builder or Python work around?


texts: Automatically reclassify grid into 5 classes using Processing

texts:
I am using the graphical modeler of QGIS to rasterize a vector, and then reclassify it automatically into 5 values, without the intervention of the user.

I have reviewed the reclassify algorithm of SAGA but need user information to tabulate old value and new value.

How could you automatically calculate a grid in 5 values using the graphical modeler of QGIS?


texts: IS it possible to run GRASS commands in the QGIS graphical modeler?

texts: I have a geoalgorithm set up in the graphical modeler and everything runs except for the GRASS commands. It is really frustrating because when I execute these commands on their own they work just fine but as soon as I run them as part of a larger process they simply do not produce the desired results. I am using the mapcalc command to extract raster values below 10, for example.


texts: QGIS Python ImportError: No module named builtins

texts:
When I am executing a python script for doing a geoprocessing task using QGIS Desktop 3.8.3 (executed using a windows .bat file), I am getting an import error.

from builtins import zip
ImportError: No module named builtins

In C:\OSGeo4W64\apps\, I can find folders namely Python27 and Python37. I have installed Python 3.7 version before installing OSGeo4W64.

The .bat file (from which I am executing this .py file) is :

@echo off
call C:\OSGeo4W64\bin\o4w_env.bat
call C:\OSGeo4W64\bin\qt5_env.bat
call C:\OSGeo4W64\bin\py3_env.bat
REM Change OSGEO4W_ROOT to point to the base install folder
set QGIS_PREFIX_PATH=C:\OSGeo4W64\apps\qgis
REM Gdal Setup
set GDAL_DATA=C:\OSGeo4W64\share\gdal\
REM Python Setup
set PATH=C:\OSGeo4W64\bin;C:\OSGeo4W64\apps\qgis\bin;%PATH%
SET PYTHONHOME=C:\OSGeo4W64\apps\Python27
set PYTHONPATH=C:\OSGeo4W64\apps\qgis\python;%PYTHONPATH%

REM Launch python job
python C:/V/V.py
pause

The .py file starts like :

from osgeo import ogr, gdal
from gdalconst import *
from qgis.core import *
import qgis.utils, sys
from qgis.gui import *
from PyQt5.QtGui import *
from PyQt5.QtCore import Qt
from PyQt5.QtCore import QFileInfo
app = QApplication([])
QgsApplication.setPrefixPath("C:/OSGeo4W64/apps/qgis", True)
QgsApplication.initQgis()
sys.path.append('C:/OSGeo4W64/apps/qgis/python/plugins')
import processing
from processing.core.Processing import Processing
Processing.initialize()
from processing.tools import *

In Environment Variables, in User variables,

Path : C:\Users\sreeraj\AppData\Local\Programs\Python\Python37\Scripts\, C:\Users\sreeraj\AppData\Local\Programs\Python\Python37\
OSGEO4W_ROOT : C:\OSGeo4W64

And in System variables,

Path : C:\Users\sreeraj\AppData\Local\Programs\Python\Python37\Scripts\, C:\Users\sreeraj\AppData\Local\Programs\Python\Python37\

Inside python-core.bat in C:\OSGeo4W64\etc\ini\, I can find

SET PYTHONHOME=%OSGEO4W_ROOT%\apps\Python37
PATH %OSGEO4W_ROOT%\apps\Python37\Scripts;%PATH%

But I can't find any variables like PYTHONPATH or PYTHONHOME in Environment Variables (I don't know whether this information will be useful for solving this importerror or not.)

How can I solve this ImportError : No module named builtins ?


texts: QGIS create raster from CSV and find raster index

texts:
I have a CSV with about 2 Mio. entries with X and Y coordinates WGS84 and many attributes. The X and Y coordinate are the south west coordinate of 100 meter by 100 meter grid over a country. But the grid is not complete.
Then I have a second CSV with X and Y coordinates WGS84 with point data.
For each point of the second CSV I want to find the appropriate raster data so I have to identify in which raster the point lies...

I tried many approaches but I can't create a raster with the first CSV.

If I just create a layer with the CSV I can't rasterize since vector to raster expects a shapefile.
When I save my layer first as a shapefile the vector to raster option does something. But how can I specify that the coordinate is the south west point of a raster? And how can I specify the the raster should be 100 by 100 meters?


texts: ST_Transform not converting from utm to lat/lon

texts:
This line:

ST_AsText(ST_Transform(the_geom, 4326))

Returns the "coordinates" as utm coordinates. Although, I thought ST_Transform was supposed to convert the coordinates from utm to lat/lon. I found a tutorial where the behavior I need is shown in an example: http://postgis.net/workshops/postgis-intro/projection.html

Is there any reason why it would not be doing this?


texts: Overlapping polygons produce rendering artifacts a low zoom

texts:
I have duplicated this problem in both QGIS and mapnik. For this question, I will show the problem using QGIS.

I have imported heightmap images along a coastal area.

I then added an overlay layer with split water polygons from OSM. This renders perfectly fine at higher zoom levels:

Once I start to zoom out, lines appear between the water polygons. These lines make it into the final rendered image.

What can I do to remove these lines?


texts:
I am trying to merge shapefiles in QGIS it stops each time at 19%.
The program doesn't show me any error messages so I don't know what is going on.

Both shapefiles are 1.9 GB. Also I noticed that coordinate information in one of the shape files is not complete (missing data)

I am wondering, how merging of shapefiles work, on what criteria are they merged?

I suspect that there is missing data in those shapefiles.
Is there any strategy how to deal with this situation?


texts: Save a table in ascending order in order to avoid overlapping proportional circles (QGIS)

texts:
Is there a quick way to sort an attribute table in QGIS in ascending order permanently and save it? When you make a map with proportional circles in QGIS the circles overlaps between each other according to their sorting order. I want the smaller circle above the bigger which is why I think sorting can solve this issue.

Note: If I click on the header, the table will be sorted in ascending order using the clicked column. However, if I save it that way and I open it again nothing has changed about row order.


texts:
I'm trying to intercept the mouse event on a QgsMapCanvasItem.
I create a new class that extends the QgsMapCanvasItem, and override the mouseReleaseEvent (i put a qDebug inside).
As long as I tried a lot of time but the event is never fired.
There is something I'm missing?

I'm using Qt Designer and the compiled source of Qgis 2.99.

class PlatformMarker : public QgsMapCanvasItem
{
    ...
    private:
        virtual void mouseReleaseEvent(QGraphicsSceneMouseEvent *event) override;
};

In cpp:

void PlatformMarker::mouseReleaseEvent(QGraphicsSceneMouseEvent *event)
{
    qDebug() << "Custom item clicked.";
    QgsMapCanvasItem::mouseReleaseEvent(event);
}


texts: Python pyshp reads some field values as None

texts:
Reading the records of AL082017_pts.shp in a recent hurricane path file with pyshp results in a lot of None.

Here's the first record:

sf = shapefile.Reader('AL082017_pts.shp')
print sf.record(0)

Output:

['GENESIS014', None, None, '08', None, '0600', None, 'al', None, 'DB', None, None, None, None]

The correct field values, as can be seen in the dbf, are:

GENESIS014  2017080206  2017    08  2   0600    1012    al  8   DB  20  0   9.5 -13

It seems that the non-string fields are printed as None. How can I read these values correctly?

Here's the output of sf.fields:

[('DeletionFlag', 'C', 1, 0), ['STORMNAME', 'C', 50, 0], ['DTG', 'F', 50, 0], ['YEAR', 'F', 50, 0], ['MONTH', 'C', 50, 0], ['DAY', 'F', 50, 0], ['HHMM', 'C', 50, 0], ['MSLP', 'F', 50, 0], ['BASIN', 'C', 50, 0], ['STORMNUM', 'F', 50, 0], ['STORMTYPE', 'C', 50, 0], ['INTENSITY', 'F', 50, 0], ['SS', 'F', 50, 0], ['LAT', 'F', 50, 0], ['LON', 'F', 50, 0]]


texts:
i work with arcview 10.3 on 20 mxd's (each mxd has 25 layers) and i try to do select by location (with arcpy) and want to see a layers that are visible in the current map extent only. All Other layers that do not get in the current map extent i want to remove them from the data frame. I try this code:

import arcpy,os,sys,string
import arcpy.mapping
from arcpy import env

env.overwriteOutput = True
env.workspace = r"C:\Project"
for mxdname in arcpy.ListFiles('*.mxd'):
     print mxdname
     print '----Removed lyr----'
     mxd = arcpy.mapping.MapDocument(r"C:\Project\\" + mxdname)
     df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
     for lyr in arcpy.mapping.ListLayers(mxd, "" ,df):
         arcpy.MakeFeatureLayer_management(lyr.dataSource, "lyr")
         dfAsFeature = arcpy.Polygon(arcpy.Array([df.extent.lowerLeft,   df.extent.lowerRight, df.extent.upperRight, df.extent.upperLeft]),
                        df.spatialReference)
         arcpy.SelectLayerByLocation_management("lyr", "INTERSECT", dfAsFeature)

         if int(arcpy.GetCount_management("lyr").getOutput(0)) == 0:
              arcpy.mapping.RemoveLayer(df, lyr)
              print lyr

     mxd.save()
del mxd

i get an error when the layers are in Group layer:

>>>
Project.mxd
----Removed lyr----

Traceback (most recent call last):
  File "C:\yaron\shonot\software\gis\tools\YARON_SCRIPTS\SelectLayerByLocation2.py",     line 22, in <module>
    arcpy.MakeFeatureLayer_management(lyr.dataSource, "lyr")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\arcobjects\_base.py", line 78, in _get
    (attr_name, self.__class__.__name__))
NameError: The attribute 'dataSource' is not supported on this instance of Layer.
>>>


texts: What address locator to use for geocoding addresses in Excel?

texts:
I've been struggling to geocode the following addresses from an Excel table both in ArcMap and ArcPro:

As you can imagine, I'm trying to geocode the addresses in the "StreetAddress" field.
When I try to just simply add the table to my document and geocode it, I'm prompted to use an "address locator", of which none of the defaults work. So I believe I have to create an address locator. I'm uncertain of the style that I should use. I've tried streets, home address, and several others and I continue to get error messages.
Can anyone identify what style I need just by looking at the above example of the table? Is there something else that I'm not doing correctly? e.g. are commas invalid characters? Should I be using a different table type that isn't Excel? Are there address locators that I can just download and use this dataset for (instead of having to create one)? I can tell I'm missing something fundamental here.


texts: Calculating distance between points within grid polygon

texts:
I'm very new to ArcGIS, and I would like to calculate the average distance between points within a grid for every grid. Is there any idea on how to implement this?

My data look like this


texts: Problems with zooming and mouse events in Openlayers

texts:
I'm having trouble with everything mouse related for an Openlayers map when the map is placed inside a DIV with position fixed. The map works fine when you're at the top of the viewport of the browser, but when you're scrolling down the page, you can see that the zooming/click effects are shifted.  You can check here http://41.32.148.178:82/#!search=KFC%2F0%2F0%2F0%2F%2F to see exactly what I mean.

Thanks.


texts:
I am going through some code and trying to implement correct COM object releasing because currently in an error situation we cannot remove some generated files as the databases are kept locked.

I have already done this successfully in most places but this IFeatureLayer is causing some issues. Below is a code snippet in which we are modifying the datasource for each layer in map document. If I try to remove the database after this code has run, it fails as the database is locked. If I comment out the line featureLayer.FeatureClass = newFeatureClass; it successfully removes the database in same situation.

            IMap map = mapDocument.get_Map(0);
        for (int i = 0; i < map.LayerCount; i++)
        {
            ILayer layer = map.get_Layer(i);
            IFeatureLayer featureLayer = (IFeatureLayer)layer;
            IFeatureClass newFeatureClass = Utils.getFeatureClass(newDataSource, layer.Name);
            if (newFeatureClass == null)
            {
                string errorMessage = String.Format(
                    "Invalid data source. Feature class {0} was not found in the new data source.");
                Logger.Error(errorMessage);
                throw new ArgumentException(errorMessage);
            }
            IFeatureClass oldFeatureClass = featureLayer.FeatureClass;
            featureLayer.FeatureClass = newFeatureClass; // Not released correctly
            IMapAdmin2 mapAdmin = (IMapAdmin2)map;
            mapAdmin.FireChangeFeatureClass(oldFeatureClass, newFeatureClass);
            Utils.ReleaseObject(featureLayer);
            Utils.ReleaseObject(layer);
            Utils.ReleaseObject(oldFeatureClass);
            Utils.ReleaseObject(newFeatureClass);
        }
        Utils.ReleaseObject(map);

I have tried to release the mapAdmin, in which case it throws an exception as the map will get released as well. Setting mapAdmin to null didn't help. Also tried to release the featureLayer.FeatureClass but it didn't help either.


texts: SAGA GIS Kriging in QGIS calculating wrong values

texts:
I have a point shapefile concerning rainfall for the whole country. The points represent the annual rainfall in mm and range from 675-3562. I have approximately 800 points.

I am attempting to use Ordinary Kriging in QGIS to create a smooth raster surface.

I have used the Saga GIS ordinary kriging tool whith the following paramaters altered from the default:

        Block Size: 1000
        Cellsize: 1000

This has successfully created a smooth surface, however the range in values does not match my point data. The Ordinary kriging output has values ranging from 798-2273 - Why is this?


texts: Draw an envelope for only selected features in many/Batch layers in QGIS

texts:
I have about 1000 layers, I did select by attributes Batch files and select by one step features I need in all 1000 layers. Now I need to draw the envelope for those selected features in all 1000 layers in one step.
There is an option in a single layer process for draw envelope in selected features only by , but in the Batch files, I can not find this option, what can I do?


texts: Measuring tool and field calculator giving wrong area in QGIS 3.10.4

texts:
I've been consistently having issues with the QGIS measuring tool in 3.10.4. When I measure an area or line (regardless of the units) it gives a wrong value when compared to the measuring tool in Google Earth. This also applies when using the $area in field calculator. I am guessing this is a CRS issue, possibly related to the base map. I am using EPSG 27700 (British National Grid) and all my layers are also set to this.
The base map I am using is a Google aerial from QuickMapServices plugin, which I have manually set to convert to EPSG 27700 as otherwise it shows up somewhere completely else. If anyone has any ideas that would be really helpful, as I've tried to fix myself to no avail.


texts: How to create generic building footprints based on attribute information?

texts:
I'm trying to create a tool that generically creates building footprints on a given parcel according to certain attribute information of the parcel like building density, max. height, utilisation, etc.

Is there a way to generically create polygons in ModelBuilder with specified dimensions?


texts: List ArcGIS geodatabase feature classes on Python console

texts:
I am trying to create a new script with a hard-coded geodatabase filename, that lists the feature classes it contains to the console.
I know there are several ways to do this, but I decided to to use ListFeatureClasses function.

I got the following sample code from ArcGIS website, but I do not know what changes I need to make in order for it to work.
I barely started studying Python, so I know very basic stuff, but I would like to get help with this.

import os
import arcpy

# Set the workspace for ListFeatureClasses
arcpy.env.workspace = "c:/base"

# Use the ListFeatureClasses function to return a list of
# shapefiles.
featureclasses = arcpy.ListFeatureClasses()

# Copy shapefiles to a file geodatabase
for fc in featureclasses:
 arcpy.CopyFeatures_management(
  fc, os.path.join("c:/base/output.gdb",
   os.path.splitext(fc)[0]))


texts: How to go about building an offline KML viewer of a specific (fixed) area of land?

texts:
I am building a desktop application that would be used to process SAR images and monitor the deformation of a specific (fixed) site. One of the modules is generating a KML file of the deformation map and displaying it over a Google Earth-kind of 3D basemap. But I don't want to use Google Earth. I want to build a completely offline software. So, for this module, I would be needing a Google-Earth kind of 3D basemap of that specific site stored in my system. Over that, I would be displaying my generated KML file.
Any idea how to get such a 3D basemap, and how to overlay a KML file over it?


texts: Rule of thumb for selecting ROIs for supervised classification

texts:
I'd like to perform a supervised classification using the Semi-Automatic Classification Plugin and a maximum likelihood classification algorithm.
I read about a rule of thumb which says "if training data are being extracted from n bands, then a minimum of >10n pixels of training data are collected for each class" (Jensen). I'm using a Landsat image with six bands, so I assume that I have to select at least 60 training pixel per class. Somehow, I'm confused because I also heard people saying that you need 10times training pixels as you have classes.

Is there a general rule of thumb for selecting ROIs? Do you count the training pixels by drawing the polygons around them?


texts: Adding new column/field into attribute table in ArcGIS Online?

texts:
Is it possible to add a new column in the attribute table using the ArcGIS Online version or do I have to return to my home computer and do this through QGIS?

A little background to help better understand exactly what I'm doing:
I am working with ArcGIS Online at work and have successfully uploaded my parcel data after using QGIS on my home computer to convert it from the geo database into a shape file. I'm being asked to list each parcel by the property owner. I figured the best way to do this would be to add a field in the attribute table named 'owner' and populate the attributes that way.


texts: Crowdsourcing spatial data with online clickable map?

texts:
Here in northern Tanzania we are experiencing a huge westward butterfly migration - originally those of us who saw it thought it was local but it turns out that it's very widespread. We're getting informal reports from all over, but would like to (quickly!) set up a way for people to report sightings:

A clickable map on a web page that lets people show the location of a sighting;
A form that collects some other basic data;
Logging of spatial points and the form data.

Ideally this could be done through a Facebook page as social media will spread this like ... well, like a butterfly migration in a high wind ;).

Any ideas?

Here's the Google Map showing sightings.


texts: Add a GeoDjango PolygonField to a Leaflet map

texts:
Is it possible to add a GeoDjango PolygonField to a Leaflet map?
I've tried to do that by adding the js code below to my template, but I get a TypeError: t is null instead:

<script type="text/javascript">
    function map_init(map, options) {
        // zoom to point
        map.setView([51.9923, -2.1580], 12);

        // get polygon
        var polygon = "{{ obj.geometry.geojson }}";
        var polygon = L.polygon(polygon).addTo(map);
    }
</script>

In the code above, obj is my Django object and the PolygonField is called geometry. So I thought maybe I should use the coords property instead. But the same error persisted.

I can add a polygon normally if it's hardcoded like in the official tutorial:

var polygon = L.polygon([
    [51.509, -0.08],
    [51.503, -0.06],
    [51.51, -0.047]
]).addTo(mymap);)

Edit:

Strangely enough, if I turn the tuple string of {{ obj.geometry.coords }} into a list string and then json.parse it the code works, although there's a problem with the reversed (lat, lon) since this area which is supposed to appear in the UK, now shows somewhere in the Indian Ocean near Kenya:

<script type="text/javascript">
    function map_init(map, options) {
        // parse geojson
        var polygon_str = "{{ obj.geometry.coords }}";
        polygon_str = polygon_str.replace(/\(/g, '[').replace(/\)/g, ']')
        var polygon = JSON.parse(polygon_str);

        // zoom to point
        //map.setView([51.9923, -2.1580], 12);
        map.setView([-2.1580, 51.9923], 12);

        // get polygon
        L.polygon(polygon).addTo(map);
    }
</script>


texts: Distances and areas in Google projection

texts:
I have a mapping website used to capture polygons using OpenLayers in Google projection (900913) and storing them in SQL server 2008 as geometry. The SRID on the SQL geometry column is set to 900913. I don't reproject from the captured data, it's just sent to the server as WKT format from the features in OpenLayers and inserted into a geometry column using geometry::STGeomFromText().

When I call STArea() on the geometry I get hugely inaccurate results. For example, we have a gridded ordnance survey tile service that projects to 900913. If I trace a polygon as a 1km grid square from this background image and save it, the coordinates of the square in 900913 come out as around 1660m apart instead of 1000m and the total area is 2.7km2 instead of 1km2. This is baffling. The units of the Google projection are supposed to be in metres but metres of what?

If I import the geometry into Quantum GIS, I can calculate the area using the measure area tool. If I have the 'Ellipsoidal' checkbox on, I get the right calculation (1km2 in example above) and with ellipsoidal off, I get the inaccurate result (2.7km2).

There's a big gap in my knowledge here, but it looks like, despite the fact Google is a projected coordinate system, you can't use geometry methods for distance and area. You have to take into account the ellipsoid/sphere to get anything accurate. Is this true? Are there any good resources to find out more about this (I've seen the discussion on Alistair Aitchison's blog but there's quite a bit of assumed knowledge that I don't have towards the end)

The polygon capture on the OpenLayers map needs to be in 900913 so I can use Bing background mapping etc. So should I re-project before saving to sql server? If so, what would be appropriate (the extent of the mapping is UK which is both OSGB and OSI grids)?

In short, how should polygons captured in Google projection be stored if you want accurate area and distance calculations?