texts:
I'm mapping on QGIS, I have a vector layer and his CRS is on EPSG:4326 - WGS 84, but also I wanted to put a satellite image, so when I do it, the image use WGS 84 / Pseudo Mercator.
The problem is, this CRS shows me differents degrees. And If I set project CRS from the vector layer, that will cut the image. Thus, I want the map with the coordinates of the vector layer.
Can I do it? How?


texts: Performing Spatial join with polygons using st_join in R?

texts:
I want to join 2 partly overlapping polygons in R the same way I'm doing it in ArcGIS Desktop . In ArcGIS Desktop I am doing it like this: Analysis toolbox - Overlay toolset - Spatial Join. The rows are matched based on their spatial locations. I want to "merge" all the columns of polygon1 with columns of polygon2 but only for rows (polygons) that exist in polygon1.

Unfortunately I can't give a reproducible example right now. This is what I did:

result_r <- st_join(x = polygon1, y = polygon2, join = st_intersects, left = TRUE, largest = TRUE)

That produced the warning: attribute variables are assumed to be spatially constant throughout all geometries. I used largest = TRUE because otherwise there would have been too many duplicated rows.

The ArcGIS Desktop result has the same number of rows as the R result. But when I check the number of rows for a certain condition, then the R result has more rows than ArcGIS Desktop :

result_arcgis <- readOGR(dsn = "path", layer = "result_arcgis")
result_arcgis <- st_as_sf(result_arcgis)
nrow(result_arcgis[which(!is.na(result_arcgis$colum_x)),])
# 30 rows
nrow(result_r[which(!is.na(result_r$colum_x)),])
# 51 rows

Obviously st_join isn't the same as the described ArcGIS Desktop procedure.

What is the correct R function to reproduce the ArcGIS Desktop procedure?


texts: How to unselect features in ONLY the current layer in QGIS?

texts:
How does one unselect features only from a certain layer?

I know there's the Deselect features from ALL layers button. But that's extreme. Is there a simple way to unselect features ONLY from (say) the current layer?

I can only think to ctrl+click everything, but this hardly seems effective.


texts: Delete files that were used in python code

texts:
I have a code that creates some tif images in .qgis2 folder, but I want, in the end, to delete these files. How do I do that?

I try with os.remove but it gives me the error32, that said that the process can´t access to file because this is used by another process.


texts: Openlayers - Different zoom levels with a WMS layer and only one resolution

texts:
I am trying to implement a homemade WMS with Openlayers. I only have one "level" of images with one resolution and I would like to add different zoom lvls with those images.

Is this possible just with Openlayers ?


texts:
I've been reading about using WMS servers to provide access to shapefile data. When a server does support shapefiles over WMS, what format does the client download the data in?

Will it typically be an image format (png, jpeg, etc) like it is for CADRG, or is there a way to extract the vector data, or will the shapefile data (geometries, etc) typically be part of the metadata available on the server?


texts:
I'm trying to download and install the Plugin "QgisGeoJSONExportPlugin" I downloaded from the "QGIS Python Plugins Repository", since I couldn't find on the "manage and install plugins" window.

So I downloaded it and added the folder in the "plugins" in the C: Directory
C:\Users\jan\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins

I didn't get an error message but the folders font color is green which I have never seen before. And even after adding the plugin in my Directory it's not in my list of installed plugins.

has anyone any idea what the problem might be?


texts: How to Set a MXD to Open With the Data Always 'Paused'? ArcMap 10.3.1?

texts:
I have several MXDs with large datasets (some with over a million features). I would like to open the MXD with the data in a 'paused' state, rather than having the MXD open and draw every feature; since this takes way to long. If I open the MXD, then try to 'pause' the drawing of the features sometimes this works immediately and the data pauses drawing, other times this attempt takes way too long.

I have to keep the layer in the MXD always on, since later the MXD will be used to create a mapping service.

Is there a setting where I can set the MXD to always open as 'paused'?

This is what I would like my MXD to look like when I open.

I am using ArcMap 10.3.1.


texts: Assigning known georeference to PNG using GDAL

texts:
I have a little experience with GDAL using the command line. However I am not sure where to start with this issue as it deals with adding projection to an image and I have never done it from scratch.

I have a PNG image (radar/rain data) that I need to convert to geoTIF. I was told what the projection and bounding coordinates are, but I am not sure if I have enough information or even if I am doing it right.

Here is the information I was given :

Longitude: -102.328, -93.7296
Latitude: 26.0468, 33.3013
Projection: AzimuthalEquidistant.

For additional clues, I know the image was generated by in python with Py-Art (weather radar tools). I found information on projection output with this tool at  this link, which states about image generation :

  The GeoTIFF can be the standard Azimuthal Equidistant projection used
      in Py-ART, or a lat/lon projection on a WGS84 sphere.

So, I know the image was generated using the first option, but is this enough information to properly georeference my image? I see projection code in that link but it seems it is for their latter option (non standard Azimuthal Equidistant Projection). My end goal is to get the PNG projected to EPSG:3857 for Mapbox, but I am not sure how to get there. Here is what I have tried :

So, from what I read online, Azimuthal Equidistant Projection is also known as EPSG:4326. Perhaps what I am trying is too simple - but here was my command using what I believe are the coordinates and projection :

gdal_translate -of Gtiff -a_ullr -102.328 33.3013 -93.7296 26.0468 -a_srs EPSG:4326 image.png image2.tif

I then used this command to get to EPSG:3857 :

gdalwarp -t_srs EPSG:3857 image2.tif image2-projected.tif

I verified that the final image was wrong (image2-projected.tif).

What am I doing wrong? I see from other examples online that there are "proj4" strings and SRS data, but I don't exactly understand what needs to be added and how to add it (if anything.)

Last, here is the image itself for completeness, and  this is a public google drive link to get the actual PNG :


texts: (.NET) Get DTED Data (Elevation) Between Two Points

texts:
I want to get DTED data (elevation) between two points (on a line).
Is there any built-in function for this purpose or I should implement it on my own?

What I want is similar to Linear Line of Sight (LLOS), but it should be on a direct line. I want to check if there is collision between two points or not.

Thank you

Sincerely


texts: Custom sort index for Data Driven Pages?

texts:
I am creating maps of underground utilities using DDP with a grid layers as my index layer. Each grid shows a segment of the run. The pages must be generated in an order that maintains the continuous run.

My problem is that I don't have a good sort field because none of the grid fields adequately align with the continuity of the run. For example, page 1 should be where the run starts and page 17 would be where it ends. But I don't have a sorting index for this  unless I create a page index field and manually assign values.

The problem is that I have over 45 MXDs to data drive, each with its own grid package, and some grid packages containing up to 15 grids apiece.

This means I am having to manually sort roughly 200 pages in total which is going to take way more time than I have.

Is there a way to rearrange the page order without having to rely on the index layers attributes? At first I accomplished this by simply reordering the PDF pages in Adobe, however this is not a good option because then my dynamic text page numbers are off (for example page 1 of 17 might be followed by page 8 of 17 followed by page 6 of 17 etc.


texts: Best Resampling Method for DEM reprojection in QGis

texts:
I need to reproject a DEM (aster file) in order to have the same CRS as all my other project data - I have different SHP for the same DEM's location but they have a different CRS.
Which resampling method would you recommend? I am using QGis.


texts: Updating cluster distance in OpenLayers 3?

texts:
I have a vector layer with clustering in OL3.4:

cluster = new ol.layer.Vector({
          source: new ol.source.Cluster({
            distance: 20,
            source: some_geojson
          }),
          style: some_style,
});

I want to change the cluster 'distance' property, after some event, particularly:

map.getView().on('change:resolution',function(){
    // change cluster distance here (and refresh it)
}

What is the best approach to do this?


texts: I'm looking for free (non-commercial use) sources of radar (SAR) satellite imagery. Does anyone know of sources?


texts: Meaning of Spatial Join ERROR 000840: The value is not a Feature Layer?

texts:
I'm using ArcGIS (ArcMap) 10.2 to determine the relative location of points of land to a regulated area. I downloaded a layer (.lyr) file that displays the regulated area, and am able to display that information with my data. However, I am unable to use the data in (or referenced by) the downloaded layer... i.e. I cannot use Spatial Join ("ERROR 000840: The value is not a Feature Layer").

To try to determine what was happening, I attempted to "Open Attribute Table" to look at the specific features, but the "Open Attribute Table" command was not visible.  That is, when I right-click the layer name in the table of contents, I don't have the normal options I do for features.

Thus, it seem I am unable to use this layer as I would a normal feature set. Why is this, and how can I still use this data?

FYI: I'm a relative novice at ArcGIS, so perhaps there are fundamental issues I'm not getting?

A screenshot of right before the error, and another right after the error


texts: Is anyone using JSONP in Web mapping applications to get around the browser same origin policy?

texts:
Javascript based Web mapping apps that use ajax requests to get data from a server are limited by the Same Origin Policy to only request the data from a server in the same domain (subdomain and port).

One method commonly used to get around this is to have a script that acts as a proxy located within the same domain.  The proxy sends the request to the server outside of the domain and then forwards the results back to the requesting browser.

I have heard of some non-spatial Web apps using JSONP "Javascript Object Notation with Padding" to get around the same origin policy.  I haven't seen it in use in any Web mapping apps.

Can anyone identify an example?  What are the pros and cons?


texts:
I want to save a KML point file with attributes and descriptions from a generic point vector layer, and I need to do this at the end of a processing chain inside QGIS graphical modeler.
Ideally, I would have the option to choose the icons of the KML points.

From here, I've started to use the "Convert Format" algorithm from GDAL asuming that this would be the simplest option.

I've wrote the following script to do so (but I'm missing how to call the -dsco NameField=RegionName option to have the descriptions correctly labeled)

##Convert=name
##Layer=vector

layer = processing.getObject(Layer)
output_path = "C:/Users/jmgc/Desktop/testes_script/" + layer.name()

processing.runalg("gdalogr:convertformat", Layer, 15, '', output_path)

Adding the -dsco NameField=RegionName gives me an error:

##Convert=name
##Layer=vector

layer = processing.getObject(Layer)
output_path = "C:/Users/jmgc/Desktop/testes_script/" + layer.name()

processing.runalg("gdalogr:convertformat", Layer, 15, -dsco NameField=RegionName, output_path)

Any thoughts on how to solve this on my prescriptive or following another one would be helpful.


texts:
Our application uses OpenLayers+PostGIS to visualize and process millions of spatial data files (GeoJSON).
We have a requirement to be able to export the processed final data as PDF/GeoPDF Vector maps.
Every solution, including the official OL example, first converts the OL canvas to a raster, then converts the image to PDF using something like jspdf. The raster conversion completely destroys the vector data, making the final PDF unusable for professional purposes (lines are pixelated, labels are blurred) (example below)

Is there a method that *preserves the Vector points/lines/labels in the map by skipping the raster conversion? Example: maps exported from AutoCAD/ArcGIS/QGIS have perfect, infinitely zoomable vector pdfs.
I've looked far and low for this, but couldn't find an answer, except for using 'gdal_translate', which seems cumbersome and poorly documented for programmatic usage.


texts:
from qgis.analysis import  QgsRasterCalculator, QgsRasterCalculatorEntry

raster = '/home/francis/QGIS/epelwrf/temp/temp0000-postwrf_d01_20180720_0600_f00000.gr2'
bohLayer = QgsRasterLayer(raster)
entries = [ raster ]

#Band
boh1 = QgsRasterCalculatorEntry()
boh1.ref = 'boh@1'
boh1.raster = bohLayer
boh1.bandNumber = 1
entries.append ( boh1 )

expression = '( ' + entries[0].ref + ' ) - 275.15 '
print expression

calc = QgsRasterCalculator( expression, "/home/francis/QGIS/tiffout/1.tif", 'GTiff', bohLayer.extent(), bohLayer.width(), bohLayer.height(), entries )

calc.processCalculation()

I'm stuck in my project. trying to convert this gr2 file (that contains tempsfc which is in kelvin) into Celsius.


texts: Build geonetwork: FileNotFoundException rebel.xml

texts:
Just tried to build geonetwork on my ubuntu machine, following this and this tutorial.

Unfortunately I ran into a problem, I don't know how to fix.
I managed to get the geonetwork clone from github and did

  mvn package install -DskipTests

This finished the build of common utils with an error, after having done all the downloads. Here's the stacktrace from running with -e switch.

[INFO] Reactor Summary:
[INFO]
[INFO] GeoNetwork opensource ............................. SUCCESS [0.341s]
[INFO] common utils ...................................... FAILURE [0.756s]
[INFO] Caching xslt module ............................... SKIPPED

...

[ERROR] Failed to execute goal org.zeroturnaround:jrebel-maven-plugin:1.1.5:generate (generate-rebel-xml) on project common: Failed writing rebel.xml: /srv/git/repos/geonetwork/core-geonetwork/common/target/classes/rebel.xml (No such file or directory) -> [Help 1]
org.apache.maven.lifecycle.LifecycleExecutionException: Failed to execute goal org.zeroturnaround:jrebel-maven-plugin:1.1.5:generate (generate-rebel-xml) on project common: Failed writing rebel.xml
        at org.apache.maven.lifecycle.internal.MojoExecutor.execute(MojoExecutor.java:217)
        at org.apache.maven.lifecycle.internal.MojoExecutor.execute(MojoExecutor.java:153)
        at org.apache.maven.lifecycle.internal.MojoExecutor.execute(MojoExecutor.java:145)
        at org.apache.maven.lifecycle.internal.LifecycleModuleBuilder.buildProject(LifecycleModuleBuilder.java:84)
        at org.apache.maven.lifecycle.internal.LifecycleModuleBuilder.buildProject(LifecycleModuleBuilder.java:59)
        at org.apache.maven.lifecycle.internal.LifecycleStarter.singleThreadedBuild(LifecycleStarter.java:183)
        at org.apache.maven.lifecycle.internal.LifecycleStarter.execute(LifecycleStarter.java:161)
        at org.apache.maven.DefaultMaven.doExecute(DefaultMaven.java:320)
        at org.apache.maven.DefaultMaven.execute(DefaultMaven.java:156)
        at org.apache.maven.cli.MavenCli.execute(MavenCli.java:537)
        at org.apache.maven.cli.MavenCli.doMain(MavenCli.java:196)
        at org.apache.maven.cli.MavenCli.main(MavenCli.java:141)
        at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
        at sun.reflect.NativeMethodAccessorImpl.invoke(NativeMethodAccessorImpl.java:57)
        at sun.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43)
        at java.lang.reflect.Method.invoke(Method.java:606)
        at org.codehaus.plexus.classworlds.launcher.Launcher.launchEnhanced(Launcher.java:289)
        at org.codehaus.plexus.classworlds.launcher.Launcher.launch(Launcher.java:229)
        at org.codehaus.plexus.classworlds.launcher.Launcher.mainWithExitCode(Launcher.java:415)
        at org.codehaus.plexus.classworlds.launcher.Launcher.main(Launcher.java:356)
Caused by: org.apache.maven.plugin.MojoExecutionException: Failed writing rebel.xml
        at org.zeroturnaround.javarebel.maven.GenerateRebelMojo.execute(GenerateRebelMojo.java:224)
        at org.apache.maven.plugin.DefaultBuildPluginManager.executeMojo(DefaultBuildPluginManager.java:101)
        at org.apache.maven.lifecycle.internal.MojoExecutor.execute(MojoExecutor.java:209)
        ... 19 more
Caused by: java.io.FileNotFoundException: /srv/git/repos/geonetwork/core-geonetwork/common/target/classes/rebel.xml (No such file or directory)
        at java.io.FileOutputStream.open(Native Method)
        at java.io.FileOutputStream.<init>(FileOutputStream.java:221)
        at java.io.FileOutputStream.<init>(FileOutputStream.java:171)
        at org.zeroturnaround.javarebel.maven.GenerateRebelMojo.execute(GenerateRebelMojo.java:220)

Any ideas how to create the missing rebel.xml or work around it? Tried mvn clean and mvn clean install as well, but without success.

UPDATE:
I updated maven version from 3.0.4 to 3.3.9, which was helpful according to this post, but without success, I still get the same error


texts: What is the name of the 3 word address system?

texts:
I remember this system(website) from a couple of years ago, but I don't remember the name.  But the general idea is that you could identify any location on the earth with 3 words.  IIRC, there was no actual data available but there was a paid web API.

What was this website?


texts: Custom Python Expression Functions (QGIS Tutorial)

texts:
I'm having trouble finishing the QGIS tutorial for Using Custom Python Expression Functions: http://www.qgistutorials.com/de/docs/custom_python_functions.html

I get to the point where I write the script in the function editor, but when I want to use the function in an expression it looks different and doesn't work.

Screenshot 1 is what it looks like in the tutorial, screenshot 2 is my result when I double click on the custom function. The code is the same, I copy pasted it into the functions editor. When I write the exact same expression as in the tutorial it says 'expression invalid' (screenshot 3).


texts: How to get vector feature by mouse location?

texts:
I'm writing a photo grouping plugin. My current problem can be broken up into "end goal" and "baby steps", because I think that's what I want.  But if there's a way to get to the end by a different means, that's what I really need.

End goal: given 2 different QGIS layers, 'nodes' and 'sections', and given a bunch of groups of photos, mouse over a 'node' > click to assign photos to node.  mouse over a 'section' > click to assign.  Repeat indefinitely, not necessarily in order.  Without toggling active layer back and forth and making selection(s) first.

Baby Step #1: find out what feature(s) on which vector layer(s) is under the mouse position.

I found the "Value Tool" plugin last night, read through the code, and it reports values from multiple raster layers based on the mouse position.  I was super excited! Unfortunately, it seems that QgsVectorDataProvider doesn't have anything analogous to the QgsRasterDataProvider's identify(position) method, which is how that plugin works.

So I've been playing around with QgsVectorDataProvider.select(), but can't get it working.  Probably doing something stupid.  This post:
https://gis.stackexchange.com/a/22919/6690
makes it sound like it's possible.

Can anyone tell me what I'm doing wrong?

Steps (in Python console for now):

  load vector layer in QGIS
  iface = qgis.utils.iface; canvas = iface.mapCanvas(); layer = canvas.layer(0)
  layer.dataProvider().attributeIndexes() >> [0,1,2,3,4,5]
  A = layer.dataProvider.select((0,1,2,3,4,5),canvas.extent())
  checks: A >> nothing;   type(A) >> 'type 'NoneType''

Once I can get a selection by position, I think I can get the list of layers I want to check, iterate through them with position, check & flag on overlaps, and we're off...


texts: slow import via osm2pgsql to postgresql database

texts:
I try to import berlin.osm (33Mb) into the postgres database 9.1 via osm2pgsql. My desktop is mem 2Gb,proc 2Ghz,hdd 60Gb, RHEL6u2x64. This process is so slowly in particular when I'm going over pending ways. I see this:

osm2pgsql SVN version 0.80.0 (32bit id space)

osm2pgsql -M -d berlin -s -S "/usr/share/osm2pgsql/default.style" -U "postgres" -W -H "localhost" -P 5432 "/opt/maps/berlin/berlin.osm.bz2"

Processing: Node(1644k 49.8k/s) Way(304k 12.17k/s) Relation(6337 52.37/s)  parse time: 179s

Node stats: total(1644720), max(1968785138) in 33s
Way stats: total(304296), max(186177757) in 25s
Relation stats: total(6337), max(2504060) in 121s
Committing transaction for planet_osm_point
Committing transaction for planet_osm_line
Committing transaction for planet_osm_polygon
Committing transaction for planet_osm_roads

Going over pending ways

Using 1 helper-processes
processing way (10k) at 0.11k/s
processing way (195k) at 0.11k/s
All child processes exited

Pending ways took 1719s at a rate of 114.01/s
Committing transaction for planet_osm_point
Committing transaction for planet_osm_line
Committing transaction for planet_osm_polygon
Committing transaction for planet_osm_roads

Going over pending relations

Using 1 helper-processes

All child processes exited

Committing transaction for planet_osm_point
Committing transaction for planet_osm_line
Committing transaction for planet_osm_polygon
Committing transaction for planet_osm_roads
Sorting data and creating indexes for planet_osm_point
Sorting data and creating indexes for planet_osm_line
Sorting data and creating indexes for planet_osm_polygon
node cache: stored: 1644720(100.00%), storage efficiency: 66.93% (dense blocks: 60, sparse nodes: 1597241), hit rate: 100.00%
Sorting data and creating indexes for planet_osm_roads
Stopping table: planet_osm_ways
Building index on table: planet_osm_ways (fastupdate=off)
Stopping table: planet_osm_nodes
Stopping table: planet_osm_rels
Building index on table: planet_osm_rels (fastupdate=off)
Stopped table: planet_osm_nodes in 0s
Stopped table: planet_osm_rels in 1s
Indexes on  planet_osm_roads created  in 5s
Completed planet_osm_roads
Indexes on  planet_osm_line created  in 21s
Completed planet_osm_line
Indexes on  planet_osm_point created  in 26s
Completed planet_osm_point
Stopped table: planet_osm_ways in 31s
Indexes on  planet_osm_polygon created  in 31s
Completed planet_osm_polygon

Osm2pgsql took 1931s overall

I used the next parameters in my postgresql.conf :

shared_buffers=512Mb
temp_buffers=16Mb
work_mem=64Mb
maintenance_work_mem=128Mb

When I checked it on the Windows desktop it acceptable worked, but if I started it on my linux system I have a problems with a import perfomance. Could someone give me advice what I should do that can fix it?


texts:
I am trying to add kml to openlayers and i am getting error in firebug

  a is undefined: Line 92 openlayers.js

here is my code :

<!DOCTYPE html>

<meta http-equiv="Content-Type" content="text/html; charset=utf-8">

<title>OpenLayers Basic Single WMS Example</title>


texts: Is it possible to set the extents in which labels appear in Leaflet?

texts:
I'm fairly new to leaflet and web mapping in general. I was curious if there is a way to set labels to only appear at a certain range of extents. The point data layer I have clusters too much preventing any legible labels. It also makes the map look really cluttered.

I have been using the leaflet.label plugin but will use others if needed.

Using Leaflet 0.7


texts: Associating Point Layer, Annotation Group and Line Layer in Map Legend?

texts:
I have an annotation group containing charts associated with a point layer in an ArcGIS 10.2 for Desktop layout. There are a number of line types on the charts.

How can I add the chart line symbology to the legend so it is displayed when the point layer is displayed?


texts: Updating (ArcMap 10) layer objects with objects from another layer?

texts:
I have layer "A" and layer "B".
Layer "A" has polygon objects and unique ID for every object in attribute table.
Layer "B" has some polygon objects from layer "A", but those polygons have been changed, unique ID is not changed.
How can I update layer "A" polygon objects from layer "B". As I say, those objects share same unique ID.


texts: QGIS Raster clipper not working on reprojected layers?

texts:
I try to clip a hillshade raster and I get this error. The hillshade was the result of a reprojected DEM from WGS84 to WGS84UTM. When I to clip the hillshade or even the reprojected DEM, I get the error. Clipping the original WGS84 DEM works fine.

Any help?


texts: Character marker symbology based on a field

texts:
I'm a bit unsure of how to express this question without being too unspecific but I hope you understand what I'm getting at...

I have a symbol library as a .ttf with geotechnical symbols. I would like the symbology to pick a symbol from the library without manually adding all the different symbols using hex or unicode values in a field.

I have found lots of information regarding applying color based on hex etc. However I have not been able to find the above example/implementation.

I know for a fact it is a possibility, at least in a web-based GIS environment. And library is based on the sum of the number of bits for each symbol translate to a hex value. I don't know if that makes any sense.


texts:
I can't find anything on google or the QGIS changelog on this issue:

We're using a style for points which contains a point, a circle and an ellipse.

Point and circle (buffer($geometry,5,80)) work fine for everyone, but the ellipse (make_ellipse($geometry,15*2,15*1.5,135,80)) don't work for a colleague of mine, neither in QGIS 3.10 nor in 3.4.

Maybe the ellipse has different software dependencies or needs a setting enabled which I forgot about? Has someone an idea?


texts: Use mapview in function to open viewer, but also do other things afterwards

texts:
The question may sound a little strange, but what I want to do is the following:
library(sp)
library(sf)
data(meuse)
meuse = st_as_sf(meuse, coords = c("x", "y"))

plotit = function(shape){
   mapview(shape)
   print("mapview is awesome")
}

plotit(meuse)

It will print "mapview is awesome", but the viewer won't open. I can remove the print statement and it will open. What am I missing here?:)


texts: Multi-labelling with different colors in QGIS 3.10

texts:
I would like to create a multi-field label based on an expression where the different fields will have different colors.
Using an icon font, this will allow me to get this kind of result (done with ArcGIS):

I tried using similar expressions in QGIS, changing format to get it "QGIS compatible", but it doesn't work (ArcGIS code: <CLR red='130' green='130' blue='130'> & [RURAL] & [URBAN] & [CAMP] & </CLR>.

It seems it was not possible in previous version, I just tried with 3.10... but still the same.

However it's possible in the print composer to use html tags to define colors. So am I missing something to make it work on data label also ?


texts:
I've just installed Geonetwork.  I'm running it via localhost:8080, but it is very slow.  Signing in can take several minutes or it disconnects because of lack of activity.  That's about as far as I've got because it's so slow.
I pinged my localhost and it responds instantly.  What else can I do to diagnose and fix this issue?

Windows7 64 bit, 32 GB RAM, Xeon processor 8 cores.
Java version 1.8.0_261
Tomcat 9.0.38
Geonetwork 3.10.3

I just got an error saying something about Java Memory Heap - I didn't write it down and can't reproduce it at the moment but it gives me something to Google at least.

I seem to have progress - I added the following line to setenv.bat:
set "CATALINA_OPTS=-Xms512M -Xmx1024M"


texts: Geotransformation, from sat coordinates to Earth coordiantes

texts:
I'm using this code: https://github.com/mcidasv/mcidasv/blob/master/ucar/unidata/geoloc/projection/sat/GEOSTransform.java
to retrieve longitude and latitude in Earth system from column, line coordinates in Satellite system.
This is my code:

public class Main {

private final static double COFF=1856;
private final static double LOFF=1856;
private final static double CFAC = -781648343;
private final static double LFAC = -781648343;

public static void main(String args[]){

    GEOSTransform geo = new GEOSTransform(0, GEOSTransform.GEOS, GEOSTransform.WGS84);

    /*
    double[] lonlat = geo.FGFtoEarth(1, 34, CFAC, COFF, LFAC, LOFF);
    System.out.println("LONLAT "+lonlat[0]+", "+lonlat[1])
    */
    for(int x=1; x<1000;x+=1)
        for(int y=1;y<1000;y+=1){
            double[] lonlat = geo.FGFtoEarth(x, y, CFAC, COFF, LFAC, LOFF);
            System.out.println("LONLAT "+lonlat[0]+", "+lonlat[1]+": xy "+x+", "+y);
        }
    }
}

and I get NaN values from FGFtoEarth. Why? For that method, the doc says "Transform fractional FGF coordinates to (longitude, latitude)." What are fraction FGF coordinates? I think there is something wrong with the passed in coordinates.
Note that const members I think are correct, but I'm not sure if CFAC and LFAC are negative like above, or positive. Anyway in both case the result is the same, NaN values.


texts:
In the win10 version of QGIS, you can install extra Python packages with pip commands directed at the QGIS python 3 environment. For example, I've got access to scikit-learn form my QGIS Python console.
What is the equivalent process for a mac?
The sticking point seems to be getting the terminal to point to QGIS's Python 3 environment. In Windows it's just a py3_env command. No luck on MacOS.


texts: Calculating Grid Convergence for Lambert Conformal Conic

texts:
What is the formula to calculate Grid Convergence angle (angle between true north and grid north) for Lambert Conformal Conic projection for a given location.

Searching online, I have been able to find the formula for calculating Grid Convergence for UTM grid only.


texts: How do I install QGIS Server and QWC2 via Gitpod?

texts:
Maybe this question is better suited to the Web Apps SE.

I want to create a self-hosted instance of QGIS Server and the QGIS Web Client 2, for times when I'm not on a computer with admin access (such as a library). Ever since C9.io shut down I've been using Gitpod for free self-hosting. However, unlike C9, Gitpod doesn't give you root access (instead everything is through Dockerfiles)

I found two Dockerfiles, one to install QGIS, and one to install the QWC2. However, I have no earthly clue how to "put it all together", especially since I'm not sure if Gitpod even supports docker-compose.

Or I could just use the Heroku CLI Dockerfile on Gitpod and host QGIS via Heroku...


texts:
I'm trying to calculate the total number of counts a raster contains using python script. When I ran the script, I didn't get any error. However, I also didn't get the number of counts for each row in the attribute table.

Code:

import arcpy

#To overwrite output
arcpy.env.overwriteOutput = True

#Set environment settings
arcpy.env.workspace = "C:/Subhasis/Test/raster-read"

#Build Raster Attribute Table
arcpy.BuildRasterAttributeTable_management("sti-9", "Overwrite")

# loop through attribute table
Rows = arcpy.SearchCursor( "sti-9" )

for row in Rows:
    count = row.getValue( "COUNT" )
    print count


texts:
I have a pivot table that has 11 ShipTypeGroup fields with "Long" field properties that can't be divided by 12, due to decimal places. Therefore I am trying to create a new field, copy attributes across from each ShipTypeGroup and recalculate. Each ShipTypeGroup must be recalculated into a field that has the same ShipTypeGroup. So for example, adding a column called Ave_STG_1 and then calculating the figures from ShipTypeGroup1 field and so on for all 10 types. However when it throws up an error or doesn't calculate the correct fields properly. I'm not sure if I have passed the newly created field into the "Calculate Field" part of the script properly, or do I add an if/else statement.

import arcpy

pt = r"C:\AIS\Density_Grid_test.gdb\PivotTable"
flist = arcpy.ListFields(pt, "*ShipTypeGroup*")
for f in flist:
    arcpy.AddField_management(pt, f.name.replace("ShipTypeGroup","Ave_STG"),"DOUBLE")
    newflist = arcpy.ListFields(pt, "*Ave_STG*")
    for nf in newflist:
        arcpy.CalculateField_management(pt, nf.name, "!{}!/12".format(f.name), "PYTHON_9.3")
        print f.name + " changed into " + nf.name + " and recalculated"


texts:
I am working on a web map using Leaflet and am using purchased satellite imagery as the base map. We want to limit the extent of the high resolution image to the area of interest, so I have resampled the image from 2.4 to 20 meters, and cut out the relevant area from the 2.4 meter image (Using QGIS). I created tiles with gdal2tiles for both of these images, and have overlayed them in my map. Here is a screenshot of the result:

On some of the edges where the high resolution imagery ends, a black line appears. Looking at the tiles in GIMP there is no such border, and whether or not there is a black line is zoom dependent, so I suppose it's due to the rendering in the browser.

I tried to solve the problem by doing the resampling again at 21.6 meters and lining up the rasters so they would "nest", but the random black lines still appear.

Now I'm thinking of ways to blend the edge of the fine resolution image, or of the generated tiles directly. I could graduate the resolution so there isn't a "drop" from 2.4 to 21.6 meters, or use a transparency gradient. I'm not quite sure how to go about the first idea, but tried out the second in GIMP. I blended the edges of one tile and was happy with the result in the map. I don't want to do the tedious work of manually editing each tile, though. Blending the image before cutting into tiles doesn't work because of the resampling of the resolution.

What is causing this black line? Is there a known solution for this?


texts: Dealing with NA-Values when plotting rasters with folium in python

texts:
I have some tiffs that I'd like to plot interactively in python using folium (or any other package, but I don't know of any). My problem is, that I have many na-values in the raster. At the moment I read in the raster with rasterio, read the data into a 2d-numpy-array, and then I write it back out as .png.
Finally I use folium with:

rasterlayer = folium.FeatureGroup(name = "raster")
rasterlayer.add_child(folium.raster_layers.ImageOverlay(
                        image=ndvi,
                        bounds=[[bottom, left], [top, right]],
                        interactive=True,
                        cross_origin=False,
                        zindex=1,
                        colormap=matplotlib.pyplt.cm.Greens
                            ))
m.add_child(rasterlayer)
m.add_child(folium.LayerControl())
m

But I get all the NA-Values printed in some colour. And I didn't find any solution to that problem. I had a look here:

https://stackoverflow.com/questions/50689068/display-raster-data-in-folium-handling-no-data-values

but also did not really help.
I also did set all NA-values in the "numpy"-array to 0, but this also was not the solution. The NA-Values always got plotted. Is there some, relatively easy way, to omit the plotting of na values. And more general: What is a good way of plotting raster with folium? It seems that in R it's so easy with packages like mapview or tmap and in python it feels way more complicated...


texts: Add row if attribute doesn't exist in feature class (ArcGIS/ArcPy)

texts:
I have 150 feature classes in a file geodatabase. Each consist of a few polygons with a different label class. I want to make a selection on Label "1E-06", but some feature classes do not have this label, only for instance label 1E-01 to 1E-05. Is it possible to make some if statement, so if label 1E-06 does not exist in feature class, add row with label = "1E-06". It doesnt have to have geometry. I want to merge all 1e-06 labels later. For statistics purposes I need the empty ones too.

EDIT:
I tried it myself (with help of answer below), but still get stuck. Right now, it can find all the shapefiles and put it in a database. It can search for rows with my expression, but I want to add a row for all features without the expression. I tried this with for row not in scursor: but this gives an error. How can I add a row in the features who do not have a row with the stated expression?

import arcpy
import os
count = 0
ws = "N:\Documents\Project\data"
dst654 = "N:\Documents\Project\Contour654.gdb"
dst800 = "N:\Documents\Project\Contour800.gdb"
expression654 = "Label = '1E-06 654'"
expression800 = "Label = '1E-06 800'"
for dirpath, dirnames, filenames in arcpy.da.Walk(ws,datatype="FeatureClass"):
    count += 1
    for file in filenames:
     filepath = os.path.join(dirpath,file)
     file_nice = file.replace('-', '_')
     file_nice = file_nice.replace(' ', '_')
     file_nice= file_nice.replace('.', '_')
     if '654' in file_nice:
         outpath1 = os.path.join(dst654,file_nice[:-4])
         arcpy.CopyFeatures_management(filepath,outpath1)
         arcpy.AddField_management(outpath1, "gid", "DOUBLE", "", "", "", "", "NULLABLE", "", "")
         arcpy.CalculateField_management(outpath1, "gid", count, "PYTHON", "")
         with arcpy.da.SearchCursor(outpath1, ['Shape', 'Label'],
                            where_clause=expression654) as scursor:
                                for row in scursor:
                                    print "yes"
                                    print outpath1
     else:
         outpath2 = os.path.join(dst800,file_nice[:-4])
         arcpy.CopyFeatures_management(filepath,outpath2)
         arcpy.AddField_management(outpath2, "gid", "DOUBLE", "", "", "", "", "NULLABLE", "", "")
         arcpy.CalculateField_management(outpath2, "gid", count, "PYTHON", "")
         with arcpy.da.SearchCursor(outpath2, ['Shape', 'Label'],
                            where_clause=expression800) as scursor:
                                for row in scursor:
                                    print "yes"
                                    print outpath2


texts:
To continue from this question about the format of a .csv file with network data.

I have two layers now: 1) layer with the coordinates (X,Y) of some nodes, and 2) layer that represents the lines connecting this coordinates (the .csv that contains that lines is in WKT with LINESTRING representing the lines). Both layers are displayed well within QGIS.

Now, I want to create a shapefile with these two layers. I tried to select the two layers and Save As> in ESRI Shapefile format. However, the new shapefile only contains the coordinates of the nodes, so I can only visualise the nodes, not the lines.

I don't see any option to do that. What do I need to do (if possible)? Is there any plug-in to do so? Any Python API that I can use?


texts: Filtering MSSQL layers in QGIS doesn't work

texts:
In my QGIS Plugin I am filtering layers with this Python code:

FeatureIterator = layer.getFeatures( QgsFeatureRequest().setFilterExpression ("Fieldname='SomeValue'") )
for feature in FeatureIterator :
    attrs = feature.attributes()
        print "Feature ID: " + str(feature.id()) + " " + str(attrs)

This works for most layer-types (shapefile, PostGIS, personal GDB,etc).
But when I tried it with MSSQL Spatial layers, no results are returned. The iterator is empty.

Edit: It works if I use the same filter clause (Fieldname='SomeValue') as a Feature subset at Layer-Properties. Then the layer shows only the intended subset.

I work with QGIS 2.2 and 2.4. MSSQL Server is 2008R2.
My test data was exactly the same in PostGIS and MSSQL, so I know that the filter expression is correct.

Could this be a bug?

Edit:
I was able to track the SQL queries that are executed on the Database. Instead of the filter expression that I intended, QGIS adds an Intersect clause and an ID filter to the query:

select .... from ... where
[shape].STIntersects([geometry]::STGeomFromText('POLYGON((0.00000000 0.00000000, 0.00000000 0.00000000, 0.00000000 0.00000000, 0.00000000 0.00000000, 0.00000000 0.00000000))',25832)) = 1 and [objectid] = 30681189080760422

The Coordinates of the intersect-rectangle are always 0 and the objectid-value doesn't exist.


texts: How to Extract data from OSM within a polygon

texts:
Although I am a Openstreetmap's regular user for quite some time, I am pretty new to GIS, so I assume this is an easy question: I would like to know how can I extract data (eg: all schools) from OSM within a certain area defined by a polygon (eg: a city boundary or a city's neighbourhood) using QGIS.

Extracting data

I have searched a bit about it and I found out (correct me if I am wrong) that there are different ways to perform queries on OSM, although I don't understand how/when do I have to use each one:

Osmosis
Overpass
XAPI

Using polygons

The problem is that as far (or as little) as I know all these methods rely on a square area, and I would like to know how can I do that within two types of polygons:

A city's boundary (this polygon exists in OSM, since is the one you see when searching for an specific city)
An arbitrary area, which would require to draw a new polygon using an external software (possibly QGIS? JOSM?).


texts: High resolution x,y dataset for Mollweide projection convert to latitude and longitude coordinates

texts:
I have a 2D array that's 144000 x 64000 that has some values that are irrelevant to the question. Let's say the array can be used like an X, Y coordinate system. I somehow need to get the Latitude and Longitude based from the Mollweide projection.

I took a look at the Wiki inverse transformation formula, but I believe I'm doing the math incorrectly. Let's say I have some X,Y coordinate (130000, 50000), how would I get the latitude and longitude from it?

Note: the actual array is composed of several arrays that are 4000 x 4000, with the number of arrays being 36 x 16.

My current math is based on this SO question:

R = (36 * 4000) / 2 / sqrt(2) = 50911.68825 radians
theta = arcsin(50000 / R / sqrt(2) = 43.9829 radians
latitude = arcsin[(2 * theta + sin(2 * theta)) / pi] = DOMAIN ERROR
    => since (2 * theta + sin(2 * theta)) / pi = 28.31853329
longitude = 130000 * pi / (2 * R * sqrt(2) * cos(theta)) = 3.941594775

So neither latitude nor longitude make ANY sense. Am I supposed to convert theta to something below 2 * Pi? What math am I doing incorrectly?


texts: Web site with: google map like or GIS, one million record database and a feature to run advanced search

texts:
What is the GPL way  - without or with little developing - to build a web site runnning on apache which has to feature: one million record database (MySQL; which may be editable), an advanced search on it with results (pins and/or images) to display on a google like map?

Example: http://livelovely.com/search


texts: Activate 'View Item Description' programmatically?

texts:
I have created an ArcMap add-in for browsing data layers and adding them to the TOC. I would like to add a button that when clicked will open the selected layers item description (before it is added to the TOC) however I can't seem to find any documentation on how to accomplish this.

Basically I want the same functionality as the 'View Item Description...' in the 'Data' menu when right clicking a layer in the TOC. So when the button is pressed the item description window is opened. Any one have any ideas I thought for sure there would be an easy way to do this?


texts: Determining average distance from points to line using ArcGIS Desktop?

texts:
I have a polyline road dataset and a study area (polygon) covering the the road dataset. I want to determine the average distance for every given point within the study area to the nearest road. i don't have a dataset with the points. Also I don't necessary want to create a new dataset. All I am interested is the average distance.

Is there a tool that can do that job in ArcGIS Desktop?


texts: Importing CSV file to create point layer

texts:
I have a CSV file with following fields:

Nodename
Nodeid
Latitude
Longtitude
Elevation

I am using Python console in QGIS to automate this. How should I write so that point with latitude, longtitude and elevation are maked as point layer?
I have tried with below code, but I  am getting error.

uri = "/home/priti/Desktop/MTP work/nodeinput.csv?\
type=csv&xField=Longtitude\
&yField=Lattitude\
&spatialIndex=no&subsetIndex=no&watchFile=no"

vlayer = QgsVectorLayer(uri, 'Nodes', "delimitedtext")


texts: setHeaderData not working in a QgsAttributeTableModel?

texts:
I want to display the content from a layer in a new view

I have so far a code like below (you should have 2 columns at least, then, just try in the QGIS Python console replacing mylayer_id)

from PyQt4.QtCore import Qt
from PyQt4.QtGui import QTableView
from qgis.gui import QgsAttributeTableModel

layer = QgsMapLayerRegistry.instance().mapLayer(u'mylayer_id')

cache = QgsVectorLayerCache(layer, 10000)
model = QgsAttributeTableModel(cache)
print model.setHeaderData(0, Qt.Horizontal, u"pk")
print model.setHeaderData(1, Qt.Horizontal, u"name")
model.loadLayer()
table = QTableView()
table.setModel(model)
table.show()

The problem is the fact that header in the view are not updated with model.setHeaderData: it returns false (meaning that the header was not updated) whereas with other models (from Qt), I don't encounter this issue.
You may wonder why I make a reference to other Qt models? It's because the setHeaderData method is inherited from Qt QAbstractItemModel.

Is is due to the QGIS cache? How can I solve this issue? Any hints is welcome


texts:
objective :

Save time in setting geoserver configuration specially in case like updating configuration from development environment to production or quickly create a copy for experimentation

What i know :

I can deploy geoserver by copying the war file in the tomcat webapps folder then tomcat automatically creates geoserver folder

What i want to know :

what if i copied the geoserver folder from say tomcat on server 1 to tomcat on server 2 .
Will it work ?
Any precautions ?

Note : I am for now NOT setting geoserver data dir in path


texts: Setting the Python interpreter in PyCharm 2020.2.3, using QGIS 3.10.11. Error "No module qgis exists"

texts:
I'm trying to find out where I'm going wrong in trying to set up PyCharm to make standalone QGIS Python scripts.
The very first step in the tutorial is "from qgis.core import *"
I've tried each Python installed in the QGIS folder, they all have no module named qgis.

The error is the same


texts: Having a problem trying to get WMSGetFeatureInfo to work in openlayers - running off of Geoserver (http://gisnexus.com/cmha/cmha_site_test.html) you click on any feature, it shows as if the layer is being read from opengeospatial.net - but nowhere in the code does it call out to it!


texts: How to calculate walking distance between multiple addresses? - QGIS

texts: I am working with historical data and I would like to calculate the walking distance between multiple addresses. For instance, in my excel dataset I have someone who lived on 1 Laurierstraat in Amsterdam and went to the pub on 10 Rozengracht, Amsterdam. I have hundred of cases like that, combining the source and the destination. Is there an easy way (for someone with no programming knowledge) to import these data and calculate the walking distance for each set?


texts:
I'm having trouble getting a gray basemap of Philadelphia in a vector scalable format using qgis. I've tried using the openstreets basemap plugin, and then taking that to print composer, but that doesn't seem to work as it aught to. I've also uploaded shapefiles of rivers, roads, and railways to try to make my own basemap, which looks nice, but doesn't really work when exported to svg either.

Are there collections of svg basemaps available?

Is there a better way to solve this problem than the ways I attempted?


texts:
I have installed Mapnik and successfully tested it. Now I am facing the problem of creating my own customized tiles. I have googled it but everywhere it was telling about OSM tile generation.

My vision is simple: I have a huge map data of my target area. I want to create my own custom map in Mapnik/Tilemill. Then create tiles and upload them into server. A step by step instruction would be helpful to me since I am a newbie in Mapnik.

The main problem: How can I create tiles with my own data?


texts: Map netncf-file with different years and save the different layers in shapefiles

texts:
Download from this website the following netncf-file: cru_ts3.23.1991.2000.pet.dat.nc.gz. This is a monthly raster with data from all over the world from january 1991 to december 2000

I am using ArcGis (multidimensional tool) to plot the map. However, I don't manage to plot the maps per month. I only receive one map, without time indication. I followed the instructions here to build a time slider but the time slider doesn't work in the sense that it is showing time periods (1900 for instance) which is not in the data file I downloaded, and it doesn't slide.

A follow up question on this is that, once the time slider works, is it possible to save the different time layers seperately in different shapefiles such as in this question. I would like to use the separate shapefiles afterwards in R.


texts:
I am working on a new project where I need to highlight census tracts by boundaries. Google Maps can do this as long as I give it GeoJSON (data layers).

I did some research and it seemed like Tiger is a good solution. I was able to install and setup PostGIs and a Tiger Geocoder database (125 gigabytes!) https://postgis.net/docs/Extras.html

After going through this incredibly long process, I could finally geocode an address using the built-in function geocode. I could even return back GeoJSON from an address:

SELECT ST_AsGeoJSON(geomout)
FROM geocode('8615 Chestnut Oak Rd, Baltimore, MD 21234', 1)

and it spits out GeoJSON:

{"type":"Point","coordinates":[-76.5617440647409,39.3951950048467]}

Finally, using the get_tract function, I can get back a census tract like this:

SELECT get_tract(geomout, 'tract_id') As tract_id
FROM geocode('8615 Chestnut Oak Rd, Baltimore, MD 21234', 1)
As g;`

Which returns a tract.

The Question:

Now that I've shown what I've been able to do, what I can't do is get the GeoJSON for this tract. Would anyone be able to nudge me in the right direction?

I am wondering if maybe there is a function to return back a geometry object from a census tract (which I can then pass into the ST_AsGeoJSON function, or some other solution).


texts: grass gis i.atcorr input should be reflectance or radiance

texts:
My question is what should be the input image in i.atcorr for grass 6.4.4 version?

I tried giving both reflectance (DN) and radiance as input but values are coming properly.  I read some where, it internally converts to radiance and uses.  Is it like that? If I give radiance as input, the values are not coming properly in shadow region effect is more.  So I used i.landsat.toar and given as input to i.atcorr without marking the flag input converted to reflectance(default is radiance).  But the output values are not correct for flat terrain and hilly terrain (both tested)

Also one more I need to know is if I would compare the output of i.atcorr with surface reflectance product (from landsat website) will it be proper or not?  I have used srtm dem and atmospheric parameters from modis website. So all the parameters are given, but is it correct way to verify our final output or any other way is there?  I will be highly thankful to you if any solution is provided.


texts:
My map is not exporting the right coordinates.
As you can see it makes a local sistem.
It should be around 46 N , 20 E.

Do you have any idea what am I doing wrong?


texts: Which SRS to use for global distance queries?

texts:
I have a database with a lot of points in WGS84.
Now I am building an cache that performs NN and points in range queries using a KDtree. The point [sic] is that the search radius will be provided in meters and that lat/lon is not a nice SRS for these geometric queries.

I am looking for a geometric SRS that is applicable to the whole world and that preserves distances. I don't care about errors of a couple of 10ths of meters.


texts: What's wrong with my Planet API geometry filter?


texts: Updating job of workflow manager using ArcGIS API for JavaScript?

texts:
I'm trying update the name of a job:

updateJob: function (jobId, workArea) {
    let loParameters = new JobUpdateParameters();
    loParameters.jobId = jobId;
    loParameters.name = "Updated";

    let loJobTask = new WMJobTask(poAppVars.urlWFMServer);
    loJobTask.updateJob(loParameters, goHelper.getValueByKeySO("UserName"),
      goHitch(this, function (response) {
        alert('¡¡¡Updated!!!');
      }),
      function (err) {
        console.log('err');
        console.log(err);
      });
  }

But I get an error saying:

  Unable to complete operation

Can anyone tell me what's the problem?


texts: Seeking interactive maps and associated technologies?

texts:
The display of maps and data within an interactive framework is becoming more prevalent online. I'm thinking beyond merely the ability to pan/zoom/control layers, but more along the lines of displaying spatial and non-spatial data, together, in a unique format (e.g. animated charts).

For example, this website goes beyond mapping and presents spatial data as a learning and exploration tool with graphical data summaries, text summaries, animations, etc.

What are some good examples of well-designed interactive maps using modern technologies, and what are the advantages/disadvantages of these tools? What technologies should be avoided?


texts: Triangulation: Estimate a location using lat/long coordinates and bearings

texts:
I would like to get an estimated location using lat/long coordinates and bearings at each of those coordinates using R.
In this case I would be approximating the location of a turkey nest that I have marked using radio telemetry by taking 5-8 bearings around the nesting hen so that the nest can be found after the hen has completed the incubation period and left.
Here is an example of what this might look like when plotted in google earth (northern most point not visible):

Full list of coordinates and bearings for the above example:
   Coord.(UTMs)......Bearing.(Aizmuth).

1....473789, 4485636.......26
2....473755, 4485659.......50
3....473728, 4485717.......80
4....473730, 4485788.......134
5....473848, 4485866.......196
6....473926, 4485796.......236
7....473925, 4485702.......291
8....473850, 4485660.......326


texts: Text or tooltip for QGIS drag and drop form designer

texts:
Is there a way, to add tool tips or a simple text next to a field widget in QGIS forms drag and drop designer?

I would like to provide more information about some fields (restrictions, expected value etc.) to the user to help him, filling in the form.


texts: Openlayers 3 setvisible not working until map clicked

texts:
In this application layers are constructed dynamically from values stored in a JSON object. There is a layer control function that is called from buttons which toggle the visibility of the groups and the individual layers within the groups.

The problem - When an individual layer inside a group is toggled, the change isn't applied to the map until the user clicks on the map view-port. So when the user switches a layer on/off nothing happens.

It's fine if layer constructors are hard coded. It's only when building the layers from JSON that this happens. Also, if map.updateSize(); is called just after the visibility is set then the map visibility changes as required but this isn't how setvisible is meant to be used.

Is there a way to fix this issue with setvisible other than calling map.updateSize() to force it to apply the change?

I have created a fiddle demonstrating the issue here:
https://jsfiddle.net/Rob_H/wvdmex8m/5

HTML

           <input type="button" onclick="toggleLayer(basemaps)" value="Group"/>

                <input type="button" onclick="toggleLayer(osm_baselayer)" value="Layer">

Javascript

var response = {
  lyrs: [{
    subLayer: "osm_baselayer",
    name: "Open Street Map",
    visible: 1
  }],

  lyrGroups: [{
    groupName: "basemaps",
    visible: 1
  }]
};

var mapLayers = [];

function makelayerGroup() {
  var lyrGroups = response.lyrGroups;
  for (i in lyrGroups) {
    var layerGroup = lyrGroups[i].groupName;
    var vis = lyrGroups[i].visible;

    this[layerGroup] = new ol.layer.Group({
      visible: vis,
      name: layerGroup,
      layers: []
    });

    mapLayers.push(this[layerGroup]);

  };
  return makelayers();
};

function makelayers() {
  var lyrs = response.lyrs;
  for (row in lyrs) {
    var group_name = lyrs[row].groupName
    var sub_layer = lyrs[row].subLayer;
    var vis = lyrs[row].visible;

    //layer constructor
    this[sub_layer] = new ol.layer.Tile({
      visible: vis,
      name: lyrs[row].name,
      source: new ol.source.OSM()
    });

    var a = mapLayers[i].getLayers().getArray();
    a.push(this[sub_layer]);
  }

  return mapLayers;
};

//intialise the map
var map = new ol.Map({
  target: 'map',
  layers: makelayerGroup(),
  view: new ol.View({
    center: ol.proj.transform([-0.7, 52.1], 'EPSG:4326', 'EPSG:3857'),
    zoom: 11
  })
});

//function to toggle layers
function toggleLayer(layerName) {
  if (layerName.getVisible() == 1) {
    layerName.setVisible(0);
    console.log(layerName.get('name'));
  } else {
    layerName.setVisible(1);
    console.log(layerName.get('name'));
  }
};


texts: Is there a tool to shrink/grow raster data edge by n pixels

texts:
Looking for an open source tool to grow/shrink raster edge by n pixels. I noticed Grass has r.grow but hoping there might be a way to do it with GDAL, PostGIS or python library.

More information:
I want to shrink it around the edge where I have data not the entire rectangular region.
Consider this raster:


texts: leaflet Uncaught TypeError: Cannot read property 'addLayer' of undefined

texts:
I'm getting the following error:

 Layer.js:48 Uncaught TypeError: Cannot read property 'addLayer' of undefined
        at NewClass.addTo (Layer.js:48)
        at ResGarbGeoJ63.html:84

From this code:

var map;

        function init() {
          map = new L.map('mapid');
          map.setView([37.39,-122.086],15);

             // Add the tiled layer
            var tiles = L.tileLayer('http://{s}.tile.osm.org/{z}/{x}/{y}.png', {attribution: "Data copyright OpenStreetMap contributors"});
            tiles.addTo(map);

            var mondayLayer;
            var tuesdayLayer;

            var selection;
            var selectedLayer;

            // add layer

              map.on('moveend', function(e) {
                //recalculate visible points
                recalculateIncludedPoints();
              })

          }

          // Set up symbols for stops
            // Purple symbol for unselected stops
            function mondayStyle(feature) {
              return {
                fillColor:'#ff00ff',
                fillOpacity: 1,
                color: '#B04173',
              };
            }

            // Yellow symbol for unselected stops
            function mondaySelectedStyle(feature) {
              return {
                fillColor: "00FFFB",
                color: '#0000FF',
                fillOpacity: 1
              };
            }

            // handle click events on the garden features
            function garbageOnEachFeature(feature, layer){
              layer.on({
                click: function(e) {
                    if (e.target.feature) {
                      //this scope will run whenever the clicked target is a valid marker point
                      resetStyles();
                    }

                    //set this style
                    e.target.setStyle(mondaySelectedStyle());
                    selection = e.target;
                    selectedLayer = mondayLayer;

                    // Insert some HTML with the feature name
                    buildSummaryLabel(feature);

                    L.DomEvent.stopPropagation(e); // stop click from being further propagated
                  }
                });
              }

              var mondayLayer = new L.geoJSON(mondayData, {
              style: mondayStyle,
              onEachFeature: garbageOnEachFeature
            });

            mondayLayer.addTo(map);
            tuesdayLayer.addTo(map);

I move the var from the var mondayLayer = new L.geoJSO... but to no avail


texts: Projection parameters different for British National Grid between QGIS and ArcMap?

texts:
I got shape files in WGS84 and need to project it to British_National_Grid (27700)
I used QGIS to save them as a new shp file with ESPG:27700

Throwing the shape files into ArcMap with the Data frame Coordinate System set to British_National_Grid
WKID: 27700 Authority: EPSG

When I try to edit the shape files I get a warning that the data is in a different coordinate system than the dataframe.

I was expecting both to be in 27700? Do they use different transformation paramenters to do the transformation which ArcMap catch or why do I get this warning?


texts: Copying/Inserting features AND user-defined Attributes as input to existing feature class in SDE using ArcPy?

texts:
import os, arcpy, zipfile, glob

# INPUT ZIP FILE
infile = arcpy.GetParameterAsText(0)
path = "D:\Database_Connections"

NAME = arcpy.GetParameterAsText(1)
CODE = arcpy.GetParameterAsText(2)
SOURCE = arcpy.GetParameterAsText(3)

if not os.path.exists(path):
    os.makedirs(path)

# DATABASE CONNECTION PARAMETERS
arcpy.CreateDatabaseConnection_management ("D:\Database_Connections", "sample.sde", "SQL_Server", "@@@.@@.@@@.@@@", "DATABASE_AUTH", "dbname", "dbpassword",  "SAVE_USERNAME", "sample","#","TRANSACTIONAL","#","3/20/2018 6:45:41 PM")
arcpy.env.workspace = "D:\Database_Connections\sample.sde\sample.DBO.samplefeatureclass"

# EXTRACT ZIP FILE TO THE Database_Connections FOLDER
zip = zipfile.ZipFile(infile, 'r').extractall(r'D:\Database_Connections')

# CODE TO SELECT POLYGON SHAPEFILE OR POLYGON FEATURE CLASS
path = "D:/Database_Connections"
fc = []
for dirpath, dirnames, filenames in arcpy.da.Walk(path,
                                              datatype="FeatureClass" or "ShapeFile",
                                              type="Polygon"):
    for filename in filenames:
        fc.append(os.path.join(dirpath, filename))

    # Open a search cursor on a feature class / table to read from
    scur = arcpy.SearchCursor(read_fc)
    # Open an insert cursor on the fc / table to write to
    icur = arcpy.InsertCursor(write_fc)

CODE TO INSERT NAME, CODE AND SOURCE VALUES TO THE samplefeatureclass

I am trying to add polygon feature and the attribute (taken as an input from the user) to the existing feature class in the SDE (MSSQL database) using arcpy. I am able to add an empty row to the feature class table in SDE, can someone guide/help me to pass the user-defined attribute information (this is taken as an input from the user) and insert in the feature class table in SDE.

INPUT : Polygon shapefile & attributes defined by user as parameter.

OUTPUT : Copying/Inserting polygon feature alongwith user-defined attributes to
         the existing feature class in SDE.

These can be seen in the attached screenshot.


texts:
I have QGIS 3.4.14-Madeira installed. I'm trying to access OpenStreetMap using QGIS.

The instructions at Using OSM Data in QGIS go to Vector → OpenStreetMap but this is not in the drop list that appears under Vector.

Instead I see:

Geoprocessing Tools
Geometry Tools
Analysis Tools
Data Management Tools
Research Tools

How do I access OSM data from QGIS?


texts: Find a proper path between two nodes in a grid

texts:
This question is about finding a proper path (route) between two nodes (P,Q) in a grid (i.e., graph, network, ...). In our concept, we define a path proper (nice) if it is the most similar* to the straight line connecting the two nodes, as shown in the following figure. So, from 10 possible short paths (all are 5 segments in length,) between green and red nodes, number 5 is chosen the best one.

most similar: the similarity between two polylines (here a polyline and a line) can be defined in various ways (see this) for example. So vertices of a polyline can be evaluated with regards to the line.

Suppose we have locations as grid nodes, (1) we are looking for implementations of the method explained above in any programming languages. Source code, coding hints, ideas etc are more than welcome. We prefer not to use existing graph libraries but are interested in digesting the solution with examples.

The method needs to be flexible to (2) adding some conditions such as "not visiting a node" as shown in the following figure (yellow node). That is, the best choice from above needs to be slightly adapted to satisfy the required condition.


texts:
I recently installed GRASS 6.4 in a Unix environment. The application opens fine, but I would like to also be able to run the GUI from the command line. The tutorials says that the
grass64 command should bring up the GUI from the command line, but the command is not recognized, so I assume I need to add the path variable to my .bash_profile. I have looked everywhere and cannot figure out in which directory grass64 is located. Does anyone know what path I need to add to make it work?

I tried adding the /Applications/GRASS-6.4.app/Contents/MacOS/ path to my .bash_profile, since this seems to be the path that the app uses, but it did not work.


texts: I am struggling to come up with a proper layer style for a rail network extracted from OSM in QGIS. I just want active rail so I choose the "rail" tagged features (rule-based rendering). But since every track is a line and escpecially within cities, there a sometime 5 or more tracks close to each other so the symbology looks pretty ugly like in the screenshot. What could I do besides, deleting the features I do not want to be shown (which would mean a lot of work)?


texts: Use lapply and writeOGR to write multiple shapefiles within a list

texts:
I have a nested list with multiple shapefiles and I want to use lapply with writeOGR ('rgdal') in order to export all the shapefiles. Sorry for not giving a reproducible example.

# My list
shp.list
$e_ang
$e_ang$e_ang_g
class       : SpatialPolygonsDataFrame
features    : 1291
extent      : -1265736, 869263.7, -3051349, -1022949  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=eqdc +lat_0=0 +lon_0=25 +lat_1=20 +lat_2=-23 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
variables   : 1
names       : e_ang_g
min values  :       1
max values  :       1

$e_ang$e_ang_l
class       : SpatialPolygonsDataFrame
features    : 1291
extent      : -1265736, 869263.7, -3051349, -1022949  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=eqdc +lat_0=0 +lon_0=25 +lat_1=20 +lat_2=-23 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
variables   : 1
names       : Value
min values  :     1
max values  :     1

Here is the code I've used

# Using lapply to write shapefiles within list
lapply(unlist(shp.list), function(x) writeOGR(object = x,
                                         dsn = paste(fruit.bats, "output", sep = ''), # folder to store shapefiles
                                         layer = names(x), # shapefile names as in shp.list
                                         driver = 'ESRI Shapefile',
                                         layer_options = 'RESIZE=YES',
                                         overwrite_layer = T))
# this returns

$e_ang.e_ang_g
NULL

$e_ang.e_ang_l
NULL

I believe the problem is in the argument 'layer' in writeOGR. Probably I can create a vector to pass the names to the function, but I was wondering if there is a more direct way to do it. I would like the shapefiles to have the name of the nested elements, 'e_ang_g' and 'e_ang_l' in the example.

I was able to do the same with a loop, but still seek a solution using lapply.

Here is the code for the loop.

n <- substring(names(unlist(shp.list)), 7,13) # vector with shapefile names

for (j in 1:length(n)){
  writeOGR(unlist(shp.list)[[j]],
           dsn = paste(fruit.bats, "output", sep = '/'),
           layer = n[j],
           driver = 'ESRI Shapefile',
           layer_options = 'RESIZE=YES',
           overwrite_layer = T) }


texts: Change Ordering of features within one layer in leaflet

texts:
My Question is exactly the same as mentioned in this Question

The main difference is that I can't use bringToFront() because I have points inside the highlighted polygon that I want to show a tooltip on hover.

So I want the polygon feature just to be the uppermost, in his feature layer.
Something like bringToFront() for within the layer only.


texts: How to check if a node is in route computed with pgRouting?

texts:
I'm new to the GIS system and I'm following the pgRouting workshop.

I've imported the locations of nodes (points, highway=street_lamp) using osm2pgsql.

now I'm calculating a few paths between two points using pgr_ksp (K-Shortest Path) to compute K shortest paths.

How can I count how many street lamps there are for each path?

There's a function that given a path returns all the nodes in that path?


texts: I used projected (Hungarian EOV) orthophoto for segmentation and classification in eCognition (I defined EOV in ecognition as well), the result was exported into shapefile. I opened orthophoto and shapefile in QGIS and ArcGIS, but exported shapefile doesn't overlap with orthophoto even they have the same coordinate system. Can anyone help me to deal with this issue.


texts: Creating simplified QGIS version with custom interface?

texts:
At my institute we are revising and updating a methodology for participatory mapping of land degradation, that involves creating a simple land use map, subdividing each land use unit in a number of sub-categories (e.g north  vs. south exposed areas) and assigning quantitative and qualitative values to a certain number of pre-defined attributes (e.g cause of degradation, degradation level, counter measures implemented, ...).

We are exploring the idea of providing users with a custom made version of QGIS with a very simplified interface, and custom buttons to guide the users to the different steps of the evaluation.

Is this technically possible and allowed by the license terms?


texts:
I just (December 12, 2018) downloaded the training manual and data. I'm trying to get started but I cant find the first part to download the OpenStreetMap part.

"In the Vector menu dropdown, select OpenStreetMap -> Download Data. You can then manually enter the co-ordinates of the region you wish to use, or you can use an existing layer to set the co-ordinates."

I see the Vector menu on the top ribbon, but the dropdown does not have the right option. I cruised around but didn't see a way to add it.


texts: Connecting circles with lines cover all circles using PostGIS?

texts:
I'm really new to postgis.

I have trouble to find out the coding for connecting circles.
You can figure out what I'm trying to do from the image I uploaded.
I want to connect the circles which have the same numberings by outerline.

My first trial was using st_convexhull, st_concavehull but both of them does not make the covering polygons as I wanted.

like the following:

create table public.cover_line as
SELECT ST_Boundary(ST_ConcaveHull(geom,0.99)) as convex
FROM circles;

Is there someone can help me?


texts: WMS/WFS layers from QGIS 2.12 projects fail when loading the project in QGIS 2.18

texts:
I've been usign QGIS Server 2.12 and QGIS Desktop 2.12 for several months without problems.

Now I installed QGIS 2.18.3 to use some new features I consider them interesting but when I open my QGIS project the WFS layers (served by QGIS server 2.12) fail to load.

Is there any WFS-T protocol update in >2.12 versions?

Is it because my WFS urls have the VERSION=1.0.0 parameter and maybe QGIS Desktop 2.18 does not handle 1.0.0 requests?

EDIT 18.01.2017:

If I add the WFS layer in QGIS 2.18 it adds ok. The problem comes when trying to open a QGIS project saved in 2.12 with QGIS 2.18.

I added a layer to canvas and saved the project in QGIS 2.12, the layers' datasource is as follows:

<datasource>http://eudala2.getxo.eus/qgis/qgis_mapserv.fcgi.exe?map=D:\OSGeo4W\apps\qgis\bin\test\vialidad_borrar.qgs&amp;SERVICE=WFS&amp;VERSION=1.0.0&amp;REQUEST=GetFeature&amp;TYPENAME=getxo_turismo_pois_old_borrar&amp;SRSNAME=EPSG:4230</datasource>

If I do that in QGIS 2.18 with the exact same WFS layer I can see that the data in datasource is splitted into different params:

<datasource> restrictToRequestBBOX='1' srsname='EPSG:4230' typename='getxo_turismo_pois_old_borrar' url='http://eudala2.getxo.eus/qgis/qgis_mapserv.fcgi.exe?map=D:\\OSGeo4W\\apps\\qgis\\bin\\test\\vialidad_borrar.qgs' table="" sql=</datasource>

What could I do to open my QGIS 2.12 projects with hundreds of layers (without having to create the projects in QGIS 2.18 again) ?


texts:
(I've been researching this for hours and have yet to find a straight answer, so at this point I don't care if this gets marked duplicate.)

I'm a complete noob to QGIS. I have a floorpan SVG that I cleaned up in Inkscape and exported as DXF. The coordinate system of the DXF is in point space and I have the (lat,lon) transformation of the 4 corners. When I import the DXF as a vector layer, I don't see any way to set a custom transformation to WGS84. Ultimately I want to be able to export the geometry as Shapefile or GeoJSON.

Most of the other responses I've researched are either out of data (plugins that no longer exist, tutorials wherein the instructions are out of date, etc.) This shouldn't be this difficult.


texts: Custom Attribute Display image errors in ArcGIS online popups

texts:
I am trying to add images in the Custom Attribute Display of arcgis.com map editor as opposed to the conventional Pop-up Media area because I want to organize them in a custom way.

However no matter what I try, it always says 'The Value entered is not valid'.
I want the photo in the popup so I don't want to use the link option.

I have tried having "http://" in my field and not. When using related table fields, the add image box would append it, but doesn't seem to when using non-related table fields.


texts: I have an existing attribute table that I'd like to add values to the fields within it while preserving the geometry, how can I do this?

texts:
I've tried to use fiona and shapely and ogr, but the tutorials seem to be insufficient to describe how to do such a seemingly simple task.  In arcpy this could be done simply with an update cursor.  I'm trying to use open source tools here but I'm having some difficulty with editing the attribute table.  I'm not interested in creating a new shapefile.

Any suggestions?

(Edit) For example how can I implement the following code using open source tools?

cursor = arcpy.da.UpdateCursor(shpfilepath, ["OID@", "Elevation"])
for row in cursor:
    row[1] = 15
    cursor.updateRow(row)
del cursor

This would change the values of the "Elevation" attribute of all features of shpfilepath to be 15.

(Edit 2)
I'm trying to run this as an external script, not as part of QGIS.
(Edit 3)
I ended up doing it like this:

def addfields(shapefilename, fieldNames):
    #adds fields and values to shapefile
    #using fieldName, a dictionary that has the structure:
    #["ColumnName":"Value"]

    #filehandling to open layer
    shape = shapefilename
    driver = ogr.GetDriverByName('ESRI Shapefile')
    dataSource = driver.Open(shape, 1) #1 is read/write
    layer = dataSource.GetLayer()

    #add fieldnames using ogr
    for field in fieldNames:
        fldDef1 = ogr.FieldDefn(field, ogr.OFTString)
        fldDef1.SetWidth(16) #16 char string width
        layer.CreateField(fldDef1)

    #iterate through features using fid and add values with ogr
    for fieldName in fieldNames:
        for fid in range(layer.GetFeatureCount()):
            feature = layer.GetFeature(fid)
            feature.SetField(fieldName, fieldNames[fieldName])
            layer.SetFeature(feature)


texts: QGIS export png with accurate hex color values

texts:
I have a 32 bit TIFF raster layer in QGIS with one band of integer values between 49 and 72 with some no data values. I'd like to end up with an 8 bit PNG image with matching pixel values. I don't care what no data ends up as as long as it's consistent. I've tried using GDAL translate, something like:

gdal_translate -ot Byte -scale 0 255 -of PNG "in.tiff" "out.png"

but the values I get back in the png range from 38 to 59 instead. Same range, but offset -13. Any ideas why, or how to do this properly?

Edit:
If I specify -scale 49 72 62 85, accounting for the offset I'm seeing, the values are correct. So that works, but it would be great to know why it's shifting -13.


texts: What might cause invalid syntax from ArcPy with Raster Calculator and Multiple Rasters?

texts:
I am having trouble for writing script in Python (ArcGIS 10.1). The task is generally divide each raster by 100. I have many of these raster files...

So what I've done is the following:

from arcpy.sa import *
  arcpy.workspace="C:/Users/Vaio/.../Monthly_Original.gdb"
  rasterlist=arcpy.ListRasters(p)
  for k in range():
     OutRaster=(Raster(rasterlist[k])/100
     OutRaster.Name='p_2000_*.*'
     OutTaster.Save(OutRaster.Name)

I have an Error in line 6. I am very new to any sort of programming.


texts: Looking for smallish raster of global elevation data

texts: I'm trying to find elevation data for the whole globe for a simple mapping application and I'm not having any luck.  I really don't want anything fancy, just something like a 4000x2000 raster of elevation values across the globe would be fine.  I plan to colormap and interpolate it myself.  Does anyone know where I might find something like this?


texts: Determine Los Angeles neighborhoods based on coordinates

texts:
I'm new to the GIS world. I've read several sites and previous posts but couldn't find a satisfactory answer.

I have a data set of ~400k coordinates from Los Angeles and would now like to determine their respective neighborhood.

Lat            Lon
-118.4382362   34.2530534
(...)

What would be the best way to determine the coordinates neighborhood?

EDIT: I started working with ArcMap and plotted some test coordinates on a shapefile.

How would I now match the neighborhood names from the shapefile to my coordinate data? When I use the join function, it only counts how many coordinates occur in a specific neighborhood. However, I'd like to know the opposite i.e. the respective neighborhood name per coordinate pair.


texts: Mac terminal and r show different gdal versions

texts:
I am trying to read a .hdf file into R using:
gdalUtils::get_subdatasets()

but I am left with this error:
ERROR 4: `/var/folders/24/9_h4zhq564g3zpvqp6jmbm0c0000gn/T//RtmpibWWKQ/MODIStsp/HDFs/MCD12Q1.A2014001.h18v04.006.2018146020544.hdf' not recognized as a supported file format.
gdalinfo failed - unable to open '/var/folders/24/9_h4zhq564g3zpvqp6jmbm0c0000gn/T//RtmpibWWKQ/MODIStsp/HDFs/MCD12Q1.A2014001.h18v04.006.2018146020544.hdf'.
Error in h(simpleError(msg, call)) :
  error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': subscript out of bounds
In addition: Warning message:
In system(cmd, intern = TRUE) :
  running command '"/usr/local/Cellar/gdal/3.2.0/bin/gdalinfo" "/var/folders/24/9_h4zhq564g3zpvqp6jmbm0c0000gn/T//RtmpibWWKQ/MODIStsp/HDFs/MCD12Q1.A2014001.h18v04.006.2018146020544.hdf"' had status 1

I have had issues with gdal before, and was forced to download and use an older version of gdal to get the above function to work properly. In mac terminal, running
gdalinfo --version

returns
GDAL 2.1.3, released 2017/20/01

I assume there is a build problem here. When I run the following command in the R Console:
getGDALVersionInfo()

I am returned:
"GDAL 3.2.0, released 2020/10/26"

Is there a method to "unlink" gdal from R and then relink it to the older version noted in my Mac terminal?
There is a similar problem posted here, but the (unaccepted) answer is ambigious: "You need to uninstall gdalinfo on your mac environment and use R gdalinfo and dont forget to restart your R after uninstall it"


texts: Zoom to spatial coordinate in OpenLayers load

texts:
I am working on a project that uses OpenLayers to show shapefiles from the GeoServer and then show them on Google Maps. The only problem is when the page loads, I want to go to the spatial coordinate and zoom to the shapefile that is loaded on map. Is there any way to do this? I want to go to the coordinates around 32 and 50 in zone 39N. Can you help me about this problem?

I tried this code:

var dataExtent = states.getDataExtent();
map.zoomToExtent(dataExtent);

But it does not work, where is my mistake?

This is the result:

h++p://173.244.165.194/

It loads many maps of world.  I have to zoom to go to my location :(
Here is exact place that I want to zoom when page loads:

http://173.244.165.194/?zoom=12&lat=31.62211&lon=51.48708&layers=B00T


texts:
I have read pixel values of an raster image using GDAL libraries in visual studio 2010(vc++).

Now I have to crop the image (subset) according to the grid given in shape file.

I have got a code for cropping image on internet as given below, but its having some error in it.
Please review it and help .

GDALDataset* crop(GDALDataset* poDS, int xoff, int yoff, int xsize, int ysize)
{
        VRTDataset *poVDS = (VRTDataset*)VRTCreate(xsize, ysize);

        // Copy dataset info
        const char* pszProjection = poDS->GetProjectionRef();
        if (pszProjection != NULL && strlen(pszProjection) > 0)
                poVDS->SetProjection(pszProjection);

        double adfGeoTransform[6];
        if (poDS->GetGeoTransform(adfGeoTransform) == CE_None)
        {
                 // Adapt the geotransform matrix to the subarea
                    adfGeoTransform[0] += xoff * adfGeoTransform[1]
                        + yoff * adfGeoTransform[2];
                adfGeoTransform[3] += xoff * adfGeoTransform[4]
                        + yoff * adfGeoTransform[5];

                poVDS->SetGeoTransform(adfGeoTransform);
        }

        poVDS->SetMetadata(poDS->GetMetadata());

        // Here I also copy metadata from my own domain
        char **papszMD;
         papszMD = poDS->GetMetadata(MD_DOMAIN_MSAT);
               if (papszMD != NULL)
                 poVDS->SetMetadata(papszMD, MD_DOMAIN_MSAT);

         for (int i = 0; i < poDS->GetRasterCount(); ++i)
         {
                 GDALRasterBand* poSrcBand = poDS->GetRasterBand(i + 1);
                 GDALDataType eBandType = poSrcBand->GetRasterDataType();
                 poVDS->AddBand(eBandType, NULL);
                 VRTSourcedRasterBand* poVRTBand = (VRTSourcedRasterBand*)poVDS- >GetRasterBand(i + 1);
                 poVRTBand->AddSimpleSource(poSrcBand,
                                 xoff, yoff,
                                 xsize, ysize,
                                 0, 0, xsize, ysize);
                 poVRTBand->CopyCommonInfoFrom(poSrcBand);

                 // Again, I copy my own metadata
                 papszMD = poSrcBand->GetMetadata(MD_DOMAIN_MSAT);
                 if (papszMD != NULL)
                         poVRTBand->SetMetadata(papszMD, MD_DOMAIN_MSAT);
         }

         return poVDS;
 }


texts: Leaflet: Are Bing Maps tiles exempt from restriction to use only with Bing Maps?

texts:
Leaflet has a plugin that can fetch Bing Maps tiles. However, I see in the Bing Maps terms of use the following prohibition:

    8.2 d) Use Content, including geocodes, other than via an authenticated call to the Services and/or in conjunction with a Bing Map.

Would this not make any use of the leaflet plugin a violation of the terms of use? Are tiles specifically different from other Bing data? Or is there a specific exemption that has been offered to Leaflet?

Further, if other services are used, e.g. geocoding with OpenStreetMap data, is that a breach of the above?


texts: Using custom Python function with QGIS raster calculator

texts:
I have a function in a python class that takes two values and returns a third value. I would like to apply the function to two overlapping rasters on a per pixel level and come up with a third raster. It doesn't appear that the QGIS or raster calculator can accept a custom function. I'm sure this can be done with GDAL, but don't have any experience with it. GDAL  is certainly an option if an easier system doesn't exist.
For a little context I'm using an NDVI raster and a vegetation height raster to calculate hydraulic surface roughness based on a table I've already created. Worse case situation will be to write a quick script to convert the table to a raster calculator equation. As the table has 49 "bins" the equation will be quite long and it seems like there should be a more elegant way to do this.

Simplified code for returning n value based on NDVI and vegetation height, a lot of error checking has been removed for clarity:
class RoughnessTable(object):
    """ Extracts predicted mean and std Manning's n values for a XS based on lidar and NDVI v0.01b """
    def __init__(self, table):
        """ table is the lidar/ndvi/n-value dataframe"""
        self.table = table

    def mean(self, l, n):
        """
        Return mean value from roughness table for lidar height 'l' and NDVI 'n'
        :param l: mean lidar height of xs segment
        :param n: mean ndvi of xs  segment
        :returns: mean predicted surface roughness (float)
        """
        rt_row = self.table.query('l_bin_min <= @l < l_bin_max & n_bin_min <= @n < n_bin_max')
        return float(rt_row.mean_)

The table passed to RoughnessTable.__init__() is a pandas dataframe that represents the NDVI/vegetation height bins. A sample is shown below.

I've had success using the raster_calculator from the pygeoprocessing library, which handles opening the raster, and splitting it into 256 x 256 pixel ndarray chunks. I'm still having to iterate over the chunk from raster_calculator, and it seems rather hacky, but it works. It is rather slow. I haven't profiled it yet, but I'm guessing the self.table.query(...) line is causing a large portion of the problem.


texts: Cannot publish Composite Locator to ArcGIS Server 10.2.2

texts:
I am unable to publish a composite locator from a 10.0 SP5 ArcSDE Geodatabase (Oracle) to ArcGIS Server 10.2.2. To narrow down the problem I created a fresh composite locator with only one component locator to keep it as simple as possible.  I am able to successfully publish the component locator but not the composite.

When publishing, the packaging is successful and then at the very end of the process I get the message "Failed to publish service".  In the ArcGIS Server logs I get the errors below:

Failed to create the service.: ERROR: code:404, Could not find resource or operation
'Composite_Address_Streets.GeocodeServer' on the system., No resource could be found
at that address.

Error executing tool.: ERROR 001487: Failed to update the published service with the
server-side data location. Please see the server's log for more details. ERROR
001369: Failed to create the service. Failed to execute (Publish Service Definition).

I am however able to create the .sd file.  If I try to upload it to server using the Upload Service Definition GP tool I get the same error as above.

A few more notes:

I am publishing the service through the "owner" connection to SDE so the privileges should not be an issue. Also, the component locator is referenced in the composite via the owner connection as well.

I have registered the SDE data source correctly and do not get any suspicious messages during the analysis step.  The only message is that the composite locator will be copied to the server which is normal since they are always copied to the server.

I am able to publish the same composite locator with a copy of the data from a 10.2 file geodatabase so it looks like it's just an SDE/ArcGIS 10.2.2 issue.


texts:
I am new to Arcade and I am trying to symbolize a feature based on a domain value in multiple fields.

The feature class is named Zones and there are 3 fields Phase1, Phase2, Phase3 with the same Domain assigned to each field. The Domain values are [Active, Complete, Future]

I want to create an expression when Phase1, Phase2 or Phase 3 = Active

I am using a hosted feature class on Enterprise 10.6.1


texts: I was wondering if someone can help me in the right direction to making a geospatial report. I used a drone to take aerial photos and used MME to stitch the photos together. I then uploaded the Geotiffs to QGIS with the Google Earth background, created a new print layout of the location map and exported the map to PDF. However, when I use the Geospatial Location Tool, the coordinates cannot be determined. I even tried to create a PDF from a .TIFF file and still no luck.


texts: g.renderer is null. And sometimes is not

texts:
I have an TypeError: g.renderer is null while I'm loading my map application. The worst is that problem sometimes happens, sometimes not.

require([
    'dojo/parser',
    'dijit/layout/ContentPane',
    'dijit/layout/AccordionContainer',
    'dijit/layout/BorderContainer',
    'app.js',
    'dojo/domReady!'
    ], function(parser, ContentPane, AccordionContainer, BorderContainer){
    $(document).ready(function(){
        mapApp = new MapApplication();
        mapApp.prototype = MapApplication.prototype; //dziedzicz

        //access to constructors - no need to load it twice
        mapApp.ContentPane = ContentPane;
        mapApp.AccordionContainer = AccordionContainer;
        //
        parser.parse();
        $.getJSON(mapApp.CONFIG_FILE, function(json){
            if (json.app)       mapApp.initApp(json.app);
            if (json.map)       mapApp.initMap(json.map, json.widgets);
        });
    });
});

Maybe it can be a clue that inside initMap I use another require, but anyway code does not reach this place because I would know (I have debug logs in initApp, initMap and so on).

Another reason I guess may be that I don't use basemaps. I load WMTSLayers and treat them as baselayers (but when loading is done, spatial reference is set properly)

My suspect is parser.parse(), but I have no idea what could be wrong with that.


texts: Fixing pgRouting error "starting vertex not found"

texts:
I've posted on here before about my battles with drivetime polygons and got some sensible suggestions back.
I'm on a Windows 7 PC. I've followed the excellent guide at underdark.wordpress.com and got postgreSQL installed, with PostGIS and pgRouting. I've downloaded OSM data for the UK from Cloudmade (and Geofabrik for that matter) and I've used the wonderful osm2po to get that data into the postgreSQL database.
So far, so good. No errors, no problems, looks like a database to me. I've got a table called network with obvious column names like id, source, target, kmh, cost etc. Just what I would have expected if I'd defined them myself.
So I've fired up the SQL query panel, banged in the SELECT * FROM Driving_Distance statement and fed it the id, source, target and cost columns. I was trying to find 30 minutes around Northampton and the output did indeed some roads within Northampton(shire) - but the rest were in North Wales, which can't possibly be right.
At 3:30am this morning it occured to me that the table also had columns named osm_source_id, osm_target_id and osm_id so perhaps it might be a bit more sensible if I used those. Made sense to me as most things do at that time of the morning, so I ran it again and now I just get "starting vertex not found" every single time.
So to summarise. This works, but doesn't find the right roads:
SELECT network.id,network.x2,network.y2,network.geom_way
 FROM network
 JOIN
 (SELECT * FROM driving_distance('
  SELECT id,
source,
target,
     cost/60::float8 AS "cost"
  FROM network',
  1516905,
  0.025,
  false,
  false)) AS route
 ON
 network.id = route.vertex_id;

This doesn't work at all:
SELECT network.id,network.x2,network.y2,network.geom_way
 FROM network
 JOIN
 (SELECT * FROM driving_distance('
  SELECT osm_id::int4 AS "id",
osm_source_id::int4 AS "source",
osm_target_id::int4 AS "target",
    cost/60::float8 AS "cost"
  FROM network',
  1516905,
  0.025,
  false,
  false)) AS route
 ON
 network.id = route.vertex_id;

The same message happens whether or not I'm using Cloudmade or Geofabrik data, so I'm at a bit of a loss.


texts: Calculating polygon's length along a line?

texts:
I need to calculate the length of polygons along / projected on a line. An image to show the idea:

I have thousands of polygons to be "projected" on lines to verify how many meters of the line they belong to. The polygons are irregular, with different orientation, part of them intersect the line, part of them do not.

I've tried:

Locate features along line tool in ArcGIS Pro, but it only
gives the lenght of line for overlapping part of a polygon (tool
description here ESRI help)
Minimum Bounding Geometry with 'convex hull' option, but it
only gives distance between antipodal points of a polygon, which is
not parallell to the my line (tool description here ESRI help

Am I missing a tool for this in ArcGIS  or QGIS? Would you have a solution for this?


texts: Obtaining land surface temperature in Google Earth Engine using GEET and TOA error is showing

texts:
I am very new to using GEE and only using for one unit I'm currently doing (only using for two months and I have no coding experience). I am trying to calculate the land surface temperature estimate and I found a code from E. Lacerda but the only thing is that I get an error saying
"band_to_toa.expression is not a function
in toa_radiance(), line 938
..."
I've spent past  few hours trying to get my way around this by googling but nothing was useful for me. I'm really not sure what to do at this point and really trying to get this done asap for due dates.
Here's the link to my script
https://code.earthengine.google.com/14a01e61105994465263d252f0c501a1
var raw = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA")
  .filterDate ('2016-01-01', '2016-12-31');

var surf_temp = function(image) {
  var geet = require('users/elacerda/geet:geet');
  var new_toa_radiance = geet.toa_radiance(raw, 11);
  var brightness = geet.brightness_temp_l8c(new_toa_radiance, true);
  var ndvi = geet.ndvi_l8(brightness);
  var veg = geet.prop_veg(ndvi);
  var lse = geet.surface_emissivity(veg);
  var surface = geet.surface_temperature_oli(lse);
  var lst = surface.select('LST');
  var img = lst
  return img;

};
Map.addLayer(raw.map(surf_temp));


texts: Is there a way to have custom functions automatically load in QGIS on start up?

texts:
I have a custom function for the print composer that prints out copyrights, but every time I start QGIS up and go to the Print Composer for the first time, I have to the Insert Expression > Function Editor > Select the function > Ok

This is awkward enough for me, but it's only going to get worse when I roll out QGIS out (to mostly replace another GIS) to my users.


texts: Rendering graduated symbology for dynamically called layers using ArcPy in ArcGIS Pro

texts:
I am trying to render graduated symbology for dynamically called layers using arcpy in ArcGIS Pro 2.7.  The first code block below adds the layers to the map, but when I try to call these layers and customize the symbology I receive errors as included below in the 2nd and 3rd code block.
Is the problem that I am not calling the layers correctly?
#List features from Results feature dataset,
#Makes feature layer from features, adds feature layers to map
m = aprx.listMaps()[0]
datasets = arcpy.ListDatasets(feature_type='feature')
for ds in datasets:
    ResultList= arcpy.ListFeatureClasses("*", "ALL", ds)
    for fc in ResultList:
        OutFile = os.path.join(wkspc, str(fc))
        arcpy.MakeFeatureLayer_management(fc, fc, None, None)
        lf = arcpy.management.SaveToLayerFile(fc, OutFile)

#render graduated symbology for the layers referenced above;
Line 2: for lyr in l():     TypeError: 'layer' object is not callable
l= m.listLayers()[0]
for lyr in l():
  if lyr.isFeatureLayer:
     sym = lyr.symbology
     if hasattr(sym, 'renderer'):
       if sym.renderer.type == "GraduatedSymbolsRenderer":
#modify graduated symbol renderer
        sym.renderer.classificationField = "Total_Deli"
        sym.renderer.breakCount = 6
        sym.renderer.minimumSymbolSize = 10
        sym.renderer.maximumSymbolSize = 25

        lyr.symbology = sym

Different approach--
Line 2: sym = l.symbology     AttributeError: 'list' object has no attribute 'symbology'
    l = m.listLayers()[0]
    sym = l.symbology
    if hasattr(sym,'renderer'):
        if sym.renderer.type == "GraduatedSymbolsRenderer":
            sym.renderer.classificationField = "Total_Deli"
            sym.renderer.breakCount = 6
            sym.renderer.minimumSymbolSize = 10
            sym.renderer.maximumSymbolSize = 25

            l.symbology = sym