texts: Measuring flood mitigation effectiveness in river at different time using ArcGIS Desktop?

texts:
I am trying to measure the effectiveness of flood mitigation in a river using ArcGIS 10.2. Data available to me are the left and right bank polyline of 4 years , settlement data and point data where the mitigation is applied. We have applied river mitigation measures at the end of first year so that its effect is seen at the end of 2nd, 3rd and 4th year.(Like the shift of river, change in the area of river and effect on the settlements near the banks of river).

Is there any way I can measure the effectiveness of the mitigation applied using ArcGIS Desktop tools?


texts:
Can we include an iframe in tooltip in Marker?

import folium

var=open("yolo.html", "r")

var2 = var.read()

map = folium.Map(location=[19.179180,72.972580],zoom_start=6, tiles="Mapbox Control Room")

fg = folium.FeatureGroup(name="first layer")

iframe = folium.IFrame(html=var2, width=500, height=300)

popup = folium.Tooltip(iframe)

fg.add_child(folium.Marker(location=[19.279180,72.972680], tooltip=popup))

map.add_child(fg)

map.save("yolo1.html")


texts: Avoiding schema lock cannot be acquired because of existing lock on File Geodatabase?

texts:
I am having a problem with a Schema Lock in ArcMap for Desktop Advanced 10.2. The license that I have is Concurrent use for two computers. On Computer 1, I created an MXD Map File and a File Geodatabase for my map in the ArcCatalog Window. I have everything saved on the a shared folder in the E Drive over the Network.

I close ArcMap on this computer then open up ArcMap on Computer 2.
I open the MXD Map File and then open the ArcCatalog Window to add more Feature Classes to my created File Geodatabase.
When I try to expand the File Geodatabase on Computer 2 an error comes up stating "A schema lock cannot be acquired because of an existing lock." and the GeoDatabase will not expand for editing feature classes.

Are there any solutions to this problem?


texts: Need a Benchmark GPS Trajectory data-set?

texts:
I'm looking for a benchmark GPS data set, available for research purposes for free. I've found GeoLife GPS Trajectories dataset from Microsoft Research but I find it a little incomplete.

What I need is the GPS activity data of a person, as in (latitude, longitude, date) tuples, tracked for at least several months, preferably continuously. I would also like the recordings to be non-sparse; at most 1 minute between each record.

I'd really appreciate it if you can point me towards such a reliable data-set.


texts: Shows Netherlands region on map on every zoom level not other region in Open Street Maps

texts:
I want to show/display only Netherland region on Openstreet map, dont want to show other region on maps. Same thing apply for zoom in or zoom out also.
I have tried with applying bounds but does not work,when user decrease the zoom level it shows the whole world map.

Can you please help me or guide me in that.
Waiting for your response.


texts: GeoTools TileServer implementation for Google Maps

texts:
GeoTools TileServer has implementations for OSM and Bing.

Is there one for Google Maps?


texts:
Is there a way to reload a previous transfer attribute definition that I made?

In case my source and target schemes are the same.

Here's an example:

I've already defined field mapping between Makleb 4601 to water line, and now I want to do the same field mapping to Makleb 4603 and water Line

Is there a way to save the field mapping and reload it?


texts:
I'm trying to create a "heatmap" style visualization from a number of geopoints and values associated with them. I have been following along with this blog post.

I've got a CSV file of points like the following:

POINT_Y,POINT_X,Value
-41.157785,174.8343481,5
-41.1575385,174.8347111,4
-41.1574106,174.8347659,9
-41.1572644,174.8348964,2
-41.1570742,174.8351378,7
-41.1569067,174.835225,1

I've converted them into a GeoTIFF file with interpolation using gdal_grid. This was the command I ran:

gdal_grid -a invdist:power=3.0:smoothing=1.0:radius1=90000:radius2=90000:max_points=0 -a_srs WGS84 -of GTiff -ot Float64 -l test test.vrt test.tiff

That seems to work and generates a GeoTIFF file that looks right (I can paste the gdalinfo output, in case that helps)

I'm wondering how I style the points to get something similar to what is shown at the bottom of that blog post? I'm using TileMill to style my maps, and when I add the GeoTIFF to my project as a layer I get a white box that's in the correct place:

But of course I can't see any of the individual points on the grid. I'm fairly sure I'm missing an important step to colourize the points, I'm just not sure what it is!

I don't have any Carto styles applied to the layer in TileMill.

Thanks in advance!

Nik


texts: How to convert a shapefile polygon to raster in C#

texts:
Does anybody know a way to convert a polygon (vector) in a shapefile to raster data WITHOUT using ArcGIS Explorer Desktop? Is there a way to do this 'progmatically' in C# (WPF with ArcGIS runtime library)?

Any ideas? Suggestions?


texts: Move CartoDB-Map to show full InfoWindow

texts:
In CartoDB, if I customize the styles of an InfoWindow, the map moves not enough to always display the full InfoWindow.

For example: If I position the div.cartodb-popup to left:26px and bottom:16px to get it (the bottom left corner) open exactly where the user clicked, larger InfoWindows will not be displayed fully inside the map (Screenshot left), require the user to drag the map to fully see the infowindow (Screenshot right).

How can I tell CartoDB to move the map X pixels more to the left or right and top or bottom, so the user can see the full InfoWindow (Screenshot right)?


texts: Why do existing map services run slow after updating to ArcGIS 10.2.1 for Server?

texts:
Has anyone had any issues with existing map services after updating to ArcGIS Server 10.2.1?  After updating they are running slow and causing my ArcGIS Javascript app to crash.  All of our GIS services are on 10.2.1 now (Desktop, Server, SDE).

Any ideas?  I can provide more information if necessary.


texts: How can I check if two attributes match between a point and the area within which the point is located? in QGIS

texts: I would like to verify a geocoding process. I have used googles API and have geocoded about 35.000 adresses within Austria. By inspecting the locations it is obvious that a few lie outside of Austria (5 total) and that around another 30 are not coded in the right counties. This I could confirm visually, but I´m not sure how much I can trust the rest of the results and can´t check 35000 points visually.  Therefore I would like to see if the adresses lie in the right local communities. I have a point layer containing the street adresses and the respective community and a vector layer containing all the administrative boundaries. So I would like to do a sort of spatial query by location AND attribute, to see if the local community contained in my point layer matches the name of the community that is contained in my administrative boundary layer to root out obvious mistakes. Does anybody have an idea how I could do this?


texts: Changing layout output DPI and dimensions using QGIS Python

texts:
I would like to use Python to change the DPI of my map image files, I am not sure where in my code I need to add this line.
from qgis.core import QgsProject
from qgis.core import QgsPrintLayout
from qgis.core import QgsLayoutExporter
from qgis.PyQt.QtCore import QCoreApplication
from qgis.core import (QgsProcessing,
                       QgsFeatureSink,
                       QgsProcessingException,
                       QgsProcessingAlgorithm,
                       QgsProcessingParameterFeatureSource,
                       QgsProcessingParameterFeatureSink,
                       QgsLayoutSize
                       )
from qgis import processing

project = QgsProject.instance()
manager = project.layoutManager()
layout = QgsPrintLayout(project)
layoutName = "PrintLayout"

layouts_list = manager.printLayouts()
for layout in layouts_list:
    if layout.name() == layoutName:
        manager.removeLayout(layout)

layout = QgsPrintLayout(project)
layout.initializeDefaults()
layout.setName(layoutName)
manager.addLayout(layout)

from datetime import datetime

today = datetime.today().strftime('%Y-%m-%d')

manager = QgsProject.instance().layoutManager()

countrylist = ['Africa', 'Asia', 'Caribbean', 'World map', 'Europe', 'Middle East', 'Oceania', 'South America', 'North America']

for i in countrylist:
    layout = manager.layoutByName(f'{i}')
    exporter = QgsLayoutExporter(layout)
    exporter.exportToImage(f'D:/Sanch/Work/Automate/{today} {i}.png', QgsLayoutExporter.ImageExportSettings())


texts: arcgisscripting module not found in Python 3.5 in Spyder IDE

texts:
I am running ArcGIS Pro 2.2 and I have cloned my conda environment from Python command prompt. When I activate this environment and attempt to run Python 3.5 from the terminal, I can load arcpy as the cloned environment runs the python located at 'C:\\Program Files\\ArcGIS\\Pro\\bin\\Python\\envs\\arcgispro-py3\\python.exe'. However, Spyder links to the python installed in the cloned environment and I get the following traceback:

File "<ipython-input-4-f7f77bfeaa96>", line 1, in <module>
    from arcpy.geoprocessing import gp

  File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\__init__.py", line 66, in <module>
    from arcpy.geoprocessing import gp

  File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\geoprocessing\__init__.py", line 14, in <module>
    from ._base import *

  File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\geoprocessing\_base.py", line 14, in <module>
    import arcgisscripting

ImportError: DLL load failed: The specified module could not be found.

I have followed the directions located here and here, and the paths mirror a previous installation of ArcGIS Pro 2.0 from which I can use Spyder and arcpy together.

Does anyone know how to resolve this?


texts: Identifying multiple polygon features that contain point features and adding attributes to table in ArcGIS Desktop?

texts:
I have a point layer and a polygon layer, both containing hundreds of records. Each point feature can be within multiple polygon features (max. 4), and each polygon feature can contain multiple point features.

I would like to

know which points are inside each polygon feature and based on that analysis,
have a new field in the point attribute table identifying which polygon it's within. In case of a point being within multiple polygons the point attribute table would have multiple new fields.

I have used Near-tool which would work perfectly if a point could be only within one polygon (distance to the nearest polygon is 0 and the tool identifies one of the polygons), but in my case it's missing the other polygons a point is within.


texts:
I'm trying to convert an equirectangular map into a sinusoidal interrupted projection (with 10 or 20 degree gores) with the intention of making a globe with the gores once printed out. However, I'm running into a problem with the lines of latitude - sinusoidal projection projects the latitudes as parallels, however this results in lines of latitude that are geometric rather than circular, particularly close to the poles (I hope I'm making sense).

Is there another projection, similar to sinusoidal interrupted that distorts the lines of latitude in such a way that when printed and attached to a globe would result in circular lines of latitude, as desired? And if so, how can I (easily) go about changing an equirectangular map into such a suitable projection?

(NB: as though my questions werent a giveaway as it was, I'm in no way a professional when it comes to GIS! I'm using photoshop to create the equirectangular map and G. Projector to convert it into Sinusoidal interrupted projection, and that's about the extent of my knowledge! Also, the map is not of earth so I can't just find one and use that (though I might be able to take a blank set of such gores and convert my existing map to it)


texts: Exporting USGS contour line shapefile to DWG or DXF?

texts:
How do I assign the Z value based on the Elevation attribute using ArcGIS 10.1 for Desktop?

I've been messing with this for an hour or two now and every time the exported DWG has ZERO for z value for all lines.


texts: Exporting raster values into ASCII file format?

texts: How do I get all pixel values of a time series NDVI image (raster) and export its values into an ASCII file?


texts: Determinating direction of profile graph

texts:
I'm creating profile graph for craters. Each of polylines was created same way:

bearing distance to line in two opposite direction from crater's center, for example with bearing angle 30° and 210°
merging together each pair of lines
dissolve each pair
interpolate them with DEM

Now I have 8 polylines through each crater center. Lets say, now I want to create profile graphs. But can I be sure that graph was created from "left" to "right", not from "right" to "left"? (on screenshot) I don't have time and resourses to manually check each graph and compare it with DTM. My created polylines only have X and Y coords of center of crater. I do know though value of bearings of each profile.

So, I need to know direction of created profile graphs. Maybe there is a basic algorithm for arcmap to create them, like "start from top left then to bottom right"?


texts: 95% contour from kernel density estimates

texts:
I am looking to calculate the 95% volume contour from a kernel density estimate to estimate an animal's home range.

Normally, I work in R, but I calculated the "kernel density estimate" in ArcGIS so that I could incorporate barriers to movement. To do this, I used the "kernel interpolation with barriers tool". The output of this tool is a raster, where the value for each grid cell is the estimated number of records per km2, given the number and distribution of the actual records it was calculated from. values where the animal is not expected to occur have a value of 0 of less.

My question is, given the output of this tool, how can a calculate the 95% contour? Ideally I would like a solution in R. I know how to import the raster in R, but I don't know of any functions that can accomplish what I want.

I tried the function getVolumeUD from the move package, but got an error:

Error in (function (classes, fdef, mtable)  :
unable to find an inherited method for function ‘getVolumeUD’ for signature ‘"RasterLayer"’


texts: Label toolbar not showing up even though it is checked in toolbar menu

texts: I have recently discovered that even though I chose to display the label toolbar through the toolbar menu that it is nowhere to be seen on the monitor. Has anyone encountered this or know a fix to it?


texts: Can I import vector data into a SQL Database with GDAL using a stored procedure?

texts:
I am trying to understand if I can use a GDAL utility like ogr2ogr to load data from a CSV file that contains Point features, into my PostGIS-enabled database.
I need to be able to pass the point features in a stored procedure that inserts data into multiple tables as a single transaction for each file. From my understanding this is not possible with ogr2ogr but maybe I miss something or there is another utility for my use case?


texts: Web service to get near by road on giving Latitude Longitude

texts:
Is there a web service available which on providing a Latitude and Longitude can give the nearest road with its type i.e. whether the road is a highway or district road?

Google reverse geo-coding just gives type route but no further classification.


texts:
I am using QGIS 3.4.11 on Windows 10 x64.

I have three large DWG files that are the output of a survey (each roughly 5MB).  No problem to import them into QGIS as geopackages using the built in DXF/DWG import at geopackage tool under the Project menu.  When I do that, the result is a Group of layers (probably 20 or so per geopackage).

Problem is that the layers are not properly georeferenced when I do that (ie they appear in the middle of the Indian Ocean).

Is there some easy way to either edit the extent property of the entire geopackage, specify the extent at import time, or to otherwise georeference the entire geopackage all at once?

I have plenty of control points, but the georeferencing tool that I use for vector files doesn't seem to work for geopackages.


texts: Agent Analyst - Arc 10.5 Refresh has stopped working, config file issue?

texts:
I'm going through the Agent Analyst EsriPress exercises and cannot get the model to change the ArcGIS agent layer (positions) during a simulation run through the updateDisplay() function within Agent Analyst. I'm using ArcGis 10.5.

I've attempted the fixes described in Agent Analyst tool interface will not open on ArcGIS for Desktop?
but have not had any success getting past exercise 2f (moving agents). Tool will run once and then I get a popup saying "Refresh has stopped working." When I change the Refresh.exe.config file to 10.5.0.0, Refresh appears once. I tried changing it to 10.5.0.6491 (my current version), the refresh error flashes quickly.

Steps Executed:

Created Refresh.exe.config file with a text editor and saved it along with Refresh.exe (in C:\Repast 3\Agent Analyst\Refresh). The Refresh.exe.config contents are identical to the contents of the configuration file described here ArcGIS 10.1 using Repast's Agent Analyst Extension - refresh.exe does not work.  Replaced newVersion="10.1.0.0" with newVersion="10.5.0.0"
Edited the agent_analyst.bat and agent_analyst_run.bat files in the installation folder (C:\Repast 3\Agent Analyst) so that the folder paths contain Desktop10.5 in the folder path strings.
Edited Environment Variables to make sure that ARC_JAVA, ARC_JAVAW, and ARCOBJECTS_JAR have a fully specified paths entered for each variable.
In my Environment Variables, I moved ARC_JAVA (C:\Program Files (x86)\ArcGIS\Desktop10.5\java\jre\bin\java), ARC_JAVAW (C:\Program Files (x86)\ArcGIS\Desktop10.5\java\jre\bin\javaw), ARCOBJECTS_JAR (C:\Program Files (x86)\ArcGIS\Desktop10.5\java\lib\arcobjects.jar), AGENT_ANALYST (C:\Repast 3\Agent Analyst), and AGENT_ANALYST_VOL (C:) from my User Variables to my System Variables.


texts: SAR GRD preprocessing workflow necessity

texts:
I was working on Sentinel 1 SAR GRD (EW mode) imagery from Copernicus when I came across this workflow here.

Is it necessary for GRD products to pass through all steps of this workflow before they can be used?
For context, I'm looking into applying segmentation and classification algorithms on SAR imagery and was wondering if this workflow is necessary.
I'm new to SAR image processing.


texts: Clip a shapefile based on raster extent with high accuracy

texts:
I am trying to clip a shapefile based on the extents of a smaller raster. I find the extent of the raster using this code:

Extents = arcpy.sa.Raster(MyRaster).extent`

Then, I use these lines to make a polygon out of these extents:

pnt_array = arcpy.Array()
pnt_array.add(Extents2.lowerLeft)
pnt_array.add(Extents.lowerRight)
pnt_array.add(Extents.upperRight)
pnt_array.add(Extents.upperLeft)
poly = arcpy.Polygon(pnt_array)

Finally, I cut clip my shapefile using this code

arcpy.Clip_analysis(shp, poly, Shp_clip)

The problem is that I am using accuracy in two places. Please take a look:

Extents.XMax =                          -86.22649016383502
poly.extent.XMax =                      -86.22650146484375
arcpy.Describe(Shp_clip).extent.XMax=   -86.22650146517043

Is there anyway that I can get a clipped shapefile with the extent exactly as Extents?


texts: RTKLIB PPP algorithm against Single point poisitioning

texts:
I've been looking for a standalone single frequency carrier-phase based positioning method and I found out PPP (I know it's not the same). According to literature, the main issue of the carrier-phase based positioning is the resolution of the integer ambiguity and this is usually done by:

Differencing the observables from the user receiver (or rover) with the measurements from a reference receiver that is spatially close in order to remove common biases.
Removing the biases directly by either using a linear combination between observables, or estimating them or obtaining their values from an external source.

Now, I've recorded data with a single frequency ublox receiver and plotted it with RTKNAVI (real-time GUI-API from RTKLIB) as static PPP and single point positioning. The result is much better accuracy using the PPP method than with SPP one even if no reference or correction data is used.

How can this happen?

How is the integer ambiguity computed just using the raw-data from the ublox?

How can this be understood taking into account the previously mentioned points?


texts: Using schema other than public with PostGIS?

texts:
I've read (too late) that it's best to use another schema to store data tables in than the public schema when using PostGIS.

e.g. here:
Using a schema other than public in PostGIS?

Considering this is all quite old (e.g. the above accepted answer is from 2012), I'm wondering if this is still as valid as it was back then for current versions of PostgreSQL and PostGIS ?


texts: Overwrite existing map service in ArcGIS Server through ArcPy?

texts:
I am using ArcGIS Desktop and Server 10.1. I want to over write existing map service though ArcPy and going through this link  and getting an errorRuntime error
Traceback (most recent call last):
  File "<string>", line 22, in <module>
NameError: name 'draftPath' is not defined and I am searching for this at my end.

I just want to know is there any alternative  to overwrite map service automatically after 1 hour in ArcGIS Server 10.1 ?


texts: Get popup info in multiple layers from one click

texts:
I have multiple layers, say a site location as a marker and a town location as a polygon. I want to be able to click on the marker and push it's attributes to a dialog, and also send the polygon's name attribute to the same dialog.

So I need to get both layers information from the same click. The marker is easy and in it's onEachFeature I can run an on-click function, but how do I fire off the click for the second layer by only coordinates?

End result would be a dialog with a tab for each layer, to show it's information.


texts: How to perform accuracy assessment in eCognition?

texts: I have already done the classification using Nearest neighbourhood, and how do I continue to complete AC with TAA?


texts:
I have a shapefile with one point (starting point of a rockfall) and want to rasterize this to do further calculations.

When I rasterize it into a Zero-Raster-File (called "startgrd", which is DEM*0), it gets me a new raster with only one value for all cells (0).

Do you know what might be the problem?


texts:
I have a problem when I create a contour line from DEM, it's become sharp and not correct ( image attached ) .

DEM cell size = 5m
Contour interval = 5m


texts:
In ArcGIS and SAGA GIS there is the "identity" tool, which creates a vector overlay and shows area contributions of polygons in shapefile B to the polygons in shapefile A (in the attribute table).

Which tools would be appropriate for this in Grass GIS? I tried "v.overlay", but somehow the result is not equal.


texts: In OpenLayers, is there any way to convert polygons to tiles?

texts:
I have a OpenLayers map which by enabling a feature by user, a very huge table of polygons (lots of small areas, with information and metadata for each) will appear on the map, which rendering this big amount of data is very slow and zooming in and out is kind of miserable.

Do you know any solution to this problem?

I was thinking of converting these polygons to tiles, which of course if possible, will not show any information for areas. But when the user reaches a certain zoom level, map starts showing information for areas, in case that area is being clicked by user. What I mean is behaving like tiles before a certain zoom level, and after that, behaving like polygons.

I will appreciate any suggestion to the problem itself, or to my (maybe crazy) proposed idea.


texts: How to avoid "Handle bad layers" message when opening a QGIS project file?

texts:
Most of the time I have to reload all the layers when opening a QGIS project file (.qgs), which is pretty annoying.

I didn't change the name of the files and folders and I didn't change the folder structure.

What am I doing wrong?

I use QGIS 2.8 on OSX Yosemite 10.10.3


texts: How to remove temporary layer source files in QGIS

texts:
I am new to QGIS and python.

Titles says it all. I'm modifying a plugin that creates some temporary files (layers) during it's run and I need to loop it a couple times. The algorithm works fine for the first 2 loops then it crashes (usually with a "Bad file descriptor" error). I believe the fact that it's creating all those extra files is causing problems, because even if I run the unmodified version of this algorithm a couple of times the 3rd run always gives me an error.

In more detail what I'm doing. I am modifying the ConcaveHull algorithm so that instead of just taking all the features from a layer and making one big hull for all of them, the user has the option to group the features from input layer by their values in a chosen field and make a separate hull for each group.
I use a separate layer (assistant_layer) and copy features from a chosen group onto it, then I feed the layer into the code of unmodified ConcaveHull as the input layer (that part works right), it does its thing, at the end I simply remove all features from my assistant_layer and move onto the next group.

Suspect parts of the code:

#dissolve all Delaunay triangles
progress.setText('Dissolving Delaunay triangles ...')
dissolved = processing.runalg("qgis:dissolve", delaunay_layer,
                                              True, '', None)['OUTPUT']
dissolved_layer = processing.getObject(dissolved)
#save result
progress.setText('Saving data ...')
feat = QgsFeature()
dissolved_layer.getFeatures(QgsFeatureRequest().setFilterFid(0)).nextFeature(feat)
writer = self.getOutputFromName(
                        self.OUTPUT).getVectorWriter(layer.pendingFields().toList(),
                                                     QGis.WKBPolygon, layer.crs())
geom = feat.geometry()

On the 3rd iteration the code crashes on either the "dissolved_layer = processing.getObject(dissolved)" line with the "Bad file descriptor" error or (in rare cases) on "dissolved_layer.getFeatures(QgsFeatureRequest().setFilterFid(0)).nextFeature(feat)" with an error message "'NoneType' object has no attribute 'getFeatures'"

At the end of each loop I delete the writer using "del writer".


texts: Project location on a (great circle) path

texts:
I've been searching this SE site for quite a few hours now, and I'm still struggling to find a solution to my question. My goal is that given a way in OSM and my location (lat/lon coordinates), I want to find the closest location (lat/lon coordinates) on that way. The point can be anywhere on the way, not limited to the points used to define the way.

So I'm thinking of the following algorithm:

Separate Path into separate edges, each edge connecting only two points.
Select the closest edge.
Project my location onto that edge.

Now there are many questions about calculating the distance between a location and a path:

WGS point to WGS line segment (great circle) distance
Calculating the distance between a point and a virtual line of two lat/lngs
How to approximate point to segment distance on sphere?

Also a very similar question of which I cannot get the calculations right or verified:

Find the nearest point on a line

There is also some info from Dr. Math about that subject. However I cannot seem to find an algorithm to calculate the location in step 3. As I haven't touched (vector) algebra in quite a while, I don't quite understand the logic in those answers.

Can someone show an algorithm to do this? A solution in any reasonable programming language is fine with me.


texts: QGIS 3 - Variables in rule-based rendering?

texts:
I'll try and lay this out clearly, I have most of it covered, but there's a very specific thing I'm looking for in the rules - whether or not it can be done and if so, how?

Key things I have done and/or plan as outcomes;

QGIS 3.4.2
A dataset of just over 556000 rows - all points representing one building each
583 groupings for each by round which are unique, the name covers region, service, round number and day like so '(SH) Food 1 Mon'
These rounds have anywhere between a few hundred and several thousand points in each
The groupings have been first set as categorised and then Rule-Based so each round is split out from the others
They're to be exported via an atlas and show one round at a time (once I've bounding geometry to use for each round) and all will have the same symbols - with only one round showing each time, there's no demand to have different symbols

I'm looking to use an atlas due to the volume of images this will produce, doing it by hand would take an age.

Previously when doing similar work on far less categories I've used this code and adapted it for each category (this is using the round example from above)

if("Round Name" = @atlas_pagename, "Round Name" = '(SH) Food 1 Mon', Null)

That works perfectly once all categories have the right round name in before the Null and also the ELSE category has been deleted.

Now here's where I've become stuck - with 583, it takes an awfully long time to change each category by hand once I've moved from Categorised to Rule-Based.

So, the key question here is this; there a way of putting a variable in where I've been typing the round name in order to have QGIS do that part for me by reading and using each round name for each category?

I've looked around on here and not found anything apart from using some Python script which I'm not very handy with so thought I'd ask this first before trying to get my head around that.

I think that covers it and hopefully it's all clear!


texts: Using more than one graphics layer with a FeatureDataGrid

texts:
I am wanting to use the FeatureDataGrid to display data. The problem is I have more than one layer that contains data that needs to be displayed.

This is how you use the DataGrid for one layer:

private void MyMap_Loaded(object sender, RoutedEventArgs e)
        {
            MyDataGrid.Map = MyMap;
            MyDataGrid.GraphicsLayer = MyMap.Layers[1] as ESRI.ArcGIS.Client.GraphicsLayer;
            MyFeatureDataForm.FeatureLayer = MyMap.Layers[1] as ESRI.ArcGIS.Client.FeatureLayer;
        }

This is the XAML for my map:

<esri:Map.Layers>
                    <esri:ArcGISTiledMapServiceLayer ID="MyLayer" Url="http://services.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer"/>

                    <esri:FeatureLayer ID="Hydrant"
                                       Url="http://localhost:6080/arcgis/rest/services/CalvertCity_Test_Server/FeatureServer/0"
                                       Renderer="{StaticResource Renderer_Hydrant}"
                                       MouseLeftButtonUp="FeatureLayer_MouseLeftButtonUp"
                                       DisableClientCaching="True"
                                       OnDemandCacheSize="0"
                                       Mode="OnDemand"
                                       AutoSave="False"
                                       OutFields="*" />
                    <esri:FeatureLayer ID="WaterValves"
                                       Url="http://localhost:6080/arcgis/rest/services/CalvertCity_Test_Server/FeatureServer/1"
                                       Renderer="{StaticResource Renderer_WaterValves}"
                                       MouseLeftButtonUp="FeatureLayer_MouseLeftButtonUp"
                                       DisableClientCaching="True"
                                       OnDemandCacheSize="0"
                                       Mode="OnDemand"
                                       AutoSave="False"
                                       OutFields="*" />
                    <esri:GraphicsLayer ID="CandidateGraphicsLayer">
                        <esri:GraphicsLayer.MapTip>
                            <Grid>
                                <Rectangle Stroke="Gray"  RadiusX="10" RadiusY="10" Fill="#77FF0000" Margin="0,0,0,5" >
                                    <Rectangle.Effect>
                                        <DropShadowEffect/>
                                    </Rectangle.Effect>
                                </Rectangle>
                                <Rectangle Fill="#DDFFFFFF" Stroke="DarkGray" RadiusX="5" RadiusY="5" Margin="10,10,10,15" />
                                <StackPanel Orientation="Vertical" HorizontalAlignment="Center" Margin="30,20,30,30">
                                    <TextBlock Text="{Binding [Address]}" HorizontalAlignment="Left" Foreground="Black" />
                                </StackPanel>
                            </Grid>
                        </esri:GraphicsLayer.MapTip>
                    </esri:GraphicsLayer>
                </esri:Map.Layers>
            </esri:Map>

I thought about trying to pull the layer for the selected graphic, but I am not sure as to how to do that. If you need any more information or want me to explain further, let me know.


texts:
I'm trying to execute a 'dissolve' operation in PostGIS using the ST_Union command.

The input layer is admittedly quite large and complex. By 'large' I mean 57,771 features, with number of vertices ranging from 4 to 758,018 per feature, averaging around 86 vertices per feature. Only about 10 of the features have >10,000 vertices. By 'complex' I mean that the polygons have lots of holes, messy overlaps, islands, etc. and that the large polygons tend to have a bounding box that covers many of the smaller polygons, perhaps rending indexes less useful.

The problem is that the the query is extremely slow to the point of being unusable. I read Paul's 2009 post here that lead me to believe that my query should still be fairly fast. I'm using the following command; am I doing something blatantly wrong or inefficient?

SELECT  ST_Union(f.geom) as geom, column1,column2,column3
FROM "inputlayer" As f
GROUP BY column1,column2,column3

Edit: I am using:

POSTGIS="2.1.4 r12966" GEOS="3.3.3-CAPI-1.7.4" PROJ="Rel. 4.7.1, 23 September 2009" GDAL="GDAL 2.0.0dev, released 2014/04/16" LIBXML="2.7.8" LIBJSON="UNKNOWN" TOPOLOGY RASTER PostgreSQL 9.3.5 on x86_64-unknown-linux-gnu, compiled by gcc (Ubuntu/Linaro 4.6.3-1ubuntu5) 4.6.3, 64-bit

The machine I'm running the db server on is a virtual machine without a lot of power. I'll try the SET work_mem=50000 idea and see how things go!


texts:
I have to merge few shapefiles into one and create a heat map out of it in QGIS. I'm a first user of the QGIS, I tried to do it there: vector-management tools-merge files into one, but it kept to make all the time just one output file out of the layer I wanted to merge, so I didn't actually merge them, but create something like copies out of them ...

I wish to have something like this:

(this pic is from my professor) but I have to have one for another group of tourists,

so I did this:

then I input a shapefile from the folder that consists all the shapeless that need to be merged, the I browsed the output file and named it Merge:

an this is what popped up, and I have no idea how to go on now.


texts: ArcGIS Viewshed's show no visibility to south east

texts:
I'm using ArcGIS 10.2.2. I've been running some viewsheds using DTM files from EDINA Digimap and data that I've imported using .csv files. These files only contain the X Y details and a code for each site. The DTM files were 'stitched' together using the 'Mosaic to New Raster' tool in the ArcToolbox.

Out of nine that I've run, six viewshed's show no visibility to the south east.

I've checked the following.

I've run the viewshed again just using the one square and nothing else. That was fine.
I've run it using the mosaic to new raster on the NE, SE, SW and NW cells as one. This was fine.
I then ran it again using the mosaic to new raster but this time included the four cells below, so its not 74 NE, SE, SW and NW 75 NE, SE, SW and NW all stitched together. This showed the same problem with no view to the south east.

Has anyone had any issues like this before? It seems odd (to me anyway).


texts: Trying to build OSRM: Could NOT find OSMPBF

texts:
I am trying to build OSRM on a Ubuntu system (14.04 LTS, up-to-date, 64 bit), but I am getting a "Could NOT find OSMPBF" error.

First I followed the Ubuntu 14.04 instructions on the Building on Ubuntu page. This appears to have installed everything successfully.
Next I followed it with the General Build Instructions:

git clone https://github.com/Project-OSRM/osrm-backend.git
cd osrm-backend
mkdir -p build
cd build
cmake ..
make

This fails at the cmake line with:

-- Building on a 64 bit system
-- Configuring OSRM in release mode
-- LTO working
-- Boost version: 1.54.0
-- Found the following Boost libraries:
--   date_time
--   filesystem
--   iostreams
--   program_options
--   regex
--   system
--   thread
--   unit_test_framework
-- Found Intel TBB
-- Looking for Luabind...
-- Found Luabind: /usr/lib/libluabind.so
-- Looking for LuaJIT 5.2
-- Could NOT find LUAJIT (missing:  LUAJIT_LIBRARIES)
-- Performing Test LUABIND_WORKS
-- Performing Test LUABIND_WORKS - Success
-- Luabind/Lua5.2 combination working with /usr/lib/x86_64-linux-gnu/liblua5.2.so
-- Looking for STXXL...
-- Found STXXL: /usr/lib/libstxxl.so
-- OpenMP support found. Linking just in case for stxxl
-- Looking for OSMPBF...
CMake Error at /usr/share/cmake-2.8/Modules/FindPackageHandleStandardArgs.cmake:108 (message):
  Could NOT find OSMPBF (missing: OSMPBF_LIBRARY OSMPBF_INCLUDE_DIR)
Call Stack (most recent call first):
  /usr/share/cmake-2.8/Modules/FindPackageHandleStandardArgs.cmake:315 (_FPHSA_FAILURE_MESSAGE)
  cmake/FindOSMPBF.cmake:46 (FIND_PACKAGE_HANDLE_STANDARD_ARGS)
  CMakeLists.txt:283 (find_package)

-- Configuring incomplete, errors occurred!
See also "/home/richard/osrm/osrm-backend/build/CMakeFiles/CMakeOutput.log".
See also "/home/richard/osrm/osrm-backend/build/CMakeFiles/CMakeError.log".

I can post the logs if necessary: CMakeError.log has a sequence of errors including undefined reference to pthread_create (doesn't appear to be related?)
CMakeOutput.log is much bigger, but the last build command reported is

Run Build Command:/usr/bin/make "cmTryCompileExec3258574548/fast"
/usr/bin/make -f CMakeFiles/cmTryCompileExec3258574548.dir/build.make CMakeFiles
/cmTryCompileExec3258574548.dir/build

...and it appears to have been successful.

Can anyone see what my problem is? Or where I should look? OSRM is new to me (hence my attempt to install & build, so that I can experiment with it), so I am not sure how OSMPBF fits in, or where to look. Also this Linux box has not been used for C++ dev work before, so it is possible a common tool might be missing? (I usually use Linux for Python dev). However the first step from Building on Ubuntu should have installed all the extras (cmake, boost, etc).


texts: Finding overlapping areas in same polygon layer using QGIS?

texts:
I'm working with QGIS to find overlapping area's in a polygon layer.

Some more info :
The Layers contains multiple circles ( buffers ) and i need to find where there is an overlapping area of 3 or more circles. And based on how many circles overlap they have to get a different color / style.

Does anybody know a fast workflow for this?


texts: Speeding up extraction time of proportion of land cover types in buffer from raster using R?

texts:
I would like to extract land cover data in a buffer of 10 km around objects of class SpatialLines and calculate proportion of each land cover type. So, I used the function extract (package raster). I also used the function crop to crop my raster. In my raster, I have 10 land cover types. Here is my code:

buf_line <- gBuffer(seg_line, width=10000) # seg_line = Lines objects
ha <-extract(x=data_raster,y=buf_line)
# The proportion of each land cover type must be in columns (one column = one land cover type)
ha_1 <-length(ha[[1]][ha[[1]]==1])/length(ha[[1]])
ha_2 <-length(ha[[1]][ha[[1]]==2])/length(ha[[1]])
ha_3 <-length(ha[[1]][ha[[1]]==3])/length(ha[[1]])
ha_4 <-length(ha[[1]][ha[[1]]==4])/length(ha[[1]])
ha_5 <-length(ha[[1]][ha[[1]]==5])/length(ha[[1]])
ha_6 <-length(ha[[1]][ha[[1]]==6])/length(ha[[1]])
ha_7 <-length(ha[[1]][ha[[1]]==7])/length(ha[[1]])
ha_8 <-length(ha[[1]][ha[[1]]==8])/length(ha[[1]])
ha_9 <-length(ha[[1]][ha[[1]]==9])/length(ha[[1]])
ha_10 <-length(ha[[1]][ha[[1]]==10])/length(ha[[1]])

I used this code in a loop lapply to extract land cover data in a buffer of 10 km around 30 000 lines. The problem is that the extraction around 30 000 lines is very slow.

How can I speed up the processing time ?


texts: How to apply a model for each value in the attribute table?

texts:
I have big data on production values of farmers (shapefile). Each farmer has a number of parcels in the attribute table. I try to apply a model for each farmer in the attribute table. The model simply following the steps below:
1- Feature to point (input: parcels)
2- Minimum bounding geometry (convex hull; input: output of feature to point' output)
4- Feature to point (input: minimum bounding geometry output)
5- Point distance (this step is to calculate the distance between the center of minimum bounding geometry and the center of parcels)
Note: The model is also attached for a farmer but we want to apply the model for each farmer. You can see the structure of the attribute table with the farmers code (tootja) as well


texts:
We're having an issue trying to display labels UNDER a heatmap-styled layer, and can't figure out why they are always defaulting to the topmost position on the map.

We're using the 'heat map' style instead of using density maps as they are easier and basically the exact same visual representation of a density map.

We're also using the Google-styled SpatiaLite OSM data for our basemaps:

Is there a setting that controls this that we can't find?


texts: Oracle Spatial: Issue with sdo_nn_distance for canada postal codes

texts:
Spatial data are available in the table st_country_postal_code with index. Inorder to get the distance between postal codes, we used the below query

SELECT sdo_nn_distance (1) distance, z.*
FROM st_country_postal_code z
WHERE z.postal_code = UPPER ('V3G1M2')
   AND z.country_id = '124'
   AND 'TRUE' =
         sdo_nn (z.spatial_data,
                 (SELECT spatial_data
                    FROM st_country_postal_code
                   WHERE country_id = '124' AND postal_code = UPPER ('V2R0N3')),
                 'sdo_num_res=999 unit=' || 'KM',
                 1
                );

For some canadian zipcodes, for the above query I am getting no rows returned like V2R0N3, V3G1M2 , J8H2A7. Is there a constraint for the sdo_nn function to not let calculate the distance?


texts:
I want to make a presentation of a museum using QGIS. I have already designed the room layout in AutoCAD and inserted it in QGIS.

Now I want to insert the images in certain boxes (museum exhibits) so when I click on them it will pop up the image.

Any ideas?


texts: Guys and gals could you please provide me with a sample code as to how to rename a field name in a Feature Layer or Table.  I've been trying to use the FieldInfo object to do so over the last couple of days but it's not happening for me.


texts:
I am trying to "physically" join non-contiguous polygons of the same layer to form one (when the polygons are spaced less than 5m), without distorting the original shape.

This is a building layer, the final goal is to simplify their representation by connecting buildings less than 5m from each other.

I know that ArcGIS offers a tool that meets my needs: "Aggregate polygons" or something, I tried several solutions like "snap geometries to layer" but It doesn't work.

Any ideas?


texts: Planet_osm_point versus Planet_osm_nodes for querying places

texts:
Hello I am using a PostGis Database with OSM Data. After importing the data with osm2pgsql in slim mode I found a table named Planet_osm_nodes that has places also. It seems really compressed and has all the metadata in the tags field.
Now I don't really know which of these two tables I should use regarding completeness of the data and performance.
If anyone of you has any thoughts and experiences on this topic it would be great!


texts: extract layers by name from DXF file to multiple destination files using ogr2ogr command line

texts:
Is this possible to extract layers of a DXF file into multiple files using ogr2ogr command line ?

Let us say my DXF file consists of 2 layers : ROAD and HOME.

I want my output to be some kml files.

So I run the following command :

ogr2ogr -f "KML" plan.kml plan.dxf -s_srs EPSG:27562 -t_srs EPSG:4326

Which works very well.

QUESTION

But how can I specify I want multiples files in output ? Or should I programmatically code this using a loop and extract each layer individually ?


texts: Use the Python console to write layer name in attribute table

texts:
The task I would like to solve is to write an individual shapefile's name into a new column of its attribute table for all shapefiles in a given directory.

My current - disfunctional - code looks like this:

import glob, os, processing
input_path = "/home/me/Folder/"
for layer in glob.glob(input_path + "*.shp"):
   processing.runalg("qgis:fieldcalculator", layer, 'columnname', 2, 25.0, 3.0, True, '@layer_name', os.path.basename(layer) + *_new.shp")

This results, so far, in an indentation error pointing to the last element of the algorithm's elements.

Using Joseph's code provided below (copied from a text editor into the Python console), the algorithm runs over one shapefile, writes the shapefile's name including the .shp-ending into a column, then ends abruptly with the following error:

Traceback (most recent call last):
  File "<input>", line 3, in <module>
  File "/usr/share/qgis/python/plugins/processing/tools/general.py", line 75, in runalg
    alg = Processing.runAlgorithm(algOrName, None, *args, **kwargs)
  File "/usr/share/qgis/python/plugins/processing/core/Processing.py", line 257, in runAlgorithm
    if not param.setValue(args[i]):
  File "/usr/share/qgis/python/plugins/processing/core/parameters.py", line 997, in setValue
    self.value = unicode(obj)
UnicodeDecodeError: 'ascii' codec can't decode byte 0xc3 in position 88: ordinal not in range(128)

The names to be written contain German Umlaute, could this cause any problems of this kind?

After removing the Umlaute, I ran this code:

import glob, os, processing
input_path = "/Some/Folder/"
for layer in glob.glob(input_path + "*.shp"):
    output = "/Some/Folder/Subfolder/" + os.path.splitext(os.path.basename(layer))[0] + "_withName.shp"
    processing.runalg("qgis:fieldcalculator", layer, 'Name', 2, 25.0, 0, 1,  '@layer_name', output)

The result is almost as it should be except for the layers' names written into the attribute table still having the .shp-ending. What's wrong with the field calculator command here?

Following Joseph's advice, the code that eventually did the trick was:

import glob, os, processing
input_path = "/Some/Folder/"
for layer in glob.glob(input_path + "*.shp"):
    output = "/Some/Folder/Subfolder/" + os.path.splitext(os.path.basename(layer))[0] + "_withName.shp"
    processing.runalg("qgis:fieldcalculator", layer, 'Name', 2, 25.0, 0, 1,  """ replace(@layer_name,'.shp','') """, output)

As far as I can tell, the replace-command works like this: replace(@layer_name, 'bit to be replace', 'bit to replace it with') - but what exactly are the triple quotes for, Joseph?


texts: QGIS 2.14 c++ standalone application crashes using QgsMapToolIdentify

texts:
lResultList.append(mpIdentifyTool->identify(lCoord.x(),lCoord.y(),QgsMapToolIdentify::TopDownStopAtFirst, QgsMapToolIdentify::VectorLayer));

leads to crash of app. It happens when layers provided by "memory" provider exist on canvas. With layers provided with "ogr" code works normally. Is there way to fix this? Or is there another way to select feature by clicking on it with mouse?


texts: Appended attribute returning "Attribute index 9 out of bounds [0;9]"

texts:
I am trying to write a PyQGIS script that copies information from one feature into an output feature and adds a field.

This is how I'm initializing my fields for my output layer

inlayer = processing.getObject(input_layer)
provider = inlayer.dataProvider()
fields = provider.fields()
fields.append( QgsField("NewField", QVariant.Double) )
writer = QgsVectorFileWriter(output_layer, 'UTF-8', fields, QGis.WKBPoint, inlayer.crs())

I'm then looping through the features and adding output features to my writer

outFeature = QgsFeature(fields)
outFeature.setGeometry(setting the geometry... nothing to see here)
outFeature.setAttributes(inFeature.attributes()) #copying attributes from the input layer
outFeature.setAttribute("NewValue", newValue) #this line is what is causing the error
writer.addFeature(outFeature)

When I run the script I am getting an steady stream of error messages in the log messages panel under the General tab

  Attribute index 9 out of bounds [0;9]

If I check the attributes table for the output layer that I just created, I have all the fields that I'm expecting. I have all the attributes from the input layer, and I have a column for the NewValue, but I do not have any data populated in the NewValue column. The NewValue column is the 10th column, so I'm confused why my script is failing to write to that column saying index 9 is out of bounds because I would think that the NewValue column would be index 9. What am I missing or doing wrong?


texts: Inverting donut polygon in ArcGIS for Desktop?

texts:
I extracted the study area polygon using raster to polygon tool in ArcGIS from a binary raster mask. However, I can only extract a polygon with fills outside of my study area. I tried using editor tool to draw along the needed area, but fail to plot precisely as extracted.

Is there any way to invert the polygon values so that it could be used as a mask for the central area?


texts:
I have recently installed QGIS 1.8

when I make a map with the adequate layer I switch to the map composition menu and try to insert a map in the white page.

I draw the rectangle but nothing appears although I tried the 3 different option

I also ticked on the show on map for all the layers (show in overview)

can you tell me if this is an installation problem (I have windows XP professional and downloaded the french version) or a trivial problem I am not being able to resolve?


texts: Using unioned feature class to quickly get attributes from multiple feature classes at once?

texts:
I've been asked to create a feature class by unioning 18 polygon feature classes together, and retaining only certain attributes. The goal is to have a single feature class that can be queried once to get multiple attributes (e.g. zip code, technician territory, municipality, county), instead of having multiple feature classes that each need to be queried.

The user plans to feed a lat/long into some .NET application and get information about that location -- he can get the ~20 attributes he needs querying individual features, but would like a "one stop" solution. Apparently the data is taking up to ten seconds to pull from a number of map services, which is unacceptably long to the group who needs the application, and sending only one request would significantly improve performance.

The technical process of creating this feature class isn't really a problem. However, I am concerned about it from a "best practice" perspective -- it's a bulky bunch of polygons, and we are duplicating data that must be constantly re-created to capture any changes in the source data. Again, that's not a technical issue, I can run this as a scheduled task... but the potential time lag between updates makes me slightly uncomfortable about the accuracy of this new layer's results.

So my concern is really whether this is a good idea, and whether this is a typical solution for GIS web applications that query multiple layers simultaneously.

Are there some alternative ideas that I can propose that would quickly provide the necessary data, while letting us maintain our ZIP code, county, territory, etc. polygons as unique feature classes that aren't duplicated somewhere else?

I am thinking of something that's sort of like a "view" in a non-spatial database, holding only the particular attributes he's interested in but including a spatial context since it's location-based -- but I don't have enough database background to know whether that's possible.

The data is in an ArcSDE database, and he's querying Esri map services through .NET (probably ArcObjects).


texts: Error using gdal_calc.py with gdal 2.2.1 on linux (undefined symbol: GDALCreateDriver)

texts:
My gdal libraries work properly with any other function used. However, using gdal_calc.py I run into troubles. I use gdal from the shell (linux) with this command:

gdal_calc.py -A /somepath/somefile.tif -B /somepath/somefile2.tif --outfile=/somepath/somefile3.tif --calc="A*B" --type UInt32 --overwrite

The error I get is:

Traceback (most recent call last):
  File "/usr/bin/gdal_calc.py", line 53, in <module>
  from osgeo import gdalnumeric
File "/usr/lib/python2.7/dist-packages/osgeo/gdalnumeric.py", line 1, in <module>
  from osgeo.gdal_array import *
  File "/usr/lib/python2.7/dist-packages/osgeo/gdal_array.py", line 28, in <module>
  _gdal_array = swig_import_helper()
File "/usr/lib/python2.7/dist-packages/osgeo/gdal_array.py", line 24, in swig_import_helper
_mod = imp.load_module('_gdal_array', fp, pathname, description)
ImportError: /usr/lib/python2.7/dist-packages/osgeo/_gdal_array.x86_64-linux-gnu.so: undefined symbol: GDALCreateDriver


texts: "Locate features along routes" tool and line features - how to use XY tolerance?

texts:
I am trying to use the "locate features along routes" tool to map guardrails to a route, but I'm seeing some very strange behavior.  As a test case, I'm selecting a single guardrail (pictured) that is about 8 meters from the route.

If I set the XY tolerance to 5m, it works perfectly and gives me a single record in the resulting event table with a from and to measure.  But if I set the XY tolerance to 10m, the resulting event table is empty.  The docs say this is a cluster tolerance when line features are used, so why would 5m work but 10m not?  I've tried several values:

2m = doesn't work

3m = works

4m = works

5m = works

6m = works

7m = doesn't work

10m = doesn't work

What is causing this?  I'm trying to process thousands of records so playing with each individual guardrail's XY tolerance is not going to work.


texts:
I am opening a QGIS project using the following snippet:

new_project = False
if not os.path.exists(self.qgis_project_path):
    prj_template_path = os.path.join(self.qgis_template_path, 'project.qgs')
    shutil.copy(prj_template_path, self.qgis_project_path)
    new_project = True    #load default layers because the user is opening a particular project for the first time

QgsProject.instance().read(QFileInfo(self.qgis_project_path))
if new_project:
    for layername in self._settings['vector_layers']:
        uri = QgsDataSourceURI()
        uri.setConnection(self._settings['db_host'], self._settings['db_port'], self._settings['db_dbase'], self._settings['db_user'], self._settings['db_pass'])
        uri.setDataSource('p_{}'.format(self._project[2].replace('-', '_')), layername, 'geom', '')

        lyr = self._iface.addVectorLayer(uri.uri(), layername, 'postgres')
        if not lyr.isValid():
            QgsMessageLog.logMessage(u'Failed to open the layer: {}'.format(layername), 'Skytek')
        styleload = lyr.loadNamedStyle(os.path.join(self.style_path, '{}.qml'.format(layername)))
        if not styleload[1]:
            QgsMessageLog.logMessage(u'Error loading style: {}'.format(layername), 'Skytek')

    for raster in self._settings['raster_layers']:
        pass

    self._iface.mapCanvas().setDestinationCrs(self._local) #qgis3: QgsProject.instance().setCrs(self._local)
    self._iface.mapCanvas().setMapUnits(QGis.Meters) #qgis3: QgsProject.instance().setDistanceUnits(QgsUnitTypes.Meters)
    QgsProject.instance().setDirty(True)
    QgsProject.instance().write(QFileInfo(self.qgis_project_path))

I don't understand why every SECOND time I call this, the loaded project shows no layers in the legend. The project file is copied to the user folder only the first time, also the layers are created only the first time. The only line that is called subsequently is:

QgsProject.instance().read(QFileInfo(self.qgis_project_path))

And like I said, layers are added to the legend only every second time the project is read. I noticed when I reopen an existing project file in QGIS with "Project->Open project" it works as expected.


texts:
I have a problem with this code :

for layer in layers :
   #surcharge avec le filtre permanent
   if QgsExpressionContextUtils.layerScope(layer).variable('layer_filter')!=None:
     sql=u"(" + sql + u") AND " + QgsExpressionContextUtils.layerScope(layer).variable('layer_filter')

When a variable is found it is added to the sql variable. It is however kept for every layer even if there is no filter variable for them.

What can be done so the SQL filter corresponds to a single layer?


texts: How to create one KML per polygon (from one Shapefile)?

texts:
I have a shapefile containing hundreds of polygons. I need to create a kml of each of the polygon. Is there any way that I may not have to deal each polygon individually and my polygon file may give me individual kml of each polygon.

I'm using ArcGIS 9.3 or QGIS.


texts: For loop for NDVI of Landsat 8 for specific path/rows with least cloud cover. TypeError: Cannot read property 'type' of null

texts:
I am trying to filter through Landsat 8 SR images to gab seven path/rows of the least cloudy scenes from the growing season in 2016 then mask out the clouds, and calculate the NDVI. Then export each image to my google drive. However, I am receiving the error

  TypeError: Cannot read property 'type' of null.

Can anyone see my mistake?

[![enter image description here][1]][1]


texts: How to conduct an NDVI and regression with ArcMap?

texts: I would ideally like to execute an NDVI on a sample of NAIP images and then use the NDVI scores per district or whatever my aggregate zone is, probably census tract.  My goal is to identify census tract parcels which are good for development based on, my NDVI score, and other socio-economic demographics such as age, population, income, businesses, etc.  I know NDVI can be calcuated in ArcMap, how would I translate the values that I receive per pixel to a vector data set?


texts:
My task is to model the surface of a lake bed, calculate volume of it and visualize it as 3D a surface.

I have acquired depth measurements in a small lake using echosounder. The precision should be pretty high (<5 cm for depths and 3-5 cm for RTK GPS position). Points are placed regularly with 2.5 m intervals between them.
I also have digitized shoreline from high resolution orthomosaic.

Which would be the best interpolation algorithm for such tasks? (QGIS, GRASS, SAGA, maybe some other free software does such things?). The biggest challenge has been to find an algorithm, in which I can define the shoreline as 0 value which has to be taken in account.


texts: What are the elements that a good map should have?

texts:
I have drawn a map with some coordinate(longitude and latitude). What are good points that a map should have eg labels, legend, etc.
 Edit:
The map if for the purpose of locating different objects and showing their interconnections through landmark ,color etc. The map should provide clue about itself.


texts: Delete a line between two polygons so I have one polygon

texts: I have two polygons that are coincident.  I want to merge them in to one polygon and delete the line that is between the two polygons; the merge didn't do it; it kept the line in between them.  How do I get rid of this stupid line in ArcGIS 10.


texts: Writing code for monthly NDVI medians in Google Earth Engine?

texts:
My goal is to get monthly NDVI medians in Google Earth Engine. This helpful post got me started: reduce image collection to get annual monthly sum precipitation

However, I am using Landsat (not MODIS). This script does run but the map layer has "no bands to visualize".

How can I visualize the bands? (Maybe the issue is that this dataset only has 1 band?).
How can I add all 12 map layers, one per month? (perhaps with stacking?)

var ndvi = ee.ImageCollection('LANDSAT/LT5_L1T_32DAY_NDVI');

var months = ee.List.sequence(1, 12);

// Group by month, and then reduce within groups by mean();
// the result is an ImageCollection with one image for each
// month.
var byMonth = ee.ImageCollection.fromImages(
  months.map(function (m) {
    return ndvi.filter(ee.Filter.calendarRange(m, m, 'month'))
                .select().median()
                .set('month', m);
}));
print(byMonth);

Map.addLayer(ee.Image(byMonth.first()));


texts: L.geoJson() - Error: Bounds are not valid

texts:
I am iterating through an array of zipcodes and passing each to Mapbox Geocoding API. With every iteration, on success, I pass to L.geojson().  After the last iteration I call map.fitBounds(jsonLayer).

Then this error:

  ...Uncaught Error: Bounds are not valid

Here is my full code:

var zipList = ["90210", "03103", "11206", "79109"];

function findMe(zipList) {
    var listLength = zipList.length;
    for (var i = 0; i < listLength; i++) {
        var myRequest = "https://api.mapbox.com/geocoding/v5/mapbox.places/" + zipList[i] + ".json?access_token=" + L.mapbox.accessToken + "&limit=1";
        $.ajax({
            method: 'GET',
            url: myRequest,
            success: (function (data) {
                jsonLayer = L.geoJson(data, {
                    onEachFeature: onEachFeature
                });
            })
        })
    }
    map.fitBounds(jsonLayer);
};

Does this have something to do with Ajax and asynchronous?  What can I do to avoid this error and fit map to comprehensive marker set?

--

Edit: I have also tried pushing lat and long coordinates to an array before calling map.fitBounds(markerArray); Same error.

function onEachFeature(feature, layer) {
    markerArray.push(L.marker([feature.geometry.coordinates[1], feature.geometry.coordinates[0]]));
}


texts: Prevent panning and zooming beyond map edges in OpenLayers 4.6.4

texts:
I have a web page with map that is using OpenLayers 4.6.4 and I want to prevent the user from panning or zooming beyond the edge of the map, like this:

There are plenty of answers for this for OL3, but that is very different than OL4. Items likes "restrictedExtent" that existed in OL3 are no longer there in OL4.

The only customizations I have added to the map are the fullscreen and scale controls and an onclick handler as I am still new to OL & GIS (though not to web development).

OL is getting the data from a GeoServer instance via WMS if that makes any difference.


texts: OGR2OGR -overwrite mssql table is appending instead

texts:
I am running

Command '['C:\\OSGeo4W\\bin\\ogr2ogr.exe', '-progress', '-skipfailures', '-overwrite', '-lco', 'UPLOAD_GEOM_FORMAT=wkt',
 '-t_srs', 'EPSG:28354', '-f', 'MSSQLSpatial', 'MSSQL:Server=xxx;Database=xxxGIS_Test;Trusted_Connect
ion=True', 'D:\\GIS_SourceData\\xxx\\HY_WATERBODY_STATE.tab', '-nln', 'DSE_HY_WATERBODY_STATE']'

But this is appending records to the table in MSSQL instead of overwriting the whole table?

Is there some switch to just tell it to update records/insert new ones or what should I do to get -overwrite to not append?

The _backup view is the original that was just backed up before running the command.

UPDATED EXAMPLE

WITH DEBUG

Is there a way in OGR to delete the table or delete the records in the table before overwrite?


texts: Finding area from overlay feature class?

texts: I want to extract the counties (green area) overlaying with another feature class (pink area) using ArcGIS 10.3.1 (Advanced). I am using spatial join and find the Join_Count > 0 to select out the involving counties. Is there a better way to do this process or what are some alternatives?


texts:
I've been provided with a dataset containing some coordinates, but I don't recognise the coordinate system used. They look like standard UTM coordinates but the values don't match up with what I expect, which are locations in Scotland, UK.

For example, a random record gives GIS_X as 349275.0 and GIS_Y as 636198.0. Using these coordinates along with zone 30, plugging these values into http://www.rcn.montana.edu/Resources/Converter.aspx results in WGS coordinates not far from the equator, which is evidently wrong.

The datasource is checked daily for various errors so I'm fairly certain the data provided to me is correct; any ideas for how these values might make sense, given the expected result?


texts: I am comparing ASTER and MODIS LST images. For this I converted the images from HDF format to GeoTIFF format using ERDAS 2013 of my study area.  MODIS image is reprojected to the same projection of aster image and the shape file of the study area also. But the two satellite images does not overlap with the shape file. How can I solve this problem?


texts: Can MOOCs Count towards the GISP course points?

texts: I am looking to obtain my GISP, and I am wondering if Massive Open Online Course (MOOCs)  can count towards the GISP course points if the MOOC is related to GIS?


texts: Can't preview a shapefile in GeoServer 2.2.4

texts: I have problem with previewing layers in GeoServer. I installed Geoserver 2.2.4 (following instructions in manual) and everything was going well. I managed to login to it and create my own new workspace. I've created a new store and published that layer but when I try to preview that layer with openLayers all I see is empty map with no features on it. I've also tried with example on http://docs.geoserver.org/2.2.4/user/gettingstarted/shapefile-quickstart/index.html  to see if the problem is with my shapefile but the same thing happens. I really don't know why because all I did was the same as it is written in manual. I hope somebody can help me with that. Thanks


texts: How to not show labels that are under polygons of another layer in QGIS 3.6.1

texts:
I am editing the polygon layer by adding new polygons and moving them to new locations. When I place a polygon in a new location I would like any labels from the two point shapefiles that happen to be under the polygon to disappear, while the labels of the points not covered by any polygons remain. This is purely for visual reasons as I am trying to cover all the points with polygons and eliminating the labels beneath existing polygons would make it more obvious which points still need to be covered. Is this possible in QGIS 3.6.1?

The closest thing I have been able to find on SE or anywhere else is the 'discourage labels from covering feature' option (Layer Properties -> Labels -> Rendering -> Obstacles) but even when I moved the slider all the way to the high weight side it didn't even move the labels, let alone make them disappear.


texts: How can I use GDAL installed with QGIS from command line?

texts:
I have QGIS 3.10 and I can run GDAL using OSGeo Shell. I want to be able to use GDAL commands from command line. For now, when I run gdalinfo --version I get message 'gdalinfo' is not recognized as an internal or external command, operable program or batch file.
How can I fix this?


texts: Geostatistical Modelling / What is the right model for my data?

texts:
I try to find a model for my data. I explored OLS regression models but I just get an R² of around 0.18 and I think I should not go with that. I also noticed that my data could be better modeled with a surface than with a line, if that makes any sense (I am a geostatistics beginner!).
I have point data of events, which represent a GPS location of a boat during dolphin watching activity. The point data was aggregated into fishnet bins of 500x500m. So Join_Cou_1=3 means that I have 3 point events in a 500x500m square.
Taking the distance to the closest port as an explanatory variable, I notice that I have a relation that could be described better with a bell shaped surface then a bell shaped curve. Even though I have a high concentration of Events at a distance range from 8-10 km, I also have areas there with a low event concentration. At very low and very high distance I have ONLY low event concentration. If I would include No-Events (I computed random points and selected those which dont share a 500x500m square with an event) it would be a continous line below in the Y-axis (I have No-Events at every distance but I ONLY have No-Events at distances below 2km at distancens higher than 16km). (Picture 1 shows the scatter Plot of the Event bins and distance)
With my second explanatory variabel, bathymetry (gridcode) it seems to be similar. Rather than drawing a diagonal line through the middle of the data I would like to draw a diagonal line in the topleft and say that everything below the line is probable and everything above not. (Picture 2 shows the scatter Plot of the Event bins and the bathymetry=gridcode)
Does that have to do with an accumulated probability distribution?
Can I model the probability with Krigin first?
Does it make any sense at all?
Any hints?


texts: Can you identify attributes of a shapefile on the gps?

texts: My company uses Trimble GeoExplorer 6000 GPS units and pathfinder office software. Is there a way once a shapefile is created to bring up the attribute table on the gps unit? We have been loading previously recorded data to the gps as a background file. When I return to a site this allows me to see where the data was collected but I can't see any information to identify the attributes. Any help would be greatly appreciated!


texts: Point markers with icon/picture symbols showing white area around them on ArcMap?

texts:
I'm trying to use some icons (pictures) but i haven't been successful yet.

This is where i got the icons: http://mapicons.nicolasmollet.com/category/markers/offices/

I want to build a 1 meter by 2 meter tourism map of a specific area of Lisbon using these icons to ilustrate where are the restaurants/bars/banks/supermarkets/etc...

What happens when i add an icon in the arcmap is this:icones2.jpg

They have this white area around them but i would really like to take it away. I tried to get some help from some photoshop geeks but whatever the solution, when i upload them to arcmap they just don't change and still have the white edge.


texts: Python error when running SpatiaLite with OSGB(27700)

texts:
I have a project open in QGIS and am trying to import an Ordnance Survey Address Point layer into SpatiaLite. I have recently installed the Qspatialite plugin.  I open SpatiaLite in the DB Manager then use the Import QGIS layer to choose the Address Point layer.  I then get the following Python error code data and have no idea of what is happening.

An error has occured while executing Python code:

Traceback (most recent call last):
  File "/Users/Nlaw/.qgis//python/plugins/QspatiaLite/importQgis/importQgis.py", line 64, in run
    if not self.parent.db.uploadQgisVectorLayer(layer, tableName=self.tablenames[i], srid=self.srid, selected=self.selected):
  File "/Users/Nlaw/.qgis//python/plugins/QspatiaLite/Classes/Database.py", line 249, in uploadQgisVectorLayer
    header,data=self.executeQuery("""INSERT INTO "%s" VALUES (%s,%s)"""%(tableName,','.join([unicode(value).encode('utf-8') for value in values_auto]),','.join('?'*len(values_perso))),tuple([unicode(value) for value in values_perso]),commit=True)
  File "/Users/Nlaw/.qgis//python/plugins/QspatiaLite/Classes/Database.py", line 45, in executeQuery
    cursor.execute(query,params)
IntegrityError: AddressPoint.Geometry violates Geometry constraint [geom-type or SRID not allowed]

Python version:
2.7.2 (default, Jun 20 2012, 16:23:33)
[GCC 4.2.1 Compatible Apple Clang 4.0 (tags/Apple/clang-418.0.60)]

QGIS version:
1.8.0-Lisboa Lisboa,

Python path: ['/Applications/QGIS.app/Contents/MacOS/../Resources/python', '/Users/Nlaw/.qgis//python', '/Users/Nlaw/.qgis//python/plugins', '/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins', '/Library/Frameworks/SQLite3.framework/Versions/3/Python/2.7', '/Library/Frameworks/GDAL.framework/Versions/1.9/Python/2.7/site-packages', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python27.zip', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-darwin', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-mac', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-mac/lib-scriptpackages', '/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-tk', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-old', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-dynload', '/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/PyObjC', '/Library/Python/2.7/site-packages', '/Applications/QGIS.app/Contents/Resources/python/plugins/fTools/tools', '/Users/Nlaw']


texts: Clipping raster mosaic changes range of cell values to -3.4E+38 --- +3.4E+38

texts:
I'm trying to clip a DSM mosaic I've made with a shapefile on ArcMap 10.7 but when I use extract by mask tool, the output is just empty and when I change the no data colour to red, all I get is a red square. The cell value ranges are -3.4E+38 --- +3.4E+38. The strange thing is that the corresponding DTM clips absolutely fine. Both the coordinate systems are the same (British National Grid).
I'm pretty new to GIS, as you may be able to tell, does anyone have any idea of what I am doing wrong?


texts: Export several configured atlases with a Python script

texts:
I'm using @Alexandre Neto's Python code (from July 8 2019) to automatically export predefined atlas layouts, see https://gis.stackexchange.com/questions/272839/export-a-configured-atlas-with-a-python-script-command-line.

I'm using QGIS 3.10.2 on a Windows 10 PC.

export_atlas.py is called by export_atlas.bat.

export_atlas.py (Edit: I have tried both / and \ in project_path and output_folder, both giving the same errors and wrong destination folder)

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# linja over bruker vi for å kunne bruke æøå uten advarsler.

import os
from qgis.core import  QgsApplication, QgsProject, QgsLayoutExporter

def export_atlas(qgs_project_path, layout_name, outputs_folder):

    # Open existing project
    project = QgsProject.instance()
    project.read(qgs_project_path)

    print('Project in ' + project.fileName() + ' loaded successfully')

    # Open prepared layout that as atlas enabled and set
    layout = project.layoutManager().layoutByName(layout_name)

    # Export atlas
    exporter = QgsLayoutExporter(layout)
    settings = QgsLayoutExporter.ImageExportSettings()
    exporter.exportToImage(layout.atlas(),outputs_folder, 'jpg', settings)

def main():
    # Start a QGIS application without GUI
    qgs = QgsApplication([], False)
    qgs.initQgis()

    project_path = 'C:\\adhoc\\lst/lstp_shp.qgz'
    output_folder = 'C:\\adhoc\\lst\\n5g'
    layout_name = 'n5 NGO48 sone G'
    print('Starter atlas-eksport')
    export_atlas(project_path, layout_name, output_folder)

    # Close the QGIS application
    qgs.exitQgis()

if __name__ == "__main__":
    main()

export_atlas.bat

@echo off
chcp 65001

REM Change OSGEO4W_ROOT to point to the base install folder
SET OSGEO4W_ROOT=C:\OSGeo4W64
SET QGISNAME=qgis
SET QGIS=%OSGEO4W_ROOT%\apps\%QGISNAME%
set QGIS_PREFIX_PATH=%QGIS%
REM Gdal Setup
set GDAL_DATA=%OSGEO4W_ROOT%\share\gdal\
REM To find proj.db
set PROJ_LIB=%OSGEO4W_ROOT%\share\proj\
REM Python Setup
set PATH=%OSGEO4W_ROOT%\bin;%QGIS%\bin;%PATH%
SET PYTHONHOME=%OSGEO4W_ROOT%\apps\Python37
set PYTHONPATH=%QGIS%\python;%PYTHONPATH%
REM For riktig PyQt med dll-er
call "C:\OSGeo4W64\bin\qt5_env.bat"

REM Launch python job
C:\OSGeo4W64\apps\Python37\python export_atlas.py

This works.

However: I get error-messages and the jpg-files are stored in C:\adhoc\lst and not in C:\adhoc\lst\n5g which is set as output_folder in the export_atlas.py (see above).

Active code page: 65001
Starter atlas eksport
libpng warning: iCCP: known incorrect sRGB profile
libpng warning: iCCP: known incorrect sRGB profile
Project in C:/adhoc/lst/lstp_shp.qgz loaded successfully
ERROR 6: The JPEG driver does not support update access to existing datasets.
ERROR 6: The JPEG driver does not support update access to existing datasets.
ERROR 6: The JPEG driver does not support update access to existing datasets.
...

Is it possible to get rid of the error and/or to save the jpg-files in the n5g-folder?


texts:
In QGIS 2.18, I had a python script that added annotations. I can't seem to get the annotations to add in 3.0:

for feature in iter:
            # retrieve every feature with its geometry and attributes
            # fetch geometry
            geom = feature.geometry()
            # show some information about the feature
            if geom.type() ==  0: #QgsWkbTypes.Point:
                print("In add annotation...")
                point = geom.asPoint()
                utmPoint = utm.from_latlon(point.y(), point.x())
                # fetch attributes
                attrs = feature.attributes()
                # attrs is a list. It contains all the attribute values of this feature
                waypointText = attrs[0]
                textItem = QgsTextAnnotation (iface.mapCanvas())
                crs = waypoint_layer.crs();
                document = QtGui.QTextDocument()
                document.setHtml("<strong>" + str(waypointText) + "</strong>")
                symbol = QgsMarkerSymbol()
                symbol.setSize(0);
                textItem.setDocument(document)
                textItem.setMapPosition(QgsPointXY(float(utmPoint[0]), float(utmPoint[1])))
                textItem.setMapLayer(waypoint_layer);
           #     textItem.setFrameSize(QSizeF(document.size().width() + 5,document.size().height() + 2))
           #     textItem.setFrameBorderWidth(2)
           #     textItem.setFrameColor(QColor("#e31a1c"))
           #     textItem.setFrameBackgroundColor(QColor("#d7d7d7"))
                 textItem.setMarkerSymbol(symbol)
           #     textItem.update()

Commented out textItem lines no longer work in 3.0, but even without setting frame specific features, the annotation doesn't show up.

Any examples of programatically adding an annotation to a map with QGIS 3.0?


texts: Display label attributes created in the GUI in QGIS Python Console for later use within PyQGIS

texts:
I want to define a label for a layer created dynamically in a Plugin with PyQGIS.
I managed this, as there is a big amount of good examples you can find in the internet.
When creating the symbology for this dynamic layer, i used this way:
1.) Create Symbology in QGIS for an existing layer
2.) Use the following command in the Python Console:
print(iface.activeLayer().renderer().symbol().symbolLayers()[0].properties())

3.) Used some parts of the output as a template, so that i know which paramters have to be changed in python to receive the same result from Python as when created in QGIS.
Does there exist a similar command for showing all the attributes for the label created before in QGIS?
I think this has someting to to with QgsPalLayerSettings in combination with the format function, but i have no idea of how to get the information.


texts: What is the sequence of select query execution in Mapinfo?

texts:
I have confusion with the following query:

Query 1:

 Select * from TABLE where vendor = "vendor1" OR
                                    "vendor2" OR
                                    "vendor3" AND site_status = "staus1"
 into sel

Query 2:

Select * from TABLE where ( vendor = "vendor1" OR "vendor2"
                            OR vendor ="vendor3")
                            AND  site_status = "staus1"
into sel

Query 3:

 Select * from TABLE where vendor = "vendor1" OR
                           vendor = "vendor2" OR
                           vendor ="vendor3" AND site_status = "staus1"
 into sel

Query 4:

Select * from TABLE where (vendor = "vendor1" OR vendor= "vendor2" ) AND
                          (site_status = "status1" OR site_status = "status2") into sel

How to find sequence of execution of query?


texts: Merging attribute tables from feature classes into single table?

texts:
I have a large number (approx. 1,100) polygon feature classes.

What I would like to do is to merge the attribute tables from each of these feature classes into a single table, in the simplest and fastest way possible.

This will be the first step in a script or model that then uses that table with all polygon values from all layers to perform some tabular joins and selections.  Each time I run this script, I will want to re-do the table merge, as the underlying polygon feature classes may change.  If it's possible to write the merged table into memory, for speed, that would be best.

I'm using ArcGIS 10.1.  I think I can do this easily enough in VBA/VB.net, but it might be simplest to set it up as a ModelBuilder model, so that it's easy to share with colleagues.


texts: Calculating Long wave radiation using Google Earth Engine?

texts:
How to perform computations on two image collections( here e0 and Ts). I have to use the formula RL_up = stfn*e0*T^4 for computing the long wave radiation.

//outgoing long wave radiation
var stfn = 5.67/100000000;

var Ts = surfacetemp.filterBounds(ROI).filterDate( '2018-01-01', '2018-12-31').sort('system:time_start', false);

var lai = LAI.filterBounds(ROI).filterDate( '2018-01-01', '2018-12-31').sort('system:time_start', false);

print(Ts);
//Map.addLayer(Ts);
var e0 = function(image){
  var e01 = image.expression(
  "0.95 + 0.01*lai",
  {
    'lai': image.select('LAI'),

  });`
  return(image.addBands(e01.rename("E0")));
};

var dataset = LAI
                  .filterDate('2018-01-01', '2018-12-31')
                  .filterBounds(ROI)
                  .map(e0);

print("dataset",dataset);

var e0 = dataset.select("E0"); // To isolate just the albedo band
print("e0",e0); // Look at it

//Long wave radiation
var RL_up = function(image){
  var RL = image.expression(
  "(5.67/100000000)*(T**4)*e0",
  {
    'T': image.select('LST_Day_1km'),
    'e0':dataset.select("E0")

  });
  return(image.addBands(RL.rename("RL_up")));
};

var dataset2 = surfacetemp
                  .filterDate('2018-01-01', '2018-12-31')
                  .filterBounds(ROI)
                  .map(RL_up);

print("dataset2",dataset2);

var RL_up = dataset2.select("RL_up"); // To isolate just the albedo band
print("RL_up",RL_up); // Look at it


texts:
I'm trying to retrieve data from remote WFS server which requires basic https authentication using GeoTools.
I  defined connection parameters for username and password like

connectionParameters.put("WFSDataStoreFactory:USERNAME", "username" );
connectionParameters.put("WFSDataStoreFactory:PASSWORD", "pass" );

However I get

"Server returned HTTP response code: 401 for URL: http://...".

What is the solution for basic https authentication in GeoTools WFS support.


texts: QGIS 2.18: Generating simple report on spatial objects falling into an area of interest

texts:
I am looking for a way to generate a simple report from spatial object falling into the selected area of interest. Report would need to be in form of a table which includes all values from the attribute tables of given polygons.

Any ideas how to go around it most effectively?

I tried Reporter plugin in Qgis but it doesn't do the work exactly. Are there any plugins or tools which could achieve the result?


texts: ERDAS Imagine signature picker colour changing RGB numbers?

texts:
I've noticed something with the signature picker in ERDAS and I was wondering if it is a bug or a feature.

When I create classes I assign a colour to make them easily readable in the classified image. For example blue for the sea, white for clouds and gray for urban areas etc.

However the RGB numbers created by merging selected signatures, change depending on what colour I pick.

Is this supposed to happen or is it a bug?


texts:
I am writing a search client, to query postgis db.

The search query will be,

SELECT tags->'name' AS name,ST_Transform(way,4326) AS pt_lonlattext FROM planet_osm_point WHERE lower(tags->'name') LIKE '%school%'

Now, I need the results in such a way that,

fully matched results should be given priority than partially matched ones.

Is there a way to re-construct the above query to get the desired results?


texts: Implementing custom IRenderer for ArcGis 10.0/WPF

texts:
EDIT: Updates in the bottom

As a followup to my earlier question, I'm trying to implement custom symbol rendering for my WPF/ArcGIS application.

I'm accessing layers that have been published on the ArcGIS server, with pre-existing symbology already in place. I want to keep this symbol, and add text to it, since FeatureLayers don't automatically display their data labels. I'm running into several problems.

I've defined my renderer thus:

public class EntityLabelRenderer : IRenderer
{
    public Symbol GetSymbol(Graphic graphic)
    {
        var objectId = graphic.Attributes["OBJECTID"]; // always set.
        var graphicFromServer = //Use a QueryTask to fetch the object;

        return graphicFromServer.Symbol; //add text
    }
}

When I load my layer, GetSymbol is called for each object, but the Graphic object I get in the call is very barebones - the Symbol property is null, and the only values I can get are the ID. So I use that ID, as I show above, to fetch the full Graphics object with all its properties using a QueryTask. But even the one I get now has its Symbol set to null. Additionally, I don't know what the default data label for this object is, as set on the server. I have all the fields, but have to decide myself which to display.

So my question is:
1) Can I get, for a Graphic object fetched from an ArcGIS server, the default data label as defined in the MXD?
2) Get I get the default symbol for an object, as defined in the MXD on the server? Can I get that, then create a composite symbol comprised of both that symbol and an additional TextSymbol?

3) Is there another, better way? The only one I could think of was to create a custom GraphicsLayer and add TextSymbols to it manually for each object of mine, but this seems insanely complicated.

----- Update -----

I've gone over this and tried some things, and here's where I am in my attempts to get this working:

If I wait until my FeatureLayer's Initialized event is fired, I see that the layer comes populated with an IRenderer for the server-configured symbols. This is great, since I now have the symbol I want to display, as well as the layer's definition of the default text (to display as a label).

So now I can create my own IRenderer which wraps this server-generated IRenderer (usually a UniqueValueRenderer or UniqueMultipleValueRenderer) and adds the label. I can see two approaches:
1) Hacking into the existing renderer, changing its ControlTemplate to add my label, or
2) Somehow returning a composite Symbol that contains both the existing Symbol and my own TextSymbol.

Option #2 is better, but it seems that unlike the Flex and Java SDKs, the WPF SDK doesn't have a CompositeSymbol object. I'd rather not go for option #1 before I'm sure that a CompositeSymbol can't be implemented somehow. Does anyone know how I can go about this?


texts: Multiple shapefiles to one GeoPandas GeoDataFrame

texts:
I'm trying to vertically concatenate multiple shapefiles that reside in one directory, into one geodataframe. They all have the same table schema and same CRS. I am trying to use the code from this Stack Overflow post, but am getting an error 'WindowsPath' object has no attribute 'split'.

from pathlib import Path
import geopandas as gpd
import pandas as pd

folder = Path("shp")

gdf = pd.concat([
    gpd.read_file(shp)
    for shp in folder.glob("*.shp") ## error message points to this line ###
]).pipe(gpd.GeoDataFrame)

I'm new to programming and have spent more time than I care to admit trying to figure this one out. I'm open to using alternative Python modules as long as it ends with a geodataframe.


texts: Same crs, different spatial reference system units, shape files / QGIS

texts:
I have two files, same crs but different reference system units and different scale. Long story short, the shp files won't overlap :(

I've read some answers to what seems to be the same problem but I just can't fix it.

Layer 1
Layer Spatial Reference System:
WGS 84 / UTM 21S
In layer spatial reference system units:
xMin,yMin 368398.86,6105694.51 : xMax,yMax 852832.93,6265658.88

Layer 2
Layer Spatial Reference System:
WGS 84 / UTM 21S
In layer spatial reference system units:
xMin,yMin -57.9891,-36.9946 : xMax,yMax -50.0697,-34.4732

Can someone help me with this problem?


texts: Use arcpy to check whether feature class geometry has ZM

texts:
I am looking for a programmatic way using arcpy to determine whether or not a feature class (FC) has Z and/or M geometry. It's easy enough to get the geometry type using the code below, but this just tells me the base type (point, line, polygon) and does not have details about ZM. I would rather not have to create a cursor and read the actual feature geometries (it should not be necessary, since the entire FC must have the same geometry).

import arcpy
desc = arcpy.Describe("D:\\MyData\\Scratch.gdb\\polyzm")
geometryType = desc.shapeType
print geometryType

which returns this result:

Polygon


texts:
I'm trying to figure out how to mathematically derive the common points of two intersecting circles on earth's surface given a center Lat/Lon and a radius for each point.

For example, given:

Lat/Lon (37.673442, -90.234036) Radius 107.5 NM
Lat/Lon (36.109997, -90.953669) Radius 145 NM

I should find two intersection points with one of them being (36.948, -088.158).

It would be trivially easy to solve this on a flat plane but I don't have any experience solving equations on an imperfect sphere such as the earth' surface.


texts: How to apply new renderer in ArcObjects?

texts:
I try to create a new renderer for a geo feature layer in ArcGIS 10. (My project is an asp.net based webservice)
There is one solution I got from a tutorial that works, but only for a ISimpleRenderer

    // Create an ArcObjects color object with the color set to blue for use by the layer's renderer
    IRgbColor rgbColor = (IRgbColor)serverContext.CreateObject("esriDisplay.RgbColor");
    rgbColor.Red = 0;
    rgbColor.Green = 0;
    rgbColor.Blue = 210;

    // Set the symbol of the GeoFeatureLayer's renderer to use the color initialized above
    ISimpleRenderer aoSimpleRenderer = (ISimpleRenderer)geoFeatureLayer.Renderer;
    ISimpleLineSymbol simpleLineSymbol = (ISimpleLineSymbol)aoSimpleRenderer.Symbol;
    simpleLineSymbol.Color = (IColor)rgbColor;

that works fine, but only for the SimpleRenderer. What I need is the ClassBreaksRenderer. Now the problem is, that if I create a new Renderer (SimpleRenderer or ClassBreaksRenderer), set the properties and apply it to the geo feature class (geoFeatureLayer.Renderer = cRenderer as IFeatureRenderer) the map is empty and no features are rendered.

Someone has an idea?


texts:
I want to move the zoom to an extent.

I have this var:

var extentGeo = esri.geometry.webMercatorToGeographic(map.extent);
console.log('extentGeo', extentGeo);
// { cache: undefined
// spatialReference: {wkid: 4326}
// type: "extent"
// xmax: -70.60487747064096
// xmin: -70.67328452935668
// ymax: -33.42183884025852
// ymin: -33.48435589231657 }

Is there a way to move the map using those extent.xmin, extent.ymin, extent.xmax, and extent.ymax ?


texts: Converting from ArcGIS layout into Google Earth as kmz with legend?

texts:
How to display an ArcGIS legend in kmz or kml file?

I am trying to do so in conversion tool of ArcGIS but while exporting it to Google Earth it is not showing the legend.


texts: Add Field giving ERROR 000464: Cannot get exclusive schema lock?

texts:
I'm doing an interface using ArcGIS and Python scripts that run at Python window.
For each time that my sript run, I need to update information, so I think that I need to delete originals tables and create news.

When I try to do this with ExcelToTable  it gives me an error to add field:

arcpy.env.workspace = "F:\Otim\inter"
arcpy.Delete_management(c_re) #delete layer
c_re=arcpy.ExcelToTable_conversion("02_Rea.xlsx", "cg1_re_otim.dbf", "Rea_L")
#arcpy.env.workspace = "F:ArcGIS\Trab\Model\Result"
arcpy.AddField_management(c_re, "Name_Ba", "Text", 50, "", "", "refcode", "NULLABLE", "REQUIRED")
arcpy.AddField_management(c_re, "Longitude", "Float", "", "", "9", "", "NON_NULLABLE", "NON_REQUIRED", "")
arcpy.AddField_management(c_re, "Latitude", "Float", "", "", "9", "", "NON_NULLABLE", "NON_REQUIRED", "")

  ERROR 000464: Cannot get exclusive schema lock.  Either being edited
  or in use by another application. Failed to execute (AddField).

How can I solve this problem?


texts:
I would like to simply rename a layer inside CartoDB.
I see buttons for the SQL Query, CSS code, and some HTML... but I don't know where to put my code.
setLayerName('Test')


texts: copying features to same layer or different layer (features didn't copy) QGIS MAPINFO TABs?

texts:
We are trying to move from Mapinfo to QGIS. We have some tables with .TAB format. What we normally do is copy feature from one layer into another layer (to the same layer). But it does not work in QGIS. Geometry type of feature are Line (WKB type: "LineString"). Properties for both layers are identical.

This is the link for two tables we want to copy
https://drive.google.com/drive/u/0/folders/1dDgYWhYe9KcTeaElJnKn43O9wlYD5SLI

Can someone tell me how to do this?
I want to copy all the feature from cable to cable2 table.


texts: Parameter to script tool gives Invalid expression ERROR 000358?

texts:
I have created python script tool in arcgis 10.1 as below.
I am trying to take input parameter and passing to select statement.

There is error when i execute this tool.

import arcpy

from arcpy import env

env.workspace = "Database Connections\sc_dev.sde"

name = arcpy.GetParameterAsText(0)

arcpy.SelectLayerByAttribute_management ("SC.STATES", "NEW_SELECTION", " \"STATE_NAME\" = name")

The error is Failed to execute....Invalid expression..ERROR 000358


texts:
I need to create a new gis application with OpenLayers and PostGIS.

I think that Geotools is a good choice for the Layer between these two.
I thought about embedding geotools in a web application

Now I need same samples to understand how geotools can be used together with OpenLayers.

Everything I found so far is good for developing a Java client application and of course for receiving and processing geo objects, but nothing that could help to build a web application.

Where should I continue to look?


texts: Iterating clip over multiple layers using QGIS?

texts:
I have multiple vector layers and a fishnet grid. I would like to clip all the layers based on all features in this fishnet grid.

I can use the iterate over layer function in the clip tool, but cannot select multiple layers. I would like this to be automatic.

So my goal is to get multiple input layers, but keep the iterate clip layer the same.

I have tried applying the knowledge in Iterating over map layers in QGIS python?, but cannot seem to understand or get it to work to iterate using the Clip tool.


texts: Building custom tooltip with hyperlink using QGIS

texts:
I am working on a project for publishing old scanned maps in jpeg at the scale of 1: 500e and converted into pdf. The goal of the project is to digitize the polygons representing the street small islands associated with the maps in order to geolocate them via a hyperlink from a tooltip when the user moves with the mouse cursor in the map canvas with QGIS.

The PDF files of the maps are located in a remote directory of our server and I would like to be able to create a tooltip in which I can display the attribute informations of polygon layers as well as the hyperlink that allows me to display the corresponding plan to my zone (location of the polygons of the small islands).

How could I achieve this in a simple way in QGIS?

All your suggestions and comments are welcome!


texts: Using the Offline Editing plugin in QGIS my project becomes 'stuck' in offline mode?

texts:
When using the offline editing plugin in QGIS to make a remotely accessed PostGIS layer editable offline, I have found that if there is a problem it is possible for my project to become stuck in offline mode. This means that the project only has the option to synchronize and the project can no longer be used to open a PostGIS layer and convert it to an offline project.

This has happened when the offline SQLite file has been moved, overwritten, created on network drive that is no longer accessible or the layer has been removed from the project while it is still offline.

I have also encountered problems when the synchronize has not worked, perhaps QGIS crashes or the connection to the database drops out. (Problems can be caused trying to Sync unsaved layer edits, or having other QGIS windows open that if interacting with them when syncing causes QGIS to hang).

Is there a way to 'repair' the project so that it's offline status can be reset?