texts: TableView has no workspacePath in ArcPy?

texts:
One of my colleagues is running a script which includes the following Python in ArcMap...

oTble = arcpy.mapping.TableView(tbleName)
Workspace = oTble.workspacePath

print "oTble:  '{}'".format(oTble)
if not Workspace:
    print "NULL Workspace!"
    sys.exit()

print "WORKSPACE:  {}".format(Workspace)

...and it exits after producing the following output:

oTble:  '<geoprocessing Standalone Table object object at 0x2B734A88>'
NULL Workspace!

Another colleague runs the same script and it continues on without this problem.

I cannot figure out how a TableView cannot have a NULL workspacePath .

Any ideas how this could happen?

UPDATE 1:

I have created a new MXD, and added the table to this new MXD.  It works fine there.  If I delete ALL layers and tables from the original MXD, and then add the table back into the original MXD in the same way as I did for the new MXD, it still fails.

UPDATE 2:

Based on the suggestion in comments, I've added the following code for debugging purposes...

d = arcpy.Describe(oTble)
print d.baseName
print d.catalogPath

This outputs (actual values modified to protect the innocent):

SDE_DBName.GISADMIN
\SDE_DBName.GISADMIN.TableName

This does not include the workspace in it anywhere.

If I run the same on a new blank document after adding in the relevant table, I get:

SDE_DBName.GISADMIN.TableName
\\host\Connections\ConnectionFileName.sde\SDE_DBName.GISADMIN.TableName


texts: arcpy doesn't save output from "Generate Points Along Lines" to "in_memory"

texts:
I'm trying to run the tool "Generate Points Along Lines" like this:

arcpy.GeneratePointsAlongLines(inputFeature,r"in_memory\temp","DISTANCE","100 meters")

when setting the output as "in_memory", an error occurs:

  000840 the value is not a workspace.

but when I'm running it from the Toolbox it works fine.

Any ideas on what I'm doing wrong or how to make it run using ArcPy?


texts: Experience with issues faced while upgrading from OpenLayers 4 to 6

texts: We are thinking of migrating from OL 4 to 6, skipping 5. I read the upgrade notes, but I wanted to know if someone here also jumped from v4 to v6 and their experience if they faced any major hurdles in upgrading.


texts: Real centroid of a polygon in QGIS 3.2.1? [¿Centroide real de un poligono en QGIS 3.2.1?]

texts:
alguien sabe como obtener el centroide real de un polígono, en la versión 2.18 había una herramienta con el mismo nombre. ahora en la versión 3.2.1 no esta, quiero obtener el centroide real de polígonos tipo peine que quede el punto dentro del polígono, que no le genere un extend. la herramienta que trae de centroide deja el punto fuera del polígono.

Google Translated

  "someone knows how to get the real centroid of a polygon, in version
  2.18 there was a tool with the same name. now in version 3.2.1 is not, I want to get the real centroid of comb polygons that is the point
  inside the polygon, that does not generate an extend. the tool that
  brings centroid leaves the point outside the polygon."


texts: Invalid geometry error when running Vector Intersection tool in QGIS 3

texts:
I am reworking an exercise that runs fine in QGIS 2.14 and 2.18 but have encountered an error in QGIS 3.2. I am using Windows 8.

I am using the Alaska training data and adding the majrivers and trees shapefiles. I select the Susitna River, which works fine and then set a buffer of 5000 around it. Here are the steps I am using in QGIS 3.2 (steps 3,4 and 5
were not needed in QGIS 2.18)

Processing Toolbox/Vector Selection/Select by Attribute (set Description = Susitna River) and it highlights
Vector/Geoprocessing Tools/Buffer ( input is majrivers, check the selected features only and dissolve, distance = 5000)
Processing tools/Vector Geometry/Fix Geometries on the buffer layer
Processing Tools/Vector Geometry/Check Validity
Processing Tools/Vector Geometry/Remove Null Geometries
Vector/Geoprocessing Tools/Intersection (Input layer is fixed geometries, intersection layer is trees)

The following error occurs when running intersection:

Feature (118) has invalid geometry. Please fix the geometry or change the Processing setting to the "Ignore invalid input features" option.
Execution failed after 0.16 seconds

Are there other tools that can correct the invalid geometry for feature 118?


texts: geojson file not interpreted properly (wrong format?)

texts:
My problem:

I am using a web application that uses geojson file to define point data - here is the application:

https://smartcville.com/racks.html

I would like to use the skeleton of the application to make my own maps,
again keeping (different) point data in a geojson file. I typed in the
json file myself, converted it into geojson with a python script, and confirmed with geojsonlint that the data are properly formatted.

Now the PROBLEM:

the geojson file I have put on my own website (just scp-ed it
over). The application does not work any more. I have tried to google for an answer, but, to my surprise, could not find any answers, anywhere.

I figured the following: the original geojson file is here:

var bikeGeoJsonUrl = "https://opendata.arcgis.com/datasets/5fe9482a4ed34551a883355db0f41164_3.geojson";

When I try to load the above geojson file in the browser I get "this is geojson file, what would you like to do, open, save, etc"

My geojson file is on my website:

var bikeGeoJsonUrl  = "https://enjoy-virginia.com/geojson/cville-bike-racks.geojson";

When I try to load it in the browser, it just loads in as a text file,
so browser thinks it is just a text.

I found that I could actually use myjson.com site to store my text
geojson file and open racks.html file with locations properly.

BUT: I would really want to know how is my geojson file on my website different in format from the way myjson.com is saving the same file?

I also remember a few months ago I had my geojson files stored on github and the applications worked fine, but it has now stopped working.


texts: Sample points randomly within raster cells in R

texts:
I would like to generate one random point for each cell in a raster, while excluding NAs. I can use sampleRandom however this just gives me the centroid of each cell:

    library(raster)
    ras <- raster(nrows = 3, ncols = 3)
    v <- c(1,2,NA,4,NA,NA,7,8,9)
    ras[] <- v
    plot(ras)
    samp <- sampleRandom(ras, ncell(ras), xy = TRUE, sp=TRUE, na.rm = TRUE)
    points(samp)

I want to generate random points so that I can then extract data from multiple other raster layers which are not necessarily the same resolution (and therefore the centroid isn't representative so I'd rather have a random point).

I am currently doing it using spsample on an spdf, but I'd rather avoid having to use shapefiles if possible as the memory requirements are getting too big.


texts: Using over() function in R with Polygons and Raster

texts:
I want to calculate the average rain by region using the over function in R but I am getting the error "unable to find an inherited method for function over for signature RasterLayer, SpatialPolygonsDataFrame". This is my code:
rain <- readOGR(dsn="/Users___.shp", layer="rain")
class(rain)
#SpatialPolygonsDataFrame

region <- readOGR(dsn="/Users___.shp", layer="region")
class(region)
#SpatialPolygonsDataFrame

s <- raster(ncol=90, nrow=90)
extent(s) <- extent(rain)
rain_raster <- rasterize(rain, s)
class(rain_raster)
#RasterLayer

test1 <- over(rain_raster, region, fn="mean")
#error: unable to find an inherited method for function ‘over’ for signature ‘"RasterLayer", "SpatialPolygonsDataFrame"

How can I calculate this average without getting this error?


texts: Concatenating Fields in ArcMap dynamically?

texts:
I added a field to an address layer in ArcMap 10 to concatenate separate address fields.  I used the field calculator and a VB Script to accomplish this.  So, for example, if the separate fields are "10" "Main" and "Street", then the new field would be "10 Main Street".

The problem is that the concatenated field needs to update dynamically.  So if the user changes "Main" in the database to "East Main" then that field will update in ArcMap.  However, the concatenated field will not update to "10 East Main Street" unless I rerun the Field Calculator script.

I need a way for this new field to always reflect the current data, but still do the concatenation.  Does anyone have any suggestions?


texts: how to draw route from start point to destination point with intermediate points on Open Source Routing Machine?

texts:
I have four points like this:
  Lat            Log

43.532284     -6.060434

43.533880     -6.065279

43.534269     -6.067436

43.533484     -6.071405

If I put these points in Open Source Routing Machine I obtain this map:

I would like to obtain the map, but with a number in each point intermediate like this:

Any idea?


texts: Can script using arcpy.mapping.MoveLayer() be made to run faster?

texts:
We're running ArcGIS 10.0 and using a python script to move layers, turn on the correct layers, zoom to the correct extent and export a pdf map.  The script normally takes 30 - 35 seconds to finish.  Moving layers takes the bulk of the processing time, usually 20 - 25 seconds using this code:

for y in arcpy.mapping.ListLayers(mxd):
    move_lyr = arcpy.mapping.Layer('Available Zips')
refLayer = move_lyr
for lyr in arcpy.mapping.ListLayers(mxd):
    if lyr.name == zip_view:
        zip_move = lyr
for lyr in arcpy.mapping.ListLayers(mxd):
    if lyr.name == zip_view:
        arcpy.mapping.MoveLayer(df, refLayer, zip_move, "BEFORE")

Is there a way to optimize this code to run faster?


texts: Colorize radar imagery so it looks like optical imagery

texts:
Here is an imagery from a radar satellite

Is it possible to colorize a radar image, so it looks like and enables the same kinds of analysis as for an image captured from an optical satellite (simple RGB imagery)?

Since radars can capture imagery even through clouds and with no visible light on the Earth, it seems I have access to a much larger archive of imagery if I can just give the "same properties" to radar imagery.


texts: Copy/ Paste from spatial attribute tables into non-spatial tables in Arcmap

texts:
Imagine I have a shapefile which has its own attribute table and I have a non-spatial stand alone table too. I am going to copy the rows from the shapefile attribute table into the stand alone table but when I start editing on the table and copy the features from the shapefile attribute table, paste doesn't work and I can't copy the values inside the stand-alone table.
I was wondering if there is a way to copy and paste only the values of spatial attribute tables rows (and obviously not the spatial data) into a non-spatial table in Arcmap.


texts:
Does anyone know how to interface a Garmin GPS 78 Marine GPS with their computer to manage waypoints (upload, download, save, etc.) Previously I was using the MN DNR Garmin tool, but it no longer recognizes the GPS.
 Any solutions are welcome
 Thanks


texts: I am currently using Mapbox GL JS and pull multiple GeoJson files from the style.json file. Each file in about 10MB and takes too long to load. Is there a way to compress the files and uncompressed them in Mapbox?


texts: Editing and Deleting Query Tables in ArcMap

texts:
How do I edit the query for a Query Table after it has been created?

How do I permanently delete a Query Table? Even after I do right-click -> remove, it still appears in dropdown lists.


texts: Eliminating unused libraries from standalone Python-QGIS application package by PyInstaller?

texts:
We have a standalone PyQGis windows application that is packaged into a multi-file folder by PyInstaller (3.2.1). The OSGeo4W64 python environment is used in the application and the spec file is very similar to the one posted in this other question thread:

Working Windows Standalone .EXE From QGIS Python Script Using PyInstaller?

The resulting package works fine, except its size is very large. In particular, it automatically pulled in libraries such as the following, among other things:

LIBEAY32.dll
libifcoremd.dll
libifportmd.dll
libiomp5md.dll
libiompstubs5md.dll
libmmd.dll
libmysql.dll
libpng16.dll
LIBPQ.dll
libsodium.pyd
libtiff.dll
wxbase28uh_net_vc.dll
wxbase28uh_vc.dll
wxmsw28uh_adv_vc.dll
wxmsw28uh_core_vc.dll
wxmsw28uh_html_vc.dll
mfc90.dll
mfc90u.dll
mfcm90.dll
mfcm90u.dll
Microsoft.VC90.CRT.manifest
Microsoft.VC90.MFC.manifest
msvcm90.dll
MSVCP100.dll
msvcp90.dll
MSVCR100.dll
msvcr90.dll

We removed these libraries to reduce the package size, but how to prevent the Pyinstaller from including these and other non-essential libraries or binaries in the first place? I think one way is to explicitly exclude these using the 'exclude' clause, but I am hoping there is a more systematic way to do this, for example, is there a way to signal to PyInstaller to only go after the first level of dependencies (not to include others recursively).


texts:
I "upgraded" from 3.2.3 to 3.4.0 and it's so unresponsive it's unusable. The startup time is the same as 3.2.3, then it needs half a minute to build its UI.

When I click on the UI, the reaction is the same as if the computer is under 100% load; the task manager however does not report any spikes. After another half minute of being the active window, the UI is responsive (can be dragged, smooth zoom, hovering over the top context menu pops up the dropdowns, etc.) until something is done that makes the window lose focus again, e. g. dragging a file from the explorer, or clicking on a layer to access its properties window. Is anyone else experiencing this problem? Version 3.2.3 runs buttery smooth.

// edit:

I'm running Windows 10, 3.1 GHz hexacore, 16 GB RAM, 8GB video card. As suggested below, I tried reinstalling, as well as the --noplugins flag. Neither led to any improvement.


texts: Calculating wait time on one-lane road with two-way travel using ArcGIS Network Analyst?

texts:
I want to know how to calculate the wait time on a single-track/one-lane road (only room for one car at a time) for one car waiting for another car to pass by.

For example, two cars drive on the one-lane road in opposition to each other. When the two cars meet in the road, one car should wait for the other car to pass.

I wonder if Network Analyst (especially, Vehicle Routing Problem) can solve this problem?


texts: How to extract a WCS layer with gdal_translate?

texts:
There is a certain WCS layer that I'd like to save it to my local computer for analysis.

According to  GDAL WCS driver  draft:

"Accessing a WCS server is accomplished by creating a local service description xml file looking something like the following, with the coverage server url, and the name of the coverage to access. It is important that there be no spaces or other content before the  element.

<WCS_GDAL>
   <ServiceURL>http://laits.gmu.edu/cgi-bin/NWGISS/NWGISS?</ServiceURL>
   <CoverageName>AUTUMN.hdf</CoverageName>
</WCS_GDAL>

I tried to create the xml file but I was unsuccessful.

Does anyone has any experience about extracting WCS  layers?
Are there any pointers (sites, tutorials, tips) on how to create the xml file?

EDIT:

The xml file im using to grab the first layer for a specific time is the following:

<WCS_GDAL>
   <ServiceURL>http://dmcsee.org/cgi-bin/mapserv?map=/var/www/tmp/dmcsee_wms/dmcsee_wcs.map</ServiceURL>
   <Version>1.1.1</Version>
   <CoverageName>SPI6</CoverageName>
   <supportedCRSs>EPSG:4326</supportedCRSs>
   <PreferredFormat>gtiff</PreferredFormat>
   <GetCoverageExtra>&amp;BBOX=10.0,31.995,48.005,50.0&amp;time=19900101</GetCoverageExtra>
</WCS_GDAL>

But for somereason returns an error


texts:
I am using ArcMap 10.8 to export Data (from the Table of Contents menu). So I chose - Data - Export Data. I then change my location from a gdb to a folder, so I can have the ability to chose the TIFF format.

But, when press save I get an error: Failed to save raster dataset, Failed to create raster data set.
Is there something here I am missing?


texts: Accessing multiple shapefiles in folder using Python and ArcGIS for Desktop?

texts:
I have been asked to create a custom script to reproject shapefiles (different projections) from a folder, but can't contain values such as filenames and path names to access them.

I have started with Batch Project script and tried to set the input parameter as 'folder' but get the error 000840: The value is not a dataset, because the input values needs to be feature class or dataset.

I have also searched the forum and found:

How do I use python to batch process the reprojection of a group of shapefiles?

But obviously this answer contains values. I have also tried setting;

Set workspace environment;

arcpy.env.workspace = arcpy.GetParameterAsText(0)

Doesn't appear to work.

It is now working after I changed my code from:

arcpy.env.workspace = arcpy.GetParameterAsText(0)

to:

arcpy.env.workspace = arcpy.ListFeatureClasses()


texts: Adding coordinates to polygon layer to sort polygons from north to south in QGIS?

texts:
I have a polygon file, where each polygon represents one property.
I want to be able to sort the polygons from north to south or vice versa. I need them to be sorted in this way when I export them to a different program. In the other program I want the properties to appear in a north to south order when I visit them one by one.

My first problem is, the polygons don't have information about there coordinates in the attribute table.

How can I add this?

I have tried using the field calculator, but I don't seem to get it right.

And, next step, if I sort them in QGIS, can I somehow "save" that order so that they appear in the same order in the other program?

I have tried the method in the answer below, but I don't get it to work on my data. Then after almost giving up, I read in the help of the field calculator that $y gives the y coordinate of a point.

But I have a polygon layer?

What do I do?


texts:
I am trying to use the Landscape Vector Overlay function in the LecoS plugin to calculate the amount of various proportions of land cover types (using an ordinary landscape raster) within buffer polygons around individual points (i.e., Landscape Proportions).

When I try to run this I get the following error:

An error has occurred while executing Python code:
AttributeError: 'NoneType' object has no attribute 'GetFeatureCount'

I have little to no Python experience so I'm hoping my issue does not involve coding. However, I have not installed anything in Python (because I don't know how) and I realize that might be involved.

I am using QGIS 3.4

I am a pretty basic QGIS user.

I have researched several other options but always end up at a dead end. The LecoS plugin seems like it is designed specifically for what I am trying to accomplish.


texts: pyshp problem with naturalearthdata-shapefile

texts:
trying to follow this example:
https://code.google.com/p/pyshp/wiki/ShapeRecords

entering the following code:

import shapefile
r = shapefile.Reader("ne_110m_admin_0_countries")
sr = r.shapeRecords()

gives this error-message under python 3,3 64 bit on windows 7:

Traceback (most recent call last):
  File "D:\stuff\beispiel1.py", line 3, in <module>
    sr = r.shapeRecords()
  File "C:\Python33\lib\site-packages\shapefile.py", line 553, in shapeRecords
    for rec in zip(self.shapes(), self.records())]
  File "C:\Python33\lib\site-packages\shapefile.py", line 525, in records
    r = self.__record()
  File "C:\Python33\lib\site-packages\shapefile.py", line 501, in __record
    value = u(value)
  File "C:\Python33\lib\site-packages\shapefile.py", line 60, in u
    return v.decode('utf-8')
UnicodeDecodeError: 'utf-8' codec can't decode byte 0xf4 in position 1: invalid continuation byte

Any ideas, what I am doing wrong (just used another shapefile from www.naturalearthdata.com) ???

UPDATE (May, 8th. 2014)
Today I received an email from the developer of pyshp, Mr. Joel Lawhead. The problem I mentioned was a bug, which only occurs, when you use python vers 3.x . python vers. 2.x was not affected. He kindly send me a patched version, which works fine for me. He also told me, that this new version will be released soon.


texts: QGIS Polygons missing after saving edits

texts:
Right, I am slightly worried. I opened QGIS this morning to discover that part of the work I did yesterday was not there (the attributes where present in the table, but were not digitised anywhere on the interface). I had assumed that I had accidentally done something last night before I went home that caused this to happen.

So I thought, not to worry, it won't take long to re-digitise the lost polygons. I finished the work, produced the map and saved the project so that I could close QGIS and was very worried to see that of my 13 polygons, only four are left! This has happened simultaneously with two layers (the two layers show the same boundary but the polygons are subdivided in one of the layers).

Is there any way I can undo, this/fix this. Whilst I should have finished with the layer, it is not beyond the realms of possibility that I will be asked to make minor changes etc.


texts: Webmap app using geoserver must have source code open?

texts:
We have a Webmap that uses geoserver + geoExt + openLayers. Even though the JavaScript code is not obfuscated or compressed, at the moment we do not intend to have a place were someone can got and download our website "source code". A friend said that since we use GPL software, we must make the source code available. Now, I believe the use of geoserver, geoExt, openLayers (with its licenses GPL, FreeBSD, BSD) does not enforce me to distribute our webmap source code, is that correct? Or am I understanding GPL wrong?

I'm basing my understanding in two sentences, one from GPL FAQ and another from Drupal FAQ.

  From GPL FAQ: If I know someone has a copy of a GPL-covered program,
  can I demand he give me a copy? No. The GPL gives him permission to
  make and redistribute copies of the program if he chooses to do so. He
  also has the right not to redistribute the program, if that is what he
  chooses.

  From Drupal FAQ: 11: Do I have to give the code for my web site to
  anyone who visits it?

  No. The GPL does not consider viewing a web site to count as
  "distributing", so you are not required to share the code running on
  your server.


texts: How to calculate centroid of a polygon defined by a list of longitude/latitude points?

texts:
I've imported some polygons (ways) from OpenStreetMap and now want to calculate the centroids of them.

For this purpose I wrote following code:

    final List<Point2D.Double> points = getPoints();
    final Point2D.Double actualCentroid = calculateCentroidOurWay(points);

    LOGGER.debug("Centroid, x: " + actualCentroid.x);
    LOGGER.debug("Centroid, y: " + actualCentroid.y);

private Point.Double calculateCentroidOurWay(final List<Point2D.Double> aPoints) {
    final Building objectUnderTest = new Building();

    objectUnderTest.setPolygon(new Polygon(aPoints));

    return objectUnderTest.calculateCentroid();
}

@Override
public Point2D.Double calculateCentroid() {
    double x = 0.;
    double y = 0.;
    double area = 0.;

    for (int i = 0; i < polygon.getPoints().size() - 1; i++) {
        final Point2D.Double point = polygon.getPoints().get(i);
        final Point2D.Double pointN = polygon.getPoints().get(i + 1);

        final double temp = point.getX() * pointN.getY() - pointN.getX() * point
                .getY();
        x += (point.getX() + pointN.getX()) * temp;

        y += (point.getY() + pointN.getY()) * temp;

        area += temp;
    }

    area /= 2.;

    x *= 1 / 6. * area;
    y *= 1 / 6. * area;

    return new Point2D.Double(x, y);
}

private List<Point.Double> getPoints() {
    final List<Point.Double> points = new LinkedList<>();

    points.add(new Point.Double(132.1764353, 43.3539477));
    points.add(new Point.Double(132.1771393, 43.3539409));
    points.add(new Point.Double(132.177145, 43.3542587));
    points.add(new Point.Double(132.1768224, 43.3542618));
    points.add(new Point.Double(132.1768183, 43.3540332));
    points.add(new Point.Double(132.1766516, 43.3540348));
    points.add(new Point.Double(132.1766564, 43.3542949));
    points.add(new Point.Double(132.176566, 43.3542957));
    points.add(new Point.Double(132.1765705, 43.3545391));
    points.add(new Point.Double(132.1769718, 43.3545353));
    points.add(new Point.Double(132.1769691, 43.3543859));
    points.add(new Point.Double(132.177131, 43.3543843));
    points.add(new Point.Double(132.1771321, 43.3544446));
    points.add(new Point.Double(132.177341, 43.3544426));
    points.add(new Point.Double(132.1773296, 43.3538152));
    points.add(new Point.Double(132.176433, 43.3538238));
    points.add(new Point.Double(132.1764353, 43.3539477));
    return points;
}

The result of the calculation is

Centroid, x: 1.8787362164666136E-11
Centroid, y: 6.162273026564668E-12

which is obviously wrong (it's not a latitude/longitude pair).

When I use Java Topology Suite, it returns same results.

private Coordinate calculateJtsCentroid(final List<Point2D.Double> aPoints) {
    final CentroidPoint centroidPoint = new CentroidPoint();

    for (final Point.Double curPoint : aPoints) {
        centroidPoint.add(new Coordinate(curPoint.x, curPoint.y));
    }

    return centroidPoint.getCentroid();
}

How should I change my code in order to get the centroid of the polygon as a point identified by latitude and longitude?


texts: Click a marker icon to open a pop-up, using OpenLayers

texts:
Background

Looking to create a pop-up that performs an Ajax request when the user clicks a marker on an OpenLayers map. The markers, however, have different icons.

The following is working code that can be used to accomplish the desired task:

var sites new OpenLayers.Layer.Markers( 'Sites' );
map.addLayer( sites );
addSiteMarkers( sites );

function addMarker( layer, lat, lng, img ) {
  var size   = new OpenLayers.Size( 32, 37 );
  var offset = new OpenLayers.Pixel( -(size.w/2), -size.h );
  var icon   = new OpenLayers.Icon( 'img/marker/' + img, size, offset );
  var marker = new OpenLayers.Marker( getLocation( lat, lng ), icon );

  layer.addMarker( marker );

  marker.events.register( 'click', layer, clickMarker );
}

function clickMarker() {
  alert( 'click' );

  // Pop-up code will go here...
}

function addSiteMarkers( sites ) {
  addMarker( sites, -113.3017895, 49.1939798, 'store.png' );
  addMarker( sites, -113.3072538, 49.2013626, 'building.png' );
}

The above code adds two markers (of hundreds) each having a different icon (from a set of fifty different icons). Each click will display a pop-up and each pop-up will have information that is retrieved from a database (using a jQuery ajax call).

Problem

I suspect this is a terribly inefficient way to use OpenLayers and map markers for opening pop-ups when a user clicks an icon.

Question

How do you use OpenLayers to present a variety of different markers that when clicked display a pop-up (with content dependent on the actual marker that was clicked)?

Update #1

Most of the examples show how to use a Vector layer, rather than a Marker layer, to handle click events. Using a Vector layer makes making markers more mind-boggling.

Related

http://dev.openlayers.org/examples/marker-shadow.html
https://stackoverflow.com/questions/7702945/adding-popup-info-in-openlayers
https://stackoverflow.com/questions/9895951/how-to-display-multiple-markers-with-individual-framecloud-popups-in-openlayers
https://stackoverflow.com/questions/6026420/openlayers-vector-features-instead-of-markers
http://dev.openlayers.org/docs/files/OpenLayers/Control/SelectFeature-js.html
http://docs.openlayers.org/library/overlays.html
http://dev.openlayers.org/docs/files/OpenLayers/Popup-js.html
http://dev.openlayers.org/docs/files/OpenLayers/Marker-js.html

Thank you!


texts:
I'm working with an attribute table in ArcMap 10.2 where there are some 30 fields mainly consisting of NULL values. I don't want to display to users a table full of NULL values as it'll be difficult to read. See below image:

So the result I'd like to create is a Field called SUMMARY that includes the "field name: Value. field name: Value" of any row that is NOT NULL e.g. A: T. D: 10/01/2011

In this way a user can select a feature with the identify tool and see what values there are rather than scrolling through loads of NULL values.

Does anyone have suggestion on the general construct of the script? Including the da.SearchCursor, da.UpdateCursor, and I guess by way of row index if statements to query if not none add to SUMMARY field with delimators ":" after field name and "." after row value.


texts: How to make a clickable label in qgis from an excel sheet?

texts:
I have location data in excel sheet where each label has corresponding data file for each financial year.
I want it that when i click on a label, years pop up, and when i click on that year, its corresponding file open up.
I prefer to keep my data in excel.
Any suggestions on how to do it?
or any other easier way/solution to process my data is also welcome.


texts: Mosaic to new raster: same number of bands issue

texts:
I am trying to mosaic 32 raster files (min value 0, max value 13) into one raster file so that I can later use this new raster to extract multiple values to points. However, when I try and mosaic the rasters, this error comes up:

<<Input and target dataset should have the same number of bands>>.

I am new to ArcGIS so I am not sure what this means - especially as I do not see a target dataset? Also, I am not sure what the band value should be? When I investigated the data source, it said:

  Reference composite imagery are median observations from a set of
  quality assessed growing season observations in four spectral bands,
  specifically Landsat bands 3, 4, 5, and 7

Therefore logically I would assume that the band value is four? I just don't understand why this isn't working. I am using ArcMap 10.2.


texts: Automatically exporting multiple maps with quantile symbology in ArcGIS Desktop?

texts:
I have a very simple shapefile to which I join a bunch of fields. I then need to have these new fields symbolized with a specific color ramp using a five quantile classification and then exported to a JPG at a certain resolution with a specific filename.

Is it possible to sort of "preload" the software manually with the layers and page layout I want and then have it loop in such a way that it merely changes the symbology field value and exports, over and over?

There is a legend but since the labels will remain the same (quantile 1, quantile 2, etc) I can just leave that as a hardcoded graphic.

I'm envisioning something like this: I add my two layers (one is a fill-free layer of US states to serve as an outline for the other layer) and then join one big file to the layer that needs symbology. This program I don't know how to write would then loop through and change the symbology field value to the next column I've added, and changes the classification back to a 5-class quantile (it always reverts to Jenks for some reason). Then since the legend is hardcoded it just needs to save this new map of field Beta1 as a jpg and then switch to Beta2, Beta3, etc.


texts: Debug QgsDataSourceUri as a QgsVectorLayer

texts:
I'm trying to figure out what is wrong with my code, sometimes it works sometimes it doesn't, I don't know how to debug it when it doesn't work.

Basically I'm working with this code:

uri = QgsDataSourceUri()
uri.setConnection(str(host), str(port), str(dbname), str(username), str(password))
uri.setDataSource(str(schema), str(table), str(geom_col), '', 'field_row_id')
uri.setKeyColumn('field_row_id')
uri.setSrid('4326')
vlayer = QgsVectorLayer(uri.uri(), str(extra_name) + str(table), 'postgres')
print(uri.uri())

Prints something like this: 'dbname=\'hortegarden\' host=geodatafarm.com port=5432 user=\'my_username\' password=\'my_password\' key=\'field_row_id\' srid=4326 table="harvest"."dyback_2015" (pos) sql='

field_row_id is a valid primary key with unique values.
I can add the table from the source manager, hence I want it to be apart of my plugin. Though the source from the working layer properties looks like this:
dbname='hortegarden' host=geodatafarm.com port=5432 sslmode=disable authcfg=0xq253x key='field_row_id' srid=4326 type=Point table="harvest"."dyback_2015" (pos) sql=

Any suggestion how I can investigate why it is not valid?


texts:
I have some contours where the z value is stored in the geometry and I want to display it as an attribute. I want to include this in a python script.

I have an expression which does this in python field calculator:

!Shape!.firstpoint.z

This works perfectly, but when I include it in my script:

arcpy.CalculateField_management( shapefile , "VALUE" , "!Shape!.firstpoint.z", "PYTHON")

this does not work, which surprises me because I thought that the Calulate Field  was exactly the same as field calcultor, but accessed through arcpy.

The error is that its a string, so has no z attribute: "  File "", line 1, in
AttributeError: 'str' object has no attribute 'z"

So I've tried:

arcpy.CalculateField_management( shapefile , "VALUE" , "!Shape.firstpoint.z!", "PYTHON")

and that tells me that "firstpoint.z" isnt an attribute of shape

I then tried:

arcpy.CalculateField_management( shapefile , "VALUE" , "!Shape.firstpoint!.z", "PYTHON")

And that seemed to view the first point as a unicode rather than a point. "File "", line 1, in  AttributeError: 'unicode' object has no attribute 'z'"

It seems that my problem is getting it to view firstpoint as a point, but I'm not not why this is.


texts: In a polyline shapefile, how can I color a label using it's corresponding symbol color? (QGIS)

texts: Maybe a pretty simple question. I have a polyline shapefile that I've symbolized through categories, and I wish for the labels to automatically be set to the same color as it's corresponding line segment. Therefore, if I end up changing the symbology in the future the labels will update accordingly. Is this possible, maybe through the expression feature in the labeling options? (I'm very new to qgis).


texts: How to create a virtual field in QGIS by Python Console?

texts: Is there any example about how to create a virtual field and subsequently change the attribute of a non-virtual field without wiping out the virtual field? Both are required to be done in the Python Console.


texts: How to pass 2 parameters in wps with jiffle in GeoServer?

texts:
I'm trying to send 2 sources to process with jiffle (adding a second input on wps:inputs, in this example both are the same layer so the result should be 0), but when I try to execute a simple task like: dest = src-src1, the result is distinct of zero

I'm not sure if it's possible to do like this, because I didn't find any similar example.

<?xml version="1.0" encoding="UTF-8"?><wps:Execute version="1.0.0" service="WPS" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns="http://www.opengis.net/wps/1.0.0" xmlns:wfs="http://www.opengis.net/wfs" xmlns:wps="http://www.opengis.net/wps/1.0.0" xmlns:ows="http://www.opengis.net/ows/1.1" xmlns:gml="http://www.opengis.net/gml" xmlns:ogc="http://www.opengis.net/ogc" xmlns:wcs="http://www.opengis.net/wcs/1.1.1" xmlns:xlink="http://www.w3.org/1999/xlink" xsi:schemaLocation="http://www.opengis.net/wps/1.0.0 http://schemas.opengis.net/wps/1.0.0/wpsAll.xsd">
  <ows:Identifier>ras:Jiffle</ows:Identifier>
  <wps:DataInputs>
    <wps:Input>
      <ows:Identifier>coverage</ows:Identifier>
      <wps:Reference mimeType="image/tiff" xlink:href="http://geoserver/wcs" method="POST">
        <wps:Body>
          <wcs:GetCoverage service="WCS" version="1.1.1">
            <ows:Identifier>cliente1:chanar_agua</ows:Identifier>
            <wcs:DomainSubset>
              <ows:BoundingBox crs="http://www.opengis.net/gml/srs/epsg.xml#32721">
                <ows:LowerCorner>338900.0 6533020.0</ows:LowerCorner>
                <ows:UpperCorner>356890.0 6542140.0</ows:UpperCorner>
              </ows:BoundingBox>
            </wcs:DomainSubset>
            <wcs:Output format="image/tiff"/>
          </wcs:GetCoverage>
        </wps:Body>
      </wps:Reference>
    </wps:Input>
     <wps:Input>
      <ows:Identifier>coverage</ows:Identifier>
      <wps:Reference mimeType="image/tiff" xlink:href="http://geoserver/wcs" method="POST">
        <wps:Body>
          <wcs:GetCoverage service="WCS" version="1.1.1">
            <ows:Identifier>cliente1:chanar_agua</ows:Identifier>
            <wcs:DomainSubset>
              <ows:BoundingBox crs="http://www.opengis.net/gml/srs/epsg.xml#32721">
                <ows:LowerCorner>338900.0 6533020.0</ows:LowerCorner>
                <ows:UpperCorner>356890.0 6542140.0</ows:UpperCorner>
              </ows:BoundingBox>
            </wcs:DomainSubset>
            <wcs:Output format="image/tiff"/>
          </wcs:GetCoverage>
        </wps:Body>
      </wps:Reference>
    </wps:Input>
    <wps:Input>
      <ows:Identifier>script</ows:Identifier>
      <wps:Data>
        <wps:LiteralData>dest=src-src1;</wps:LiteralData>
      </wps:Data>
    </wps:Input>
  </wps:DataInputs>
  <wps:ResponseForm>
    <wps:RawDataOutput mimeType="image/tiff">
      <ows:Identifier>result</ows:Identifier>
    </wps:RawDataOutput>
  </wps:ResponseForm>
</wps:Execute>


texts: QGIS Buffer distance vs Distance in Distance Matrix

texts:
I created Buffer (with Vector, Geoprocessing Tools) with 500 M distance from a point in the map (let's say Point A). I made another point (Point B) and it is not inside the 500 m Buffer.
The Distance matrix tool shows that the points are separated about 413 M from each other.
I noticed that this difference comes (using the Measure Line tool), because the buffer has Cartesian distance, yet the distance in Distance Matrix shows the ellipsoidal distance values.

How to either make the buffer in ellipsoidal distance, or let the distance matrix shows Cartesian distance?

Note: The Cartesian distance of point A and B is about 661 Meters for each other.


texts: ESRI imagery service that is WMS Compliant

texts:
I want to use the ESRI free online imagery service

http://services.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer

But I want a service that is WMS compliant.
Is there a version of this service that is WMS compliant?


texts:
I would like to make the sum of 12 different rasters to get one raster as output.
I wrote the following code but it doesn't work (or better, it works, but at the end it doesn't give me any layer, but the writing "Invalid layer, Gdal provider cannot open Gdal dataset"), any suggestions to how modifying it?
import qgis
import gdal
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
from qgis.core import QgsRasterLayer

suffix_input=['I_jan', 'I_feb', 'I_mar', 'I_apr', 'I_may', 'I_jun', 'I_jul', 'I_aug', 'I_sep', 'I_oct', 'I_nov' ,'I_dec']
inputpath = '/Users/macbook/Desktop/TESI/PROGETTO/I/'

list=[0,1,2,3,4,5,6,7,8,9,10,11]
sum=0
for i in list:
  inputrasterfile = QgsRasterLayer(inputpath + suffix_input[i] + ".tif")

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

  sum = 'sum + ras@1'

  calc = QgsRasterCalculator( sum, '/Users/macbook/Desktop/TESI/CAPITANATA PROGETTO/Itot.tif', 'GTiff', inputrasterfile.extent(), inputrasterfile.width(), inputrasterfile.height(), entries )
  print(calc.processCalculation())
  iface.addRasterLayer('/Users/macbook/Desktop/TESI/PROGETTO/Itot.tif')


texts: Create not overlapping (Thiessen) polygons in service area analysis in ArcGIS?

texts:
Using the network analyst extension, I'm trying to create not overlapping service area polygons around a point layer in ArcMap 10.6.1. Here are the selected settings for polygon generation:

The rest of the settings are default, with the impedenace being TravelTime(minutes) and the default breaks being 60, direction away from facility, and restrictions only for 'Driving an Automobile'.
Nonetheless, the results look like this:

I've also attempted to run the analysis with the Trim Polygon settings, but the polygons still appear to overlap. Can anyone point out a possible reason or a workaround?

EDIT: The polygons seen above do not overlap (it only "seems" they do, but the boundaries indicate distinct features indeed). However, the problem is rather in the incorrect allocation of polygons to points (as also noted in the picture). The result should be reminiscent to a Thiessen polygon analysis - and it isn't.


texts: Do the atmospheric conditions affect all spectral bands equally?

texts:
I would like to know if the atmospheric conditions affect equally the spectral bands of sensors, such as Sentinel-2?
Although the structure of "correctly" analyzing satellite imagery include correcting at-sensor reflectance to bottom-of-atmosphere reflectance, for analysis using spectral band ratios such as NDVI change, does one still need to correct for atmospheric conditions?


texts: Finding the boundary polygon from a set of gridded points (no hole)

texts:
Let us consider a set of a 2-D points (longitude, latitude) each of which is centre of square grid (intersect point of the diagonals of the square).

I try to make you understand what I want to do in a graphical manner.

Sample data (point set):

lon<-c(88.56630, 88.62501, 88.60013, 88.57499, 88.65879, 88.63392, 88.60879)

lat<- c(21.03517, 21.01287, 21.05434, 21.09610, 21.03207, 21.07354, 21.11531)

Plots:

Figure 1: First one is the scatter plot of the point set.

Figure 2: All the points are bilinear (this is known).

Figure 3: Each point is centre of a square of size (approx.) 4km by 4km (we need construct those squares from the given centre point).

Figure 4: This is what I want to have, The polygon obtained by the outer sides of the squares (black lines).

Note: I have lots of point sets of arbitrary shapes like the above sample with maximum cardinality of the sets as 100.

How can this get this polygon (no hole)?


texts: Difference between gdal_calc and QGIS raster calculator

texts:
I am hoping someone can help me with an issue I am having with gdal_calc. I am looking to create NDVI rasters from 4 band 8 bit ortho tiles. I have processed the tiles using the following calculation:

 gdal_calc -A DDE1.tif --A_band=4 -B DDE1.tif --B_band=1 --outfile=DDE1_ndvi.tif --type=Float32 --calc="(A-B)/(A+B)"

The results from the calculation produces integer values that are not typical of NDVI. However, when I use the same data in QGIS raster calculator using the expression below a different, more typical result is created.

 ("DDE1@4"-"DDE1@1")/("DDE1@4"+"DDE1@1")

Is there a different setting or function I need to call in gdal_calc to get values that have decimals? I thought the --type=Float32 would ensure output values as decimals, but they are integers. Below is a graphic with the two results.


texts: Create points along line and apply a 90° offset to them - pyqgis

texts:
I want to create points along a line (works!) and then additionally move the points by a certain offset (perpendicular to the line) with pyqgis 3.   See this sketch for further illustration:

My current python script does following: It basically interpolates the position at the certain length of the line and then it constructs the position of the offset point by triangulation. The script works so far, but the drawn point has a slight deviation in it's x coordinate (see figure 2).

My code so far:

def createPointsAtOffset(distance, offset, geom):
length = geom.length()
currentdistance = distance
feats = []

point = geom.interpolate(currentdistance)
fet = QgsFeature()
fet.setAttributes([currentdistance, offset])

if offset is not 0:
    xy = point.asPoint()
    alpha = (math.degrees(geom.interpolateAngle(currentdistance)) - 90) # return in degree
    # create delta x and y via triangulating
    delY = math.cos(math.radians(alpha)) * offset
    delX = math.sin(math.radians(alpha))* offset #12
    delT = math.sqrt(offset*offset - delY*delY)

    pointX = xy[0] + delX
    pointY = xy[1] + delY
    fet.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(pointX,pointY)))

else:
    fet.setGeometry(point)

feats.append(fet)

return feats

def pointsAlongLineOffset(distance, offset):
# Create a new memory layer and add a distance attribute
vl = QgsVectorLayer("Point", "distance offset nodes", "memory")
pr = vl.dataProvider()
pr.addAttributes( [ QgsField("distance", QVariant.Double), QgsField("offset", QVariant.Double) ] )
layer = iface.mapCanvas().currentLayer()
# Loop though all the features
for feature in layer.getFeatures():
    geom = feature.geometry()
    # hier muss dann eine schleife durch alles distance offset paare durchgehen
    #features = createPointsAt(distance, geom)
    features = createPointsAtOffset(distance, offset, geom)
    pr.addFeatures(features)
    vl.updateExtents()

QgsProject.instance().addMapLayer(vl)

# example for the image above:
pointsAlongLineOffset(1.2,2.5)

I think my error is that I use a wrong angle in the triangulation part. In the example the line is tilted about 6° from east but when I insert about 12° into the calculation of delX it works? Am I missing some kind of deviation from north or is my approach completely wrong?


texts: Using Custom OSM tiles with Mapbox and/or CartoDB?

texts:
For a newspaper, we would like to use CartoDB's editor to produce quick overview maps with small amounts of geometry and information such as tooltips etc. For this job, the CartoDB editor seems to be perfect, and suits our needs better than Mapbox' editor.

However, we would also like to have a custom styled base map for our "corporate identity", just as ZEIT Online recently launched it: https://www.mapbox.com/blog/zeit-launches-custom-maps/

Hosting our own, custom styled OSM tiles is not really an option, especially because we would like to have a highly detailed map of the whole globe.

Still a bit confused by the possibilities of both service providers, we don't know how to proceed: We want to use CartoDB so our journalists can produce quick and good looking maps, but we also want our custom styled OSM tiles - and CartoDB won't host these for us, right?

So we would have to use Mapbox where we would upload our tiles and use their tiles in CartoDB (the editor allows to specify a Mapbox tileserver)?

But then again, where do we get the OSM data from, that we would style in Tilemill?

Do we really have to download the whole planet, and how is it possible to upload such huge data back to Mapbox?

I am quite confused about the workflow that would be best used here.


texts:
I am using QGIS 3.8 and I am trying to merge three rasters according to the scheme pictured on the attached screenshot. I tried many syntax expressions for the Raster calculator but I was unable to get the desired result (merged raster).

Does anybody have any idea how to merge the rasters? Each raster is categorized to 4 groups based on pixel values. The range of groups is same for all rasters - 4 = red color, 3 = blue color, 2 = brownish color, 1 = yellow color (not present on the screenshot).


texts:
Is it possible to extract all 3d data (terrain + objects: buildings, vegetation, landmarks, etc.) from Google Earth in an easy way and 3d print it? I found tutorials about extracting the terrain itself, but I would like to print a section of a city which is nicely populated with detailed models.
Since it would still cover a larger area, downloading and printing the buildings separately is not an option.


texts: SQL for Select By Attribute to only select Feature Class Polygons based on values from a Table

texts:
I have a polygon feature class with one of the attributes being Parcel ID Numbers (Parcels_in_City).
I have a table with only Parcel ID Numbers (Parcels_with_Trees).
Is there a SQL query that can be used in the Select By Attribute tool that will select all of the Parcel ID Numbers from Parcels_in_City based on the values in the Parcels_with_Trees table?


texts: Configuring CodeBlocks to read shapefile using GDAL library in c++ in console application

texts:
I am using CodeBlocks on Windows10 to load vector data using console application in c++ using GDAL/OGR library,  settings are set to as follows:

compiler settings/
    search directories/
        compiler/
            C:\GDAL\bld\include (contains .h files)
        linker/
            C:\\GDAL\bld\lib(contains .lib files)
project build options/
    linker settings/
        C:\GDAL\bld\lib\gdal_i.lib
    search directories/
        [same as in compiler settings]

still I am getting ‘undefined reference to’ errors..

How to solve it? Is there any other way for setting a proper configuration?


texts: Calculate Azimuth from polygon in GeoPandas

texts:
Is there a simple way / ready-to-go function to calculate the azimuth of 2D-polygons in a GeoPandas DataFrame? My goal is to eventually compute azimuth angles of building footprints.
Most of the functions I have found explicitly calculate the azimuth for a line, which is not what I am looking for.
I have been thinking along the following lines to calculate the azimuth of a polygon using its normal vector (2 dimensions). Is my approach going into the right directions? I am a bit unsure how to adjust for the quadrants, also my results seem a bit odd.
azimuth = 90 - math.degrees(math.atan2(normal[1], normal[0]))
if azimuth >= 360.0:
    azimuth -= 360.0
elif azimuth < 0.0:
    azimuth += 360.0


texts:
Looking to find Distance along a line(road) by a user input(Mouse Click). Here is my python code TESTED and works

def onMouseDownMap(self, x, y, button, shift):
 mxd = arcpy.mapping.MapDocument("CURRENT")
 df = arcpy.mapping.ListDataFrames(mxd)[0]
 pointGeom = arcpy.PointGeometry(arcpy.Point(x, y), mxd.activeDataFrame.spatialReference)
 road_layer = arcpy.mapping.ListLayers(mxd, "narr", df)[0]
 arcpy.SelectLayerByLocation_management(road_layer, "WITHIN_A_DISTANCE", pointGeom, "50 FEET")
 arcpy.RefreshActiveView()
 output_layer = arcpy.CopyFeatures_management(road_layer, arcpy.Geometry())
 dist_along_road = output_layer[0].measureOnLine(pointGeom)
 pythonaddins.MessageBox(dist_along_road / 5280, "Address", 1)


texts: New installation of Boundless QGIS failing on start-up - several plugins couldn't load due to an error

texts:
I am fairly new to QGIS, but I had successfully installed QGIS 2.12 on Windows 10 with several plugins, and it worked OK. I wanted to try the Boundless OpenGeo Explorer plugin for QGIS but there was an error on installation. I had had a previous error trying to update the Processing plugin, so I decided to uninstall QGIS, and reinstall a clean version (2.8) from Boundless, as this would have the OpenGeo Explorer plugin pre-installed.

The installation seemed to go OK, but now QGIS will not start properly. It reports something like

  Couldn't load plugin XXX due to an error when calling its classfactory() method

for several plugins (MetaSearch Catalogue Client, qgis2web etc). Then I get

  Couldn't load plugin 'opengeo' from ...",

and the same for plugin processing.

On startups I do get to the QGIS screen, and from the plugins option I notice that there are two entries for the Boundless plugin repository. I cannot delete either of them. And it seems to have carried over the non-core plugins I installed before into what was supposed to be a clean installation. Several of the plugins are said to be invalid.


texts: How to use ST_Difference to create linestrings from non-intersecting polygon boundaries?

texts:
There are two tables with polygon features in my PostgreSQL database. I want to create linestring features where the boundries of the polygons are not intersecting.

CREATE TABLE schema.poly_a
(
gid SERIAL NOT NULL,
geom geometry(polygon, SRID),
CONSTRAINT poly_a_pkey PRIMARY KEY (gid)
);

CREATE TABLE schema.poly_b
(
gid SERIAL NOT NULL,
geom geometry(polygon, SRID),
CONSTRAINT poly_b_pkey PRIMARY KEY (gid)
);

I've tried the following SQL but it's not working. There are no features in the view.

CREATE VIEW line_difference AS
SELECT
row_number() over() AS gid,
g.*
FROM
(SELECT
ST_Difference(ST_Exteriorring(poly_b.geom), poly_a.geom)::geometry(linestring, SRID) AS geom
FROM schema.poly_a, schema.poly_b)
AS g;

A geometry collection instead of linestrings will also do the job.

EDIT #1:

I've tried the SQL query postet by chill. Unfortunately it's not possible to load the layer in QGIS ('There is no entry in geometry_columns!'). DB Manager says there are two linestrings, one multilinestring and one geometrycollection. The preview in the DB Manager is not available.

EDIT #2:

I've already solved it. Just the boundaries for 'ST_Intersects(poly_a.geom, poly_b.geom) IS FALSE' are missing.

CREATE VIEW line_difference AS
SELECT
row_number() over() AS gid,
g.geom
FROM
    (SELECT
    (ST_Dump(COALESCE(ST_Difference(ST_Exteriorring(poly_a.geom), ST_Exteriorring(poly_b.geom))))).geom::geometry(linestring, your_SRID) AS geom
    FROM schema.poly_a, schema.poly_b
    WHERE ST_Intersects(poly_a.geom, poly_b.geom))
AS g;


texts: Adding field to shp subset layer using PyQGIS?

texts:
This is an independent followup from my previous question: QGIS Python generate index field

I have two files (points.shp, polylines.shp), each with a variable amount of sublayers of the given geometry type.

I need to add a simple incremental index to each of these sublayers: Add a field, and fill it with an increasing, unique integer.
Shouldn't be too hard, but I do get more index fields than I desire:

It might be related to the fact that each layer is not an actual individual file, but instead has a source path such as this:

polylines.shp|layerid=0|subset=Layer = 'pathways'
polylines.shp|layerid=0|subset=Layer = 'hedgerows'
points.shp|layerid=0|subset=Layer = 'trees'

For each file:

The first layer to be processed gets a new index field, and the integer values as intended.
The second layer gets the index field with the values correctly, but also gets an index_1 field with NULL values.
The second layer gets the index field correct, and gets index_1 and index_2 with NULL values.
This continues for every added layer, with the last ending up with dozens of new empty fields (index21 in my case) added.

Why does this happen? How can I prevent it?

This is the processing script:

##input_layers=multiple vector

from qgis.core import *
from PyQt4.QtCore import *

def add_index_to_layer(layer):
    fieldindex = layer.fieldNameIndex("index")
    if fieldindex >= 0:
        print "Layer " + layer.name() + " already has an index field"
        return False
    layer.dataProvider().addAttributes([QgsField("index", QVariant.Int)])
    layer.updateFields()
    layer.startEditing()
    fieldindex = layer.fieldNameIndex("index")
    i = 0
    for feature in layer.getFeatures():
        layer.changeAttributeValue(feature.id(), fieldindex, i)
        i +=1
    layer.commitChanges()

input_layers = input_layers.split(';')
for i in input_layers:
    layer = processing.getObject(i)
    add_index_to_layer(layer)


texts: PostGIS Implementation of Brabyn and Skelly Road Network Travel Time Model

texts:
Is there a PostGIS implementation of the Brabyn and Skelly road network travel time model for the New Zealand road network? This model is outlined in:

Brabyn L, Skelly C. Modeling population access to New Zealand public hospitals. Int J Health Geogr. 2002;1(1):3

and

Lauder C, Skelly C, Brabyn L. Developing and validating a road travel  time network for cost path analysis. The 13th Annual Colloquium of the Spatial Information Research Centre. 2001.


texts: Can I also import road type using osm2po?

texts:
I have imported the entire UK OSM dataset using the excellent osm2po. The metatdata of the imported table is shown below and my question is whether any of the columns map to ROAD TYPE.  If not may I instruct osm2po to include road type as one additional column?

Of course, I wish to do this in order to help estimate accurate costs for different types of journey.

Many thanks,

        Table "public.result_2po_4pgr"
    Column     |       Type        | Modifiers
---------------+-------------------+-----------
 id            | integer           | not null
 osm_id        | bigint            |
 osm_name      | character varying |
 osm_source_id | bigint            |
 osm_target_id | bigint            |
 clazz         | integer           |
 flags         | integer           |
 source        | integer           |
 target        | integer           |
 km            | double precision  |
 kmh           | integer           |
 cost          | double precision  |
 reverse_cost  | double precision  |
 x1            | double precision  |
 y1            | double precision  |
 x2            | double precision  |
 y2            | double precision  |
 geom_way      | geometry          |


texts: GeoServer: Show (Random?) Subset of Features

texts:
Background

We have a GeoServer layer backed by a PostGIS table with 22 million smallish features covering the continental US.  We're using GeoServer to serve the features as a WMS layer to a client application that is only interested in the data at large scales (close-up).  To avoid running queries that return millions of features, I've used a style that shuts off display of all features at smaller scales.  This gets the job done, but it's annoying not to be able to get a preview of a larger area (e.g. in GeoServer's OpenLayers layer preview).

Question

Is there some easy way to tell GeoServer, "at small scales, just draw some features"?  I'm not really picky about which features so long as the map generally appears to have data wherever we have data.  A random subset of them would be fine so long as it was "spatially random" (if that makes sense).


texts:
I'm using libgeotiff and proj4 to read and process GeoTIFFs. After reading the file, I want to get the corner coordinates in WGS84/UTM but there are so many functions and parameters to try and understand,

My current function goes like this:

//------------------------------------------------------------------------------
BOOL CGeoTIFFImporter::GetUTMCorner(Corners cnr, double& x, double& y)
//------------------------------------------------------------------------------
{
    BOOL bRet = FALSE;
    if (m_bIsGeoTIFF)
    {
        x = -DBL_MAX;
        y = -DBL_MAX;
        double xtmp = -DBL_MAX;
        double ytmp = -DBL_MAX;
        CString spjUTM = "";
        projPJ pjUTM, pjPJ;
        char *pszProjection;

        switch (cnr)
        {
            case LowerLeft: x = 0.0; y = m_Ypx; break;
            case LowerRight: x = m_Xpx; y = m_Ypx; break;
            case UpperLeft: x = 0.0; y = 0.0; break;
            case UpperRight: x = m_Xpx; y = 0.0; break;
        }

        // Try to transform the coordinate into PCS space
        bRet = GTIFImageToPCS(m_pGTF, &x, &y);

        if (bRet && m_GTFDef.Zone != 0)
        {
            xtmp = x; ytmp = y;
            if (m_GTFDef.Model == ModelTypeGeographic)
            {
            }
            else if (!GTIFProj4ToLatLong(&m_GTFDef, 1, &xtmp, &ytmp))
            {
                return FALSE;
            }

            spjUTM.Format("+proj=utm +zone=%d +ellps=WGS84 +datum=WGS84", UTMZone(xtmp));

            pjUTM = pj_init_plus(spjUTM);

            if (pjUTM == NULL)
            {
                return FALSE;
            }

            pszProjection = GTIFGetProj4Defn(&m_GTFDef);

            if (pszProjection == NULL)
            {
                return FALSE;
            }
            pjPJ = pj_init_plus(pszProjection);
            free(pszProjection);

            if (pjPJ == NULL)
            {
                return FALSE;
            }

            bRet = (pj_transform(pjPJ, pjUTM, 1, 1, &x, &y, NULL) == 0);
        }
    }
    return bRet;
}

Somewhat unsurprisingly given my state of knowledge on the subject, this doesn't work. Can anyone give me some pointers as to what I'm doing wrong?


texts: How to insert 3d line into postgis database?

texts:
I have 3D line which I exported from GRASS in .shp format. When I try to insert this line to Postgis database with QuantumGIS, it gives an error:

  ERROR:  Invalid type name - valid ones are:
      POINT, MULTIPOINT,
      LINESTRING, MULTILINESTRING,
      POLYGON, MULTIPOLYGON,
      CIRCULARSTRING, COMPOUNDCURVE, MULTICURVE,
      CURVEPOLYGON, MULTISURFACE,
      GEOMETRY, GEOMETRYCOLLECTION,
      POINTM, MULTIPOINTM,
      LINESTRINGM, MULTILINESTRINGM,
      POLYGONM, MULTIPOLYGONM,
      CIRCULARSTRINGM, COMPOUNDCURVEM, MULTICURVEM
      CURVEPOLYGONM, MULTISURFACEM,
      or GEOMETRYCOLLECTIONM
  CONTEXT:  SQL statement "SELECT AddGeometryColumn('', $1 , $2 , $3 , $4 , $5 , $6 )"
  PL/pgSQL function "addgeometrycolumn" line 4 at SQL statement

Can someone help me? How can I insert this 3d line into postgis database?


texts: How to get the fill value of a style when using a styleFunction?

texts:
I have a cluster layer which I styled with a style function. I want to be able to have a color picker in my frame for each vector layer so that I could change the fill value of the style independently for each layer.

When I tried to set the initial value for the color picker to the actual fill color of the layer I noticed that if I use layer.getStyle() I don't get a Style object, but a function!

How could I retrieve the fill value from that? Moreover, how could I adjust that value with a color picker?

Edit:
as dvmac01 suggested I am using the forEachFeature function to change the style.

function changeColor() {
    layer.getSource().forEachFeature(function(feature) {
        // this happens when the source is a cluster
        if(feature.get('features')) {
            var resolution = map.getView().getResolution();
            var style = clusterStyle(feature, resolution, this.value);
            feature.setStyle(style);
        } else {
            var style = new ol.style.Style({
                fill: new ol.style.Fill ({color: this.value})
            });
            feature.setStyle(style);
        }
    }, this);
}
colorPicker.addEventListener('input',changeColor);

I try to check whether the current layer is a cluster layer or a regular (polygon) layer. If I change the color of a polygon layer everything works fine. However, for my cluster layer, everything appears to be working fine when I check the styles on the developer console, but for some reason the points are not redrawn with the correct color.

If I type this in the console map.getLayers().item(2).getSource().getFeatures()[0].getStyle().getImage().getFill(); I can see that the value is the one that I selected with the color picker, yet the circles keep the original color.


texts: Importing column from one shapefile to another using PostGIS?

texts:
So I have two shapefiles: one is a hexbinned shapefile with data of how many dogs, cats, and birds in it and another with county lines.

How would I merge the two shapefiles to show which county each hex bin is in? I've tried the following query

create table newhex AS (select * from hexes, counties where ST_Within(hexes.geom, counties.geom));

but I get the error Error executing user query: ERROR:  column "gid" specified more than once

I'm assuming that there are two difference geom columns and that conflicted the geom columns. So I tried another query with explicit selecting the the geom column of the hex bin and the county name with:

select shp_ohhexsm.geom, ohcounties.name from shp_ohhexsm, ohcounties where ST_Within(shp_ohhexsm.geom, ohcounties.geom)

But I get this weird result

Is there any way I can just merge the name of the counties into my hex files without the weird spaces?


texts: Possible to call Rscript from ArcMap Python console without it opening up Windows cmd console?

texts:
I am helping to put together an ArcMap Toolbox, one component of which makes a system call to R's Rscript command-line executable. Everything works nicely, except that while the call to Rscript is being executed, a black Windows command line console opens up. (The console looks much like the one you get if you click the Windows Home button and then type cmd and press return.)

To reproduce what I'm seeing, run the following code in ArcMap's Python Console (I'm using ArcMap 10.5.1.), replacing the path in that first line with the path to Rscript on your own computer.

rscript = "C:/R/R-current/bin/Rscript.exe"
rscript = rscript + ' --vanilla -e'
rcall = rscript + """ "Sys.sleep(10)" """
rcall
## 'C:/R/R-3.1.2/bin/Rscript.exe --vanilla -e "Sys.sleep(10)" '

import subprocess
subprocess.call(rcall)

Is there any way to get this (or some related call to R) to run without ArcMap opening up the distracting Windows cmd console?

(It may be worth mentioning that this behavior seems to be specific to ArcMap's Python console. Running the same Python commands from other Python IDE's (e.g. Visual Studio Code, Spyder) executes the R commands without opening up a Windows cmd console.)


texts: How to append geojson file using featureselection


texts: openlayers3 restore original style of specific feature

texts:
I have a map with LineStrings. There are different types of lines, the vector layer contains a style-class and the styles are dynamically loaded via ajax, i.e. stroke color and width.

My map is interactive, I can select some line strings to assign attributes in my data base. When I select a line string with the mouse, I set the color of the line to red to visulaize all selected line strings on the map.

Now the question: If I un-select a line string, how can I restore the original style for this specific feature (line string)?

A simplified page can be found here: http://mapit.ch/test/map_line_style.html

It uses "Geo Admin API" which is based on openlayers 3. Geo Admin API is used to show different maps from Switzerland, this is removed in the example above. The only function which is used from the API is ga.style.StylesFromLiterals(data) "data" is the json file with the styles, retrieved via ajax (the style file can be changed interactively by the user)

The vector layer is also retrieved via json, each feature has its stlye-class, as a reference to the style file.

I didn't manage to run this in a fiddle, because of the ajax requests and the API which requires a license.


texts: Regrouping concat labels with overlapping intersecting geometries in QGIS

texts:
All the fields i want to show on the map are in the same layer.
All the labels have to be showed on the map.
But the concentration is high.

So i would like to concat the labels when geometries are intersecting. I cannot get intersecting the current geometry with his layer without itself.

I tried to use Overlapping and geomoverlap with rulesbased labeling but nothing works.

The result could be  : resulting label = "label 1 || / || label 2 ||/|| label 3"
Could you show me the way?


texts: Find what polygons are not fully covered by union of polygons from another layer in PostGIS

texts:
I have two tables in a PostGIS database

>>> select postgis_version();
"2.5 USE_GEOS=1 USE_PROJ=1 USE_STATS=1"

One represents the country borders and another one is administrative boundaries. Both are of global coverage covering the whole world. Both tables have spatial indices on the geometry columns.

My goal is to find what countries are not fully covered by the administrative boundaries and, if feasible to achieve within a reasonable time, the location/geometry of part of country geometry where it is not covered by the administrative boundaries.

To give you an example.

Two countries (marked as 1 (consists of three polygons) and 2 (consists of two polygons), dark green polygons (two rows in table in this case).
Three admin boundary polygons (marked as A, B, C) (can be multipolygon/polygon), transparent blue polygons (three rows in table in this case).
Country 1 is fully covered by the admin boundary A and B. Country 2 is not fully covered by admin boundaries (C), so I would like to get this reported. If possible, the hatched area of country 2 could be reported as well.

This sounded like a job for ST_Difference function however my problem is that it is just too slow running on the global dataset.

Running this SQL query on Singapore takes 4 seconds. This is because I have just a few admin boundaries located around the country which I've extracted beforehand and I use only single country geometry. I get the results I expect.

drop table sg_missing_;
WITH SG AS (SELECT iso_cc, geom FROM countries_agg WHERE iso_cc = 'SG'),
    AD AS (SELECT ST_Union(geom) AS geom FROM admins_sample)
SELECT ST_Difference(c.geom, ad.geom) AS geom
INTO sg_missing_
FROM AD ad, SG c

Running the same query on a larger country, such as Indonesia, took 15 hours, so there is something wrong with this approach. I was trying to play smart and get into the analysis only those admin boundaries that intersect the boundary box (envelope) of the country, but it still was running for 10 hours for Indonesia.

--preparing a single country bbox to use
SELECT ST_Envelope(geom) AS geom
INTO country_bbox
FROM countries_agg
where iso_cc = 'ID';

--preparing single country table and union of admin boundary geometries that intersect the given country envelope
WITH
IND AS
(SELECT iso_cc, geom FROM countries_agg WHERE iso_cc = 'ID'),
AD AS
(SELECT ST_Union(ST_Intersection(IND.geom, AD.geom)) AS geom
                                 FROM AD AD, country_bbox IND)

--get geometry of country shape that is not covered by the admin boundaries
SELECT ST_Difference(c.geom, ad.geom) AS geom
INTO sg_missing
FROM AD ad, IND c

It is acceptable if this query would take a few days to run as long as all countries could be analysed in this time. How should such a query be written, if possible in PostGIS? I am open to solutions that involve other tools and software, it does not have to be done solely in PostGIS.

EDIT:

Running

SELECT  <id>,
        ST_Union(geom) AS geom
FROm    (
    SELECT  a.<id>,
            ST_Difference(a.geom, b.geom) AS geom
    FROM    <countries> AS a
    JOIN    <admins> AS b
      ON    a.geom && b.geom
    WHERE   ST_Intersects(a.geom, b.geom)
) q
GROUP BY
        id
;

does not report those country polygons that are located fully outside of admin boundaries, that is, they do not intersect any admin boundaries polygons at all. The bright green polygons should have been included along with the tiny red part of the country polygon that was included due to intersection with the admin boundaries polygon.


texts: I have a clipped surface using Washington State Plane North (feet) added to a QGIS project using the same CRS.  I am trying to calculate the area under this surface (a dirt pile) using SAGA Raster Volume.  The volume returned by SAGA is crazy big.  I have tried different versions of QGIS and get the same results.  The dirt pile has the correct elevation values (in feet above NAVD88).  The surface draws in the correct location.  Why does this tool return such a high number?  Does SAGA Raster Volume not return the volume in map units?


texts: Area of an irregular polygon using side lengths

texts:
I have only lengths for the sides of an irregular polygon, can anyone tell me how I can measure the area of the polygon? Remember only lengths of all the sides, no angles or coordinates.
Few forums mention about trangulation of the polygon etc. But I only have side lengths.
Does anybody have any feedback?


texts: Similar to this rather dated question, is there a way to generate map tiles on QGIS 3.6? The plugins linked from that question only support QGIS 2.x. I wonder if things have changed in QGIS 3.x rendering plugins unnecessary.


texts: Exporting map to TIF with embedded georeference/World using QGIS

texts:
The Map export from the Project menu supports export to TIFF with the World file embedded. The Print Layout supports TIFF export with a separate World file. I need the embedded version.
For automation, what is the API function to replicate: Project -> Import/Export -> Export Map to Image?
I have been unable to find it. There are related answers in the forum, but they deal with separate World files.


texts: Coordinate system for shapefile, no .prj file, using R

texts:
I am trying to use the shapefiles found on this page for oil & gas data from the Netherlands: https://www.nlog.nl/en/files-interactive-map from table at the bottom with "ARC_grid" and "Google Earth (WGS84)" columns and icons.

Specifically the "Oil and gas fields" ARC_grid files (https://www.nlog.nl/sites/default/files/nlog_velden_ed_1950_utm_31n_20170829.zip)

I am using R to download, read and plot the shapefile. I am able to do this OK, but the shapefile does not seem to have a coordinate reference system, I believe due to the absence of the .prj file.

The coordinates I get when I read the file in look like this:

 [1,] 683985.7 5931987
 [2,] 684138.5 5931975

I would like them in Lat/Long, but do not know how to find out what conversion to use initially.

I am looking for an R-based solution.

Using R, this is how I have loaded the file:

temp <- tempfile(fileext = ".zip")
download.file("https://www.nlog.nl/sites/default/files/nlog_velden_ed_1950_utm_31n_20170829.zip", destfile = temp)
filepaths <- unzip(temp)

map <- readOGR(".")

map@proj4string
# CRS arguments: NA


texts:
I am mosaicing rasters in arcpy but it is returning the error, AttributeError: 'module' object has no attribute 'MosaictoNewRaster_management'

I am using the code:

import arcpy
from arcpy.sa import *
import os
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension('Spatial')
shfolders=r'F:\Sheyenne\Final_Imagery\For_Moasic\Sheyenne'
vikfolders=r'F:\Sheyenne\Final_Imagery\For_Moasic\Viking'
out=r'F:\Sheyenne\Final_Imagery\For_Moasic\Mosaiced'
for folder in os.listdir(shfolders):
    if not os.path.isdir(os.path.join(out,folder)):
        os.makedirs(os.path.join(out,folder))
    arcpy.env.workspace=os.path.join(shfolders,folder)
    shrasters=arcpy.ListRasters('*.tif')
    for folder in os.listdir(vikfolders):
        arcpy.env.workspace=os.path.join(vikfolders,folder)
        vikrasters=arcpy.ListRasters('*.tif')
    for sh, vik in zip (shrasters,vikrasters):
            arcpy.MosaictoNewRaster_management([sh,vik],out, sh,"PROJCS['NAD_1983_UTM_Zone_14N',GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-99.0],PARAMETER['Scale_Factor',0.9996],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]","16_BIT_UNSIGNED","30","1","LAST","FIRST")
print "Done Mosaicing"

Within shfolders there are 4 folders with 80 .tif files each and within vikfolders there are 4 folders with 80 .tif files each. I am zipping the .tif files between shfolders and vikfolders and mosaicing the zipped files. The problem here seems to be with arcpy not recognizing the tool. I am using anaconda, I am not sure if that matters.


texts: Can't get layout to be in Map Package (.mpkx) in ArcGIS Pro

texts:
I'm using ArcGIS Pro version 2.7 and I'm trying to export a Map Package (.mpkx) that includes both my map and my layout. I've tried using the 'Map' option in the Package tool bar and tried the geoprocessing tool 'Package Map'. In neither case does the resulting .mpkx contain my layout. I don't see many options for these tools, so I am not sure what I am missing.
Back in ArcMap 10.7 I was able to export a Map Package (.mpk) and it contained both the data view and the layout. That is why I am expecting ArcGIS Pro to include the layout in the package.


texts:
How to use GDAL to export a layered pdf from the map canvas?
I want to export a layered pdf from QGIS. In response to a previous question, I received an answer which implies that it's possible to export a layered pdf from map canvas. Unfortunately, I can't figure out how. I followed the links provided, and read all the information. I still don't understand.
This is the answer:

This exact feature has not be implemented yet for the composer. However, you can export layered pdf from the canvas. You might want to follow the feature request as you will find some tips on the subject : https://issues.qgis.org/issues/9362
See for example :
As GDAL 1.10 supports GeoPDF (See http://www.gdal.org/frmt_pdf.html), one can export a PDF with layers (a raster layer and a polyline vector layer) from Processing inserting a OGR_DATASOURCE parameter in the "Additional creation parameters" box to add the vector layer. However, it only exports vectors without style.
...You can also be a sponsor to have this feature developed ! ;)

This answer has only led me to more questions:

Where is "processing"? Is it the processing menu, or the processing
toolbox?
In the processing toolbox, I see GDAL/OGR [48 geoalgorithms]. None of
the geoalgorithms seem to have anything to do with exporting a pdf. Should I be using one of these geoalgorithms? If so, which one?
Where is the "Additional creation parameter" box?


texts:
I managed to import .dxf into Postgis table, and now have, among others, rows with Point geometries and values in text column (which corespodent with Label values from ACAD .dxf).

What is the best way to display these geometries in Openlayers so they look the same as in any other dxf viewer (format, size, color)? Fonts, sizes and colors are imported into postgis using sql "select *,ogr_style from entities".

The rest of the geometries (lines, polygons) are correctly displayed using HTTP protocol and GeoJSON format, but with labels missing.

Also, is there a way to import Text from .dxf into Postgis as set of lines/polygons? That would also solve the problem.


texts: Can I use the PyQgis QgisPluginLayer to create heatmaps?

texts: I am working in an application to generate maps from a bunch of MIF/MID files, I already created almost all the maps that I need, but, now I want to create heatmaps. I can't find any way of creating them with the current PyQgis API. Recently, I noticed there is something called QgisPluginLayer. I was wondering if there is a way of using those kind of layers to create the heatmaps in PyQgis.


texts: Calculate field result is float with Arc tool, integer with arcpy script

texts:
I am trying to calculate a field for some points within a python script using arcpy functions. (This is calculating the percentage visibility from a cumulative viewshed, if that matters.)

# Process: Get Count
TotalFlightPoints = int(arcpy.GetCount_management(FlightPath_Pts).getOutput(0))

# Process: Add Field
Obsv_Vis = arcpy.AddField_management(Obsv_Vis, "PercentVis", "FLOAT")

# Process: Calculate Field (2)
Obsv_Vis = arcpy.CalculateField_management(Obsv_Vis, "PercentVis", "!RASTERVALU! / " + str(TotalFlightPoints), "PYTHON", "")

Even though I cast that field as a FLOAT when adding it, the result is consistently an integer (0 if less than 1, 1 if the point is 100% visible). I can get around this by "!RASTERVALU! * 100 / " + str(TotalFlightPoints), which will at least give me the percentage (without any trailing decimals). That level of accuracy may be sufficient for this application, at least for now.

But, if I use the Field Calculator within ArcMap on a "float" field (with, say, !RASTERVALU! / 121 as the equation), it gives me the result as a float. So I'm wondering why the field calculator seems to work differently in different situations, and/or what I can tweak to get what I really want.


texts: Field Calculator Python Code for Between X-Y Calculation

texts:
I am trying to identify category based on two continuous variables and using IF ELSE logic for this.

I intent to add elif statement, but even this code is throwing script error.

Both attributes are set to double, but I have also tried this by converting them to long int.


texts: QGIS 3.10: setting default CRS for new projects is ignored

texts:
QGIS 3.10 honors my change of the default CRS for new projects--until I exit... (Caveat: I've discovered that something is changing the default setting even during a session.)
Repeatably, QGIS 3.10, Ubuntu 18.04:

% qgis
Project->Properties->CRS, note default CRS: WGS 84
Settings->Options->CRS, note CRS-for-projects setting: "Use CRS from first layer added."
Change to: Use a default CRS and change default from WGS 84 to EPSG:5070 - NAD83/Conus Albers
still on CRS page, note CRS-for-layers setting: Default CRS: WGS 84 and "Leave as an unknown CRS."
Leave default CRS for layers set to WGS 84
Change to "Use Project CRS"
Unselect "Ask for a datum tranformation" (I don't know which to pick anyway.)
Click "OK" to close Settings window
Project->Properties->CRS, note default CRS still WRS 84
Project->New
Project->Properties->CRS, note CRS has changed to EPSG:5070
Save the new, empty project.
Or don't: the behavior is the same.
Exit QGIS and re-start.
GOTO 1. (that is, changes have been ignored)

Is there a step I'm missing to inform QGIS to make the change permanent?
Thank you.


texts: I want to transform a geojson in a .png image. Do you have any lib/tool/tip to accomplish this task in Java? I googled a lot but no luck.


texts: Creating mask file from shapefile using ArcGIS for Desktop?

texts:
I have a shapefile of Lakes over my study area. Each polygon has attribute information of each lake, name, shape etc. Actually I want to eliminate all lakes from analysis. In other word, I want lakes to have 0 value and the area that out of boundary of lake to have their orginal pixel values. when I do subset using this shapefile, the output is not I want because my shapefile deoesn't have any value outer boundary of lake.

How can I build a proper mask file (inner part of lake has 0 value, outer is 1) using this shapefile in ArcGIS Desktop?


texts: Difference between Ground Sampling Distance and pixel size in RapidEye

texts: I am a bit confused about reading technical specifications on RapidEye satellites. I assumed that Ground sampling distance (GSD) and pixel size express the same meaning, while GSD is 6.5 m, and pixel size is 5 m. Do you know the difference between the two?


texts:
I built an Open Street Maps tile server using the instructions found here:
https://www.linuxbabe.com/ubuntu/openstreetmap-tile-server-ubuntu-18-04-osm
It works fine.
Now I want to add objects such as placemarks to the tiles.  The preferred method would be by doing INSERTs on the underlying database.  (The placemarks will be generated by a program, not added by people.) How does one accomplish this?


texts:
I 'm trying to automate the conversion of raster to polygon and the rasters files are in "TIF" format which is a "unsigned integer". I have the code below. However, I 'm getting an error while executing it. Could you have some suggestion please?

Code:

import arcpy, os
from arcpy import env
from arcpy.sa import *

#To overwrite output
arcpy.env.overwriteOutput = True

#Set environment settings
env.workspace = "C:/Subhasis/Test/Neshanic_Python/extract"

outws="C:/Subhasis/Test/Neshanic_Python/extract/clip/"

#checkout ArcGIS spatial analyst extension license
arcpy.CheckOutExtension("Spatial")

inraster = arcpy.ListRasters("*", "TIF")

for i in inraster:
    filename=os.path.splitext(i)[0]
    outPolygons= outws + str(filename) + '.shape'
    field = "VALUE"
    arcpy.RasterToPolygon_conversion(inraster, outPolygons, "NO_SIMPLIFY", field)

  Error: Traceback (most recent call last):   File
  "C:\Subhasis\Test\Neshanic_Python\extract\Clip.py", line 23, in

      arcpy.RasterToPolygon_conversion(inraster, outPolygons, "NO_SIMPLIFY", field)   File "C:\Program Files
  (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\conversion.py", line 241, in
  RasterToPolygon raise e


texts: Points drawn outside a certain area get placed on perimeter of smaller area

texts:
I am trying to add points to a table in MapInfo (geologic map) but cannot get points to display in correct location. Are there boundaries for tables?
How to change this?

I'm using Mapinfo 10.5.


texts: Transparency parameter of colors (vector, legend,...) not working in QGIS?

texts:
It's a pity that QGIS offers to set transparency in the selecting color tool, but doesn't save it...

Explanation :
I guess that the AA (in French, sorry) in #RRVVBBAA would enable to set transparency (0% => AA = 00 ; 50% => AA = 80 ; 100% => AA = FF).
The same thing with rgba (r, g, b, a) where QGIS doesn't allow to set a <> 1.

For printing issues, I would like to have a black 50%, it would be possible by setting #00000080 or rgba( 0, 0, 0, 0.50 ) but QGIS re-set theses values back to #000000ff and rgba( 0, 0, 0, 1.00 )...

If there is no way to do this thing now, do you think it will be improved in a next version ?


texts: Finding out corner shape in vector polygon file in Qgis

texts: how can i findout and distinguish a corner  shape in Qgis as shown in image?


texts: proj: Oblique Mercator specifying two points on central line

texts:
I'm trying to use Oblique Mercator in QGIS 3.4, which uses proj to calculate coordinates conversion. I specify two points in the central line, rather than specifying the center and alpha.

Oblique Mercator: Usage explains

  In the second case, the azimuth is given indirectly by specifying two points on the central line, using the options +lat_1, +lon_1, +lat_2, and +lon_2.

Here are the two locations I'm using:

Tokyo: latitude: 35.652832, longitude: 139.839478
New York: latitude: 40.730610, longitude: -73.935242

Calling proj command (Rel. 6.1.0) specifying the two locations for the central line of Oblique Mercator and converting the two points to x-y coordinates gives:

$ echo -73.935242 40.730610 |proj +proj=omerc +lat_1=35.652832 +lon_1=139.839478 +lat_2=40.730610 +lon_2=-73.935242 +ellps=GRS80
5162164.32  14191389.18
$ echo 139.839478 35.652832 |proj +proj=omerc +lat_1=35.652832 +lon_1=139.839478 +lat_2=40.730610 +lon_2=-73.935242 +ellps=GRS80
1450567.63  3987778.87

In the output, I expect the "y" parts of the two coordinates are the same, because these two points are on the central line of projected coordinates. However, they are not (14191389.18 != 3987778.87).

Can somebody explain how omerc with two points in the central line works and provide correct parameter for proj?


texts:
I get this error when I run my script:

  Traceback (most recent call last):
  File "", line 1, in  AttributeError: 'Symbology' object has no attribute 'renderer'

Script:

p = arcpy.mp.ArcGISProject("CURRENT")
m = p.listMaps('Map')[0]
l = m.listLayers('Raster1')[0]
sym = l.symbology
sym.renderer

https://pro.arcgis.com/en/pro-app/arcpy/mapping/symbology-class.htm

According to this URL's page there shouldn't be any error.


texts: Basemap function on webmap returns AttributeError: 'Item' object has no attribute 'basemap'

texts:
I am working with an ArcGIS Notebook and have the following script (see attached image as well):

from IPython.display import display
import arcgis
from arcgis.gis import GIS
from arcgis.mapping import WebMap

gis = GIS("home")

webmap = gis.content.search(query="title:empty_webmap*", item_type="Web map")
webmap

empty_webmap1= webmap [2]
empty_webmap1

empty_webmap1.basemap

I then receive the following error message:
KeyError
Traceback (most recent call last)
/opt/conda/lib/python3.6/site-packages/arcgis/gis/__init__.py in __getattr__(self, name)
   7785         try:
-> 7786             return dict.__getitem__(self, name)
   7787         except:

KeyError: 'basemap'

During handling of the above exception, another exception occurred:
AttributeError
Traceback (most recent call last)
<ipython-input-35-5e8ee0740a1e> in <module>
----> 1 empty_webmap1.basemap

/opt/conda/lib/python3.6/site-packages/arcgis/gis/__init__.py in __getattr__(self, name)
   7786             return dict.__getitem__(self, name)
   7787         except:
-> 7788             raise AttributeError("'%s' object has no attribute '%s'" % (type(self).__name__, name))
   7789
   7790     def __getitem__(self, k): # support item attributes as dictionary keys on this object

AttributeError: 'Item' object has no attribute 'basemap'

Could this be due to my environment e.g. an older Python API? And if yes, how could I check this? I am working in an virtual environment.


texts: How to add geometry to current section in arcpy?

texts:
I looping through feature class. I want select the features that are 3 meters distance but not intersected features. I know we can do by select layer by location or attribute.

with arcpy.da.UpdateCursor("ROAD_CENTERLINE",["shape@","S_CODE","CATEGORY","SUBCAT","oid@"]) as usc:
    for row in usc:
        geom,o=row[0],row[4]
        print o
        sbl=arcpy.SelectLayerByLocation_management("ROAD_CENTERLINE","INTERSECT",geom,"3 Meters","NEW_SELECTION","NOT_INVERT")
        sbl=arcpy.SelectLayerByLocation_management("ROAD_CENTERLINE","INTERSECT",geom,"","REMOVE_FROM_SELECTION","NOT_INVERT")
        sbl=arcpy.SelectLayerByLocation_management("ROAD_CENTERLINE","ARE_IDENTICAL_TO",geom,"","ADD_TO_SELECTION","NOT_INVERT")

to do above operation any simple step is there?

input

required output selection is


texts: Unable to pull a bts.dot.gov MapLayer using Leaflet in R

texts:
I'm trying to get the noise pollution layer from the following bts.dot.gov ESRI MapServer to render in R:

https://maps.bts.dot.gov/services/rest/services/Noise/CONUS_road_and_aviation_noise/MapServer

I'm using the following code in RStudio (all necessary packages are installed):

leaflet() %>%
  setView(-117.1, 32.9, 10) %>%
  addEsriBasemapLayer(esriBasemapLayers$Oceans, autoLabels = TRUE) %>%
  addEsriDynamicMapLayer(
    url = paste0("https://maps.bts.dot.gov/services/rest/services/Noise/CONUS_road_and_aviation_noise/MapServer"))

However, none of the noise pollution data is being rendered; I only see a blank map.

The following, virtually identical code for an ArcGIS MapServer does work (wildfire data):

leaflet() %>%
  setView(-117.1, 32.9, 10) %>%
  addEsriBasemapLayer(esriBasemapLayers$Oceans, autoLabels = TRUE) %>%
  addEsriDynamicMapLayer(
    url = paste0("https://maps7.arcgisonline.com/arcgis/rest/services/USDA_USFS_2014_Wildfire_Hazard_Potential/MapServer"))

Am I missing something here? The two map services seem very similar, but I have little experience. The BTS map does support dynamic layers, while the Wildfire map does not, but I have tried pulling the map using this option (among others) as well:

options = dynamicMapLayerOptions(dynamicLayers = "0")


texts: How to set up QGIS print composer options from the python console?

texts:
If I have a print composer called "ABC", I would like to "tell" its composer  that I would like to

update "Map 0" to simulate the "set to map canvas extent" button activation
update certain text labels
activate "export to image"

Thank you.


texts: How do I use Python to edit a field in an attribute table?

texts:
I have a few data tables full of information on salt water injection wells, with injection volumes given on a monthly basis. I need to add these monthly values up in order to get an annual value for each well to use on my map. Due to the sheer number of wells (about 1100 per data table), I would prefer to do this using a Python code rather than by hand.

How would I go about populating the new AnnualVolumes field with values calculated from another field?

I'm using ArcMap 10.3.


texts: Downloading new WorldClim version 2.1 bioclimatic variables for specific region of world using R

texts:
I want to download the new (19) WorldClim Bioclimatic variables (version 2.1) for a specific region of the world, but it looks like you can only download the entire world, which is around 10 GB!
Around 8hrs to download with my internet connection.
I know there is a way to download tiles in R using the raster library. Here is the code:
library(raster)
getData('worldclim', var='bio', res=0.5, lon=-60, lat= -37, download = TRUE, path = "")

    *trying URL 'https://biogeo.ucdavis.edu/data/climate/worldclim/1_4/tiles/cur/bio_44.zip'
    Content type 'application/zip' length 5663180 bytes (5.4 MB)
    downloaded 5.4 MB*

As you can see, this works for worldclim version 1.4 (data are obsolete). Anyway to do the same but with worldclim 2.1?


texts:
I have a base 30 m.Dem . However, i also have 1m Dem for an area of interest only. I would like to overlay the 1m Dem on the 30 m Dem meaning the area of interest will have 1m and the rest 30m.

Is that possible?


texts: Change font size of all labels print composer

texts: I have to print a map in A0 paper size and also in A3. First I completed the job for the A0 size map, then I duplicated the print composer and I am changing the fonts of attribute tables and label boxes, one by one, to a smaller size for the A3 map. Is there a way to change all the fonts for every text in the print composer, at once? (QGIS 2.18)


texts: tabbed popup window for Arcgis Javascript

texts: Anyone know examples of a tabbed popup window for Arcgis for javascript, NOT an infowindow? If you could point me in the the right direction it would help.


texts: Extract OSM nodes: all surface parking spots in germany

texts:
I would like to extract all surface parking spots in germany as nodes from osm.
Gladly I have found overpass-turbo.eu, but unfortunately it is not working for my purpose...
Either the server is too busy or I get a timeout.
My request is:

/*
This has been generated by the overpass-turbo wizard.
The original search was:
“parking=surface in germany”
*/
[out:xml][timeout:500000];
// fetch area “germany” to search in
{{geocodeArea:germany}}->.searchArea;
// gather results
(
  // query part for: “parking=surface”
  node["parking"="surface"](area.searchArea);
);
// print results
out body;
>;
out skel qt;

Am I doing something wrong?


texts: Calculating Bearing/Direction via Spatial Join?

texts:
I was just wondering if it was possible to calculate a bearing field via a spatial join?

Basically I have needed in the past to calculate distances from one singular point (e.g a well) to a number of different closest attributes for different shapefiles (such as environmental sensitivities, so SPAs, SACs, Pipelines etc, therefore ranging all different vector formats.)

I found out you could do this simply by using a spatial join so I created a model using model builder that basically created spatial joins of my input point to every single environmental variable (selecting the CLOSEST parameter), and after merging all the ouput feature classes gives me an output feature classes with the distances of every closest sensitivity.

However, I also need the bearing from the input point to each environmental variable, however I cannot figure out if this is even possible using the current method I'm using, even though it seems like it should be an easy enough variable to include...

I have tried using "Add Field" and "Calculate Field" after the spatial joining but this does not seem to work, so I wonder if there is somehow a way to do it during the spatial joins?

Or any other way?


texts: Whats the unit of circle's radius in Google Maps

texts:
I need to know what is the unit of radius (e.g. miles, kilometers etc) in following question:
How to place markers on the outline of a circle in google maps?


texts: Increasing Contrast of Topographic Image in R

texts:
Brand new to R, so a little foreign to me, but I Created a topographic image using Landsat, and .hgt data. I got an image out but tried to increase contrast using the scales tool. It is still giving me out a whitewashed image, and if I try to increase it form c(0,1) it gives me an error. {See both below}

What should I do to make this look better?

install.packages(c("rayshader", "raster", "sp"))
library(rayshader)
library(sp)
library(raster)
library(scales)
---
memory.limit(size = 56000)
---
elevation1 = raster::raster("C:/Users/jacob/OneDrive/Desktop/N47W123.hgt")
elevation2 = raster::raster("C:/Users/jacob/OneDrive/Desktop/N47W124.hgt")

wash_elevation = raster::merge(elevation1,elevation2)

height_shade(raster_to_matrix(wash_elevation)) %>%
  plot_map()

------

wash_r = raster::raster("C:/Users/jacob/OneDrive/Desktop/LC08_L1TP_047027_20140915_20170303_01_T1_B4.TIF")
wash_g = raster::raster("C:/Users/jacob/OneDrive/Desktop/LC08_L1TP_047027_20140915_20170303_01_T1_B3.TIF")
wash_b = raster::raster("C:/Users/jacob/OneDrive/Desktop/LC08_L1TP_047027_20140915_20170303_01_T1_B2.TIF")

wash_rgb = raster::stack(wash_r, wash_g, wash_b)
raster::plotRGB(wash_rgb, scale=255^2)
---
wash_rgb_corrected = sqrt(raster::stack(wash_r, wash_g, wash_b))
raster::plotRGB(wash_rgb_corrected)
--
raster::crs(wash_r)
## CRS arguments:
##  +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
---
raster::crs(wash_elevation)
## CRS arguments:
##  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
---
crs(wash_r)
## CRS arguments:
##  +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
---
wash_elevation_utm = raster::projectRaster(wash_elevation, crs = crs(wash_r), method = "bilinear")
crs(wash_elevation_utm)
## CRS arguments:
##  +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
---
bottom_left = c(y=-123.704087, x=47.013711)
top_right   = c(y=-122.124594, x=47.930524)

extent_latlong = sp::SpatialPoints(rbind(bottom_left, top_right), proj4string=sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
extent_utm = sp::spTransform(extent_latlong, raster::crs(wash_elevation_utm))

e = raster::extent(extent_utm)
e

---
wash_rgb_cropped = raster::crop(wash_rgb_corrected, e)
elevation_cropped = raster::crop(wash_elevation_utm, e)

names(wash_rgb_cropped) = c("r","g","b")

wash_r_cropped = rayshader::raster_to_matrix(wash_rgb_cropped$r)
wash_g_cropped = rayshader::raster_to_matrix(wash_rgb_cropped$g)
wash_b_cropped = rayshader::raster_to_matrix(wash_rgb_cropped$b)

washel_matrix = rayshader::raster_to_matrix(elevation_cropped)

wash_rgb_array = array(0,dim=c(nrow(wash_r_cropped),ncol(wash_r_cropped),3))

wash_rgb_array[,,1] = wash_r_cropped/255 #Red layer
wash_rgb_array[,,2] = wash_g_cropped/255 #Blue layer
wash_rgb_array[,,3] = wash_b_cropped/255 #Green layer

wash_rgb_array = aperm(wash_rgb_array, c(2,1,3))

plot_map(wash_rgb_array)
---

wash_rgb_contrast = scales::rescale(wash_rgb_array,to=c(0,1))

plot_map(wash_rgb_contrast)
---
plot_3d(wash_rgb_contrast, washel_matrix, windowsize = c(1100,900), zscale = 15, shadowdepth = -50,
        zoom=0.5, phi=45,theta=-45,fov=70, background = "#F2E1D0", shadowcolor = "#523E2B")
render_snapshot(title_text = "West Washington | Imagery: Landsat 8 | DEM: 30m SRTM",
                title_bar_color = "#1f5214", title_color = "white", title_bar_alpha `


texts: Button that allows project search and zooms to site once executed

texts:
I am trying to create a button that will execute a Python script. I have limited Python experience. What we are trying to accomplish is to recreate these steps:

Query on project Site and Radius Layers (both in an SDE)
Zoom to the site
Adjust scale to 10,500

I have pasted the code I wrote while trying to figure this all out, see below:

import arcpy
mxd = arcpy.mapping.MapDocument('CURRENT')
df = arcpy.mapping.ListDataFrames(mxd, "Layers") [0]
df.zoomToSelectedFeatures()
df.scale = 10500
arcpy.RefreshActiveView()

This works, but is not anything that I can't do with ArcMap tools.