texts: Is there a way to delete a raster using arcpy / ArcGIS 10.0?

texts:
I'm trying to find a way to delete a raster dataset (in any format) using arcpy. As it is now I have been writing custom functions for each file type (eg, delete directory and pyramids for grids, delete .tif, .tfw, ...etc for geotiffs).

Seem like this should be an obvious thing but I can't find it in the function list.


texts: QGIS, 3D-view, Time Manager Plugin. How to update 3D-view automatically?

texts:
I have a set of polygons with values, and time information.
I am able to set Time Manager plugin to work with the values also with 2.5D view. So every year with new value for polygon I have a higher or lower polygon respectively according to values.
When I try to start it in 3D view, it also starts, but when Time Manager is playing with values via filter operator, it seems 3D view is not refreshing itself. If I go to properties panel of a layer, and I click apply, it is applying new value. But it is not refreshing automatically.

Is there any way to fix it?


texts: Maximum Likelihood Classification says there are 0 classes when there should be 5

texts:
I'm trying to do an image classification. I created a polygon with 5 different classes (identified in the Id column; 1-5). I then created a signature file (which took about an hour to complete). When I did the maximum likelihood classification using that sig file, it gave me an immediate error saying there are fewer than 1 class.

I set the data frame to Canada Lambert Conformal Conic. I added the scanned paper map. Added a townships shapefile. Georeferenced it using about 6 control points. Created new shapefile in ArcCatalog using Canada lambert conformal conic as the CRS. Created 3 polygons for each class, differentiating them using 1,2,3,4 or 5 in Id column. Created signature in same folder as other data. Performed MLC using that sig file but it immediately gave me that error.

How do I fix this?


texts: QGIS keyboard shortcut to launch a 'Select by expression' query

texts:
I can't find if it is possible to launch a 'Select by expression' query (not to open the dialog) with a keyboard shortcut. I can't find this option in the Settings > Configure shortcuts list...

If you have an idea on this question...


texts: Assigning row-wise value to rotated grid of points using ArcGIS Desktop?

texts:
I am working with a climate model that operates by means of a rotated grid of points. The rotation angle is about 15 degrees counterclockwise. The points are provided with (N, E) coordinates, but their IDs are nonconsecutive, and have irregular intervals.

In the image below, the bottom left point has OBJECTID=35, the one to its right has OBJECTID=40, the third to its right has OBJECTID=50, but then the leftmost point in the row above has OBJECTID=51.

My problem: I want to assign each point a new ID (integer) that starts from 0 (bottom left point) and increases towards the right hand side of the row, and row by row until it reaches the upper right point (highest ID).

The rotation causes problems because none of the coordinates are steady, but is there any way I could do this?

The grid has 400 points in the blue arrow direction, and 300 in the green arrow direction.

ArcGIS version: 10.3.1


texts: Is it possible to have the output features of a model have a specific visual style? The features come out in a random assortment of colors and very small. Is there a tool for this? or any suggestions?


texts: How to set Tomcat to completely prevent browser from caching my PBF vector tile?

texts:
I need to publish a geometry table (in PostGIS) as a layer of vector tiles (in PBF format) in GeoServer. As the table is dynamically written by a production web app, i need to set the layer's caching as disable in both server and client.
For that reason, in the layer's Tile Caching in Geoserver, I set -1 in both:

Expire server cache after n seconds
Expire client cache after n seconds

like the following picture:

However, my observation shows that the above only prevents server cache and not client cache, as newly added points will only be shown in my OpenLayers page when any of the following applies:

when the browser is cleared from cache
when the tile is newly received from server (not from cache)

So i dig in more references especially on how to set Tomcat to prevent client from caching. Based on a reference I tried to add ExpiresFilter in the Tomcat's web.xml:

<filter-name>ExpiresFilter</filter-name>
<filter-class>org.apache.catalina.filters.ExpiresFilter</filter-class>
<init-param>
  <param-name>ExpiresByType application/x-protobuf</param-name>
  <param-value>access plus 0 seconds</param-value>
</init-param>

The result is now each tile has new 4 headers in Response Headers: max-age=0, Expires date (with expected value), and 2 more headers; like the following image.

Now the problem: page refresh in browser gets these tiles from browser cache, with HTTP 304 status, and the new 4 headers are now gone; like the following image.

In this situation, any newly added points in my PostGIS table will only appear in OpenLayers when the corresponding tile never gets to browser before.

Question:

Is this normal?
How to completely prevent browser from caching my PBF vector tile?

My ecosystem:

OpenLayers 5.1.3
GeoServer 2.13.2 on Tomcat 8.5.32
Windows 10 Pro 64-bit
PostgreSQL 9.6.1 with PostGIS 2.3
browser: FireFox, IE


texts:
My goal is to convert gps coordinates to the format which the OLEX plotter has, so I can import them as marks.

So I created a point at 64 10.239N and 51 41.523W
I exported the point, and the file says 3850.2395026 -3101.5235411 1495562630
The last one seems to be the time, and is not important to me at this point.

So my question is: what format does the OLEX plotter have?
And possibly how to calculate this?


texts: Is it possible to fade the ends of  a line? I have a map the shows the city of limits and because of its shape, there are some areas outside of the city limits where we don't have major road data in  but we have centerline data. So I was wondering if there was a way to get the lines to fade out instead of abruptly ending?


texts:
My image is gray scale image which i have converted to binary by smoothing and then edge detection method. Now i want to place the circle /ellipse on edge detected. I don't know how to proceed further. I am using visual C++ and gdal libraries for image processing and Qt for GUI. I should get the output as shown in reference image i have attached.

I tried the following code to detect circle. It misses some circle and the size of the circle is same in the image. Please help me out . I am not understanding the argument (specially the threshold). Is canny edge detection is performed internally in houghcircle?

using namespace cv;

int main(int, char** argv)
{
   Mat src, src_gray;

   /// Read the image
   src = imread("D:\\hardik\\crater8.tif");
          if( src.empty() )
     { return -1; }

   /// Convert it to gray
    cvtColor( src, src_gray, CV_BGR2GRAY );

   /// Reduce the noise so we avoid false circle detection
    GaussianBlur( src_gray, src_gray, Size(9, 9), 2, 2 );

    vector<Vec3f> circles;

   /// Apply the Hough Transform to find the circles
    HoughCircles( src_gray, circles, CV_HOUGH_GRADIENT,2.0,100.0,128, 64,1, 60 );

   /// Draw the circles detected
    for( size_t i = 0; i < circles.size(); i++ )
    {
         Point center(cvRound(circles[i][0]), cvRound(circles[i][1]));
         int radius = cvRound(circles[i][2]);
         // circle center
         circle( src, center,3, Scalar(0,255,0), -1, 8, 0 );
         // circle outline
         circle( src, center, 30, Scalar(0,255,255),2, 8, 0 );
    }

   /// Show your results
    namedWindow( "Hough Circle Transform Demo", CV_WINDOW_AUTOSIZE );
    imshow( "Hough Circle Transform Demo", src );
//  imwrite("C:\\Test2\\val1.jpg",src);
    waitKey(0);

    src.release();
    src_gray.release();

    return 0;
}


texts: Creating centreline for every polygon in submarine canyon database using ArcGIS for Desktop Basic?

texts:
I have the Basic license for ArcGIS and so I cannot use the collapse dual lines to centreline tool or the river bathymetry toolset by ESSA as they both need an advanced license.

Any idea how I can find the centreline or thalweg of the canyons using another tool which works in the standard version of ArcGIS please?


texts: Programming a geoweb application for cloud. What does it imply?

texts:
I've recently started to develop a geoweb application using GeoDjango + PostGIS in the backend, HTML5 + JS in the frontend. I'm really at the beginning, so, I have time to change a lot of things.

Basically, I want to store raster/vector data, allow some basic operations over them and visualize the result. Simplified, translate the concept of a Desktop GIS (QGIS, OpenJUMP...) to the web

My big issue right now is I'd like to operate over a cloud infrastructure, but I'm not sure of the implications from a developer point of view. I'm learning about the concept "Dev-Ops", and I've read some articles, like this, but I'm afraid of doing things wrong, because I never built an application specially conceived to be deployed in a Cloud infrastructure before. And I don't want to incur in an "analysis paralysis" situation.

For example, should I first care about installing PostGIS in a Cloud infrastructure and test how to scale it?

I think my problem is I don't know what kind of things are pure "developer things" and what kind of things are "enviroment configuration". And how to mix them. What comes first, what comes next, etc.

Any advice is welcome. Many thanks in advance.


texts: Calculating mean between the intersection of two polygons

texts: I need to calculate the mean value (for a specific field of the attribute table) at the intersection of many polygons (.shp). I converted the polygons to raster based on the value I want to use to calculate the mean. I´m using Arcmap10.2. I’m open to suggestions.


texts:
I am unable to use pdb from qgis console.
If I run a function that calls pdb.set_trace() I get the error

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "D:\Scripts\routDesigner.py", line 29, in ttt
    pdb.set_trace()
  File "C:\PROGRA~1\QGISBR~1\apps\Python27\lib\bdb.py", line 53, in trace_dispatch
    return self.dispatch_return(frame, arg)
  File "C:\PROGRA~1\QGISBR~1\apps\Python27\lib\bdb.py", line 91, in dispatch_return
    if self.quitting: raise BdbQuit
BdbQuit

prefixing pdb.set_trace() with  QtCore.pyqtRemoveInputHook() does not change anything. Finally, using pdb.pm() has also strange effects, it seems that it prints (Pdb) and quits the debugging environment before I can write anything.
I'm using qgis 2.6 ob windows 8


texts: Why ArcGIS Topology "Must Not Have Dangle" can be fixed with Planarize tool (Advance Editing) even its still have dangle?

texts: I have some project fixing line features that "must have no dangle" but when I use planarize tool on error part, the error was fixed. can anyone explain this?


texts:
In the code below the user picks from one of the following "GetParameterAsText" values:

Tampa, Orlando, Denver, Cosprings

# Variables/Parameters
market = arcpy.GetParameterAsText(0)
nfid_fqn_associations = "G:\\Project_Verizon_" + **market** + "\\02GISData\\VZ3GIS\\SHP Workspace"
workspace = "G:\\Project_Verizon_" + **market** + "\\02GISData\\VZ3GIS\\File Workspace"

# Create table
arcpy.TableToTable_conversion(nfid_fqn_associations, workspace, "SHP_Report", "Where Clause")

How do I make it so if they select Tampa as their input variable it inserts the value "57" into the arcpy.TableToTable_conversion() "Where Clause"?

Associated integer:

Tampa : 57

Orlando : 1

Denver : 45

Here's how I'd attempt:

# Variables/Parameters
market = arcpy.GetParameterAsText(0)
value = 0
if market == "Tampa":
    value = 57
if market == "Denver":
    value = 45

# Create table
arcpy.TableToTable_conversion(nfid_fqn_associations, workspace, "SHP_Report", "MARKETID = " + "'" + value + "'")

This script will be turned into a geoprocessing tool so I'm not entirely sure how to take the users input and initialize the variable "value" based off the associated integer above.


texts: How to move layers in the Layer Order Panel using PyQGIS?

texts:
Related to the following question but looking for a PyQGIS method:

QGIS Layer Order Panel - add layers at top of layer order

The following is a simple setup containing a group with three layers.

The code I use adds a new layer at the end of this group:

root = QgsProject.instance().layerTreeRoot()
group = root.findGroup('Main group')

vlayer = QgsVectorLayer('LineString?crs=epsg:27700', 'vlayer', 'memory')
QgsMapLayerRegistry.instance().addMapLayer(vlayer, False)
group.insertChildNode(-1, QgsLayerTreeLayer(vlayer))

In the Layer Order Panel, the newly added layer is at the end.

Is it possible to move this to the top without moving the layer in the Layers Panel?


texts: Enabling preview pane of Windows Explorer for mxd files?

texts:
The preview pane in Windows explorer is an efficient way to have sneakpeak at files without having to open them. It works great for a number of file extensions; images, office documents, pdf, dwg files etc.

However, it does not work for mxd files. I usually work around this by having recent exports from each mxd close to it, but now I'm being forced to keep exports in a separate folder.

Is there a way to enable/activate previews for mxd files? I have found this page for example, but the method suggested seems only to work if I have an image/video or plain text.


texts: Add a timestamp to a dataframe in Python/Earth Engine

texts:
I came across an example kindly shared on youtube by Madhu Mysore. Here is the link: https://www.youtube.com/watch?v=IGUYSeoeHhg&t=2105s.

I changed the code as I am not interested in regions at the moment, but rather in daily/weekly observations. What I am struggling to do is to add a timestamp to the following code. I checked multiple forums and EE developers page but cannot figure it out.

import ee
import pandas as pd
ee.Initialize()
districts = ee.FeatureCollection('users/mdm/za/za_districts')

chirps = ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY").select('precipitation').filterDate('2018-01-01', '2018-12-31').mean()
evi = ee.ImageCollection('MODIS/006/MOD13A2').select('EVI').filterDate('2018-01-01', '2018-12-31').mean()

composite = ee.Image.cat([chirps, evi])

dist_stats = composite.reduceRegions(districts, 'mean', 5000)

dist_stats = dist_stats.select(['precipitation', 'EVI'], ['tot_rainfall', 'med_evi'], retainGeometry=False).getInfo()

dtstats_df = pd.DataFrame()
for dist in dist_stats['features']:
    df = pd.DataFrame([dist['properties']],columns=dist['properties'].keys())
    dtstats_df = pd.concat([dtstats_df, df], axis=0)

dtstats_df = dtstats_df[['tot_rainfall', 'med_evi']]

dtstats_df


texts: How to query data from an Open Data Cube (Python / Win 10)?

texts: I would like to learn how to access and manipulate satellite images from an Open Data Cube (ODC) using Python. I have tried to follow the instructions found in https://datacube-core.readthedocs.io/en/latest/user/intro.html. The first step (installation) went smoothly, but I got stuck in the second step (database setup). After downloading, installing, and executing the suggested PostgreSQL version, I wasn't able to type the necessary commands described in https://datacube-core.readthedocs.io/en/latest/ops/db_setup.html. Even copy-pasting the commands did not work. I am not sure of the reason behind this, and I would like to know if there is a simpler way to access an ODC and query data from it online (without setting up a local database and configuring an environment).


texts:
I have a polygon layer and want to change the colour of single polygons.

Usually I use a rule based symbolization and filter the polygons be their "ID". But by that I have 30 or more rules.

Is there any other way to change the colour of single polygons in the same layer? For example with select features by area or single click ?


texts: Problems using r.reclass in QGIS: invalid reclass rule error?

texts:
I need to reclass some rasters, but i am having a couple of problems, i can´t generate a new raster with the reclassification using the r.reclass of GRASS.
For the question here i´m using as example (images attached below) a slope raster that i need to reclass for exact values (declives_rules). When i run the operation, happens a error!

I wonder if i could use the raster calculator to reclass this raster?


texts: QGIS: raster image georeference and overlay

texts:
I've created a grid and points in QGIS which I used to survey. Using a test image at the moment and I'm trying to fit this image as an overlay within a boundary. When I insert the image it has a white background which I can't remove in QGIS.

Image 1 is what the area looks like, I want to get the survey results inside the red line. Image 2 is what happens when I insert the image, the inserted image of the geophysical results is a .BMP file, the same happens with .TIFF files.

The image is moved into place by georeferencing corners of the boundary area or points on the grid.

Image 1:

Image 2:


texts: How upload KML from your local disk using Google Maps API?

texts:
I am trying to get load my KML from my local disk with given below code :

<script type="text/javascript">
     function initialize() {
    var mapOptions = {
      center: new google.maps.LatLng(43.321675,-75.172631),
      zoom: 8,
      mapTypeId: google.maps.MapTypeId.ROADMAP
    };
    var map = new google.maps.Map(document.getElementById("map_canvas"),mapOptions);
    var ctaLayer = new google.maps.KmlLayer('C:/Users/Ashutosh_k/Downloads/newyork.kml');
    ctaLayer.setMap(map);
  }
</script>

Problem : When my google map is displaying its not showing up the elements of KML .


texts: Rmoving single class from ArcGIS legend for graduated polygons

texts:
I have the following legend item:

I would like to remove the class in the red circle from the legend. In QGIS you can simply delete or remove it.

How can one do this in ArcGIS without affecting the actual graduated symbology?


texts: Rotating features in memory layer using PyQGIS

texts:
I can rotate polygon features in a shapefile using:
vlayer = iface.activeLayer()
vlayer.startEditing()
for ft in vlayer.getFeatures():
    geom = ft.geometry()
    pt = ft.geometry().centroid().asPoint()
    geom.rotate(45, pt)
    vlayer.dataProvider().changeGeometryValues({ 0 : geom })

vlayer.commitChanges()

But if I use this on a memory layer, nothing happens. True values are printed on the console but the features are not rotated, even after refreshing the canvas and repainting the layer.
Is there a way to do this for memory layers?


texts:
How could I get the lat,long of an object by knowing its xy values(in meters) from my vehicle (assuming it is(0,0) the center)?
For example my vehicle lat,long is 30,31 but it's xy values are always going to be (0,0). The object xy for example is (1.3,1.2) from the vehicle (in meters). How could I then change it to a lat and long?


texts:
I'm dealing with a varierty of spatial polygon files (shapefiles, file/personal gdb feature classes, sde feature classes (some with spatial views)) and depending on the file type, I get a mulittude of Shape Area feature class field names (the required field, not a created one)

So for example, the field names vary from "Shape_Area", "SHAPE_Area", "SHAPE.AREA", "GEOMETRY_Area", "GEOMETRY.AREA"

I'm writing a python script to pull out the area of all these types of feature classes and I do something along the lines of:

for field in arcpy.ListFields(fc):
    if field.name in ["Shape_Area", "SHAPE_Area", "SHAPE.AREA", "GEOMETRY_Area", "GEOMETRY.AREA"]:
         do something

I'm wondering if there's a more clever way of pulling out the shape area field, instead of just listing out all the ways it can be named? There doesn't seem to be a field property for something like Area. is there any other way of doing this?


texts: Calculating distance (lat, lon) for geom point using PostGIS

texts:
I have point in my db (Postgres postgis):

CREATE TABLE geometries (name varchar, geom geometry);

With example data:

INSERT INTO geometries VALUES ('Point', 'POINT(0 0)');

INSERT INTO geometries VALUES ('Point', 'POINT(10 0)');

I would like to select one point closest to input - for example (1, 1)

It must be geographics distance using haversine formula (https://www.movable-type.co.uk/scripts/latlong.html).

What is the best way (fastest by performance) how can I do it?


texts: Custom TTF fonts in QGIS input app on iPhone

texts: Is it possible to use custom TTF fonts in QGIS input app on iPhone, or do I have to convert every single character to SVG instead?


texts: Replacing a feature class in a work version, or emptying it\Splitting features by mask?

texts:
I am working in a versioned database, in a specific feature class. After changes are made to this feature class (a copy outside of the versioned geodatabase) I want to be able to "update" the named feature class in the versioned database. I want this to happen with the click of a button (using Modelbuilder or Arcpy).

Further clarification:  I have a work version from the default and I "export" data from this by using clip. This clipped data is given to an external user to edit, change geometry and attribute table etc. I then get supplied with the new data and need to implement it into the work version, replacing old data with new data, before reconciling. The main issue is that almost everything done in ArcMap creates a new output.

So far I have found different options, but none of them do the job. These are the options I have tried so far:

Delete features and then append (This "works" but it takes over 3 hours, this is way too long)
Use Delete (This clears the entire database, I only need to replace 1 feature class with the same name)
The update tool (When polygons are "made smaller" this fills up the empty space with a polygon of it's own)

The delete features would work if I could split by mask (other feature class) and then delete the selection, but I do not know of any way to do this.

I am using ArcGIS 10.3.1


texts: Creating the straight line feature class representing the shortest path between two features using ArcGIS Desktop?

texts:
I have a feature class containing the district boundaries for India and a feature class of one of India's national highways. For each district, I calculated the straight line distance to the nearest highway segment using the Near tool (which calculates the distance between the highway to the nearest boundary of a given district).

How can I create the feature class representing the straight line distance created by the Near tool?

I am using ArcMap version 10.5


texts: Default Property/Attribute Template For Features On Specific Layers

texts:
The map application I'm currently designing allows users to create features on one or more vector layers. What I'd like to know is if it is possible to set default feature properties at the layer level so that any features created on the layer automatically have those properties/attributes.

For example if I setup a vector layer called 'House_Residential' when a user creates a feature on it it automatically has an attribute/value of 'house:residential'.

At the moment the only solution I can think of is to capture when the feature is added to the layer and give it a property/attribute at that point, not sure if this is the best way though.


texts:
The query I am running is taking a considerable amount of time. I was wondering whether I am missing an opportunity to drastically improve performance.

My goal is to find out how many neighboring houses a house has. I based my code on this example: How to use st_touches to get all nearby states .

ST_Touches seems to be the right function. The query seems to be correct. I just have a very large amount of rows, not like in the example where they use states. Any ideas? It has a spatial index.

with
panden as (select * from source.houses where pandstatus = 'Pand in gebruik' OR  pandstatus = 'Pand in gebruik (niet ingemeten)')
select
    pand1.identificatie, count(pand2.geovlak)
from
    panden as pand1, panden as pand2
where
    st_touches(pand1.geovlak, pand2.geovlak) AND
    pand1.identificatie <> pand2.identificatie
group by
    pand1.identificatie;


texts: PostgreSQL query to join tble1 with multiple table based on common id to validate features and attribute

texts:
I have 3 shapefiles imported into PostgreSQL with PostGIS to do query find out if any mismatch. So now I want left join the Table1 on to Table2 and Table3 on id retrieve two results

1.) Geomtery -
`Table2.geom` based on recent date if 2 or more id matches on `table2` .
`Table3.geom` if 1 id matches from Table2 and 1 id from Table3.

2.) Attribute -
`Itmtype` and other Attribute based on Older Date. however in both case I would concatenate Number of `sry_unique` and `SurveyCounts` as Number of item found. and

I tried below code but it still needed me code more so please suggest me how do I improve on it to achieve what I needed as below.

SELECT t1.id,
          count(t2.sry_unique),
          count(DISTINCT(t3.sry_unique)) ,

          string_agg(t2.sry_unique,', ' ORDER BY t2.date) AS uniqueref_att,

          string_agg(DISTINCT t3.sry_unique,', ') AS lidar_att,
          t1.geom
   FROM table t1
   LEFT JOIN table2 t2 ON t1.id = t2.id
   LEFT JOIN table3 t3 ON t1.id = t3.id
   GROUP BY t1.id

Output table should look smilar to this.

+----+--------------+----------+------------+----------+--------------+----------------+---------------------+------------------------+----------------------------+------------+----------+------------+-----------------+---------------+-----------+
| id | surveycounts | Itmtype  | sry_unique | tble2.id | tbl2.Itmtype | tble2.olderITM |  tble2.sry_unique   | tble2._RecentSryUnique | tble2.Olderecent_SryUnique | tbl2.Date  | tble3.id | tble3.Date | tble3.sry_uniqe |  Geom_match   | Attribute |
+----+--------------+----------+------------+----------+--------------+----------------+---------------------+------------------------+----------------------------+------------+----------+------------+-----------------+---------------+-----------+
|  1 |            1 | MG       | XXX1       |        1 | MG           | MG             | XXX1                | XXX1                   | XXX1                       | RecentDate |        1 | olderDate  | XXXXXTY1        | tble3geomtery |           |
|  2 |            2 | SG       | XXX2       |        2 | SX,SG        | SG             | XXX2, XX2_X         | XXX2                   | XXX2_X                     | RecentDate |        2 | olderDate  | XXXXXTY1        | tble2geomtery |           |
|  3 |            3 | MN       | XXX3       |        2 | SG,MN,MX,MYY | MN             | XXX3,XXX3_X,XXX3_XX | XXX3                   | XXX3_XX                    | RecentDate |        3 | olderDate  | XXXXXTY3        | tble2geomtery |           |
+----+--------------+----------+------------+----------+--------------+----------------+---------------------+------------------------+----------------------------+------------+----------+------------+-----------------+---------------+-----------+

Table1 - Main Data

+----+--------------+----------+------------+
| id | surveycounts | Itmtype  | sry_unique |
+----+--------------+----------+------------+
|  1 |            1 | MG       | XXX1       |
|  2 |            2 | SG       | XXX2       |
|  3 |            3 | MN       | XXX3       |
+----+--------------+----------+------------+

Table2 - Multiple Survey Data

+----+--+---------+------------+------------+--+
| id |  | Itmtype | sry_unique |    Date    |  |
+----+--+---------+------------+------------+--+
|  1 |  | MG      | XXX1       | 2015-11-04 |  |
|  2 |  | GN      | XXX2       | 2011-12-03 |  |
|  2 |  | SG      | XXX2_X     | 2010-06-03 |  |
|  3 |  | MN1     | XXX3       | 2009-05-03 |  |
|  3 |  | MXx     | XXX3_X     | 2009-01-02 |  |
|  3 |  | MX      | XXX3_XX    | 2008-02-16 |  |
+----+--+---------+------------+------------+--+

Table3 - LidarData

+----------+--------------+
| lidar_id | lidar_unique |
+----------+--------------+
|        1 | XXXXXTY1     |
|        2 | XXXXXTY1     |
|        3 | XXXXXTY3     |
+----------+--------------+


texts: Moving file flat-nodes.bin from an osm2pgsql import

texts:
I'm running PostgreSQL 9.4 and osm2pgsql 0.82.0 on Ubunt 14.04 (64bit).

I'm new-ish to postgres and less new to OSM. My goal was to import the planet file.

Unfortunately, when I started the import, I made the rookie mistake of not starting the osm2pgsql command from a sensible point in the filesystem.

Since I imported the planet file with the --flat-nodes option, the result (of course!) is that the flat-nodes.bin file is located in the folder I launched the command from e.g. /home/foo/bar.

I would really like to correct this and move flat-nodes.bin to a better location.

Is there a way to change where Postgres looks for the file?  I would really rather not reimport the planet.


texts: How to create high quality map with QGIS and Stamen tiles

texts:
I'm wanting to create a simple but high quality map in QGIS with two minimal vector layers and Stamen's landscape tiles as a basemap. When viewing my map in the QGIS print composer everything looks great, but a .png export of the map from the print composer looks a lot different. I know that this issue is somewhat common and I'm not dead set on using the QGIS print composer.

I've tried saving the map as a .png straight from QGIS (not in the print composer), and this looks pretty good but obviously does included a needed legend, scale bar, etc. I've tried importing this .png into Inkscape to manually include map elements, but the resolution always looks reduced in Inkscape. Can I use QGIS, Stamen tiles, and Inkscape to create a high quality map, or would other tools better fit this task?


texts:
I'm using the Interpolation plug-in in QGIS 2.18 to perform IDW Interpolation. I was hoping to specify a power of 2 and use the 10 nearest points to calculate the output, but I don't necessarily see a way to do so.

Is there another tool that I can use to accomplish this?


texts: Fresh QGIS 3.12 and 3.10 Installation result in unknown error when saving Spatialite layer

texts:
This is the most interesting problem to resolve.
On a brand new Windows 10 HP Laptop, I install Windows with the latest updates and then install both versions of QGIS 3.10 and 3.12.
On creating a new Spatalite layer, I add a few Point type and click Save.
The system returns with Unknown cause:

SQLite errors: unknown cause
SQL: INSERT INTI "test"("geometry",)VALUES (GeomFromWKB(?4326),)

I imported a GeoJSON file with the Point geometry and this worked as expected.
The final test is to create a row using the DB Browser for Spatialite.  This returned "error adding record no such function GeometryConstraints" which is unexpected.

Any guidance to understand what is wrong on a PC that has no other software?

Thanks


texts: Is there a procedure to determine optimal location viewpoints with QGIS?

texts: I must visually cover a 7,000 ha area from observation points. Does anybody know a procedure using QGIS to determine their locations, minimizing the number of points and maximizing the viewshed?


texts:
I am writing my first Python plugin. There was a problem during work. The plugin is designed to select a raster and calculate the statistics for the area bounded by the vector. I have created a list of rasters, but I do not know how to refer to the selected raster in the further part of the plugin to calculate statistics for it. The following code is created.

List of raster layers:

def create_layer_list(self):
    ''' Build list of layers'''
    self.dialog.chooseCombo.clear()
    self.layers = list(self.registry.mapLayers().values())
    #self.layers = self.iface.legendInterface().layers()
    layer_list = [""]+[layer.name() for layer in self.layers]
    self.dialog.chooseCombo.addItems(layer_list)

Main function for counting statistics:

def get_statistic(self):

      index2 = self.dialog.chooseCombo.currentIndex()
      rlayer = self.layers[index2].layer_list()
      provider = rlayer.dataProvider()

      for band in range(1, rlayer.bandCount() + 1):
        stats = provider.bandStatistics(band, QgsRasterBandStats.All, self.extent, 0)
        min = stats.minimumValue
        max = stats.maximumValue
        mean = stats.mean
        statsGathered = stats.statsGathered
        results= print ('min: {:.2f}, max: {:.2f}, mean: {:.2f}, statsGathered: {:.2f}'.format(min, max,mean,statsGathered))
        self.dialog.textlabelLayout.setText(results)

Unfortunately, Python throws an error when indexing ..


texts: Filtering layer in Query Builder with array from other layer using QGIS

texts:
I have the following situation using QGIS 3.4:
'Layer 1' with lines and a field ('Field X') and 'Layer 2' with points where values in a 'Field X' are available.
Now I would like to filter all points from 'Layer 2' for which the value in 'Field X' exists also in Layer 1.
I tried the following principle in query builder -> extraction of an array of all values in 'Field X' from 'Layer 1' and check whether value of 'Layer 2' is included:
array_contains(aggregate('layer id Layer 1','array_agg',"Field X"),"Field X")
But this does not filter any points and displays them all. Whereas the same (!) formula does work perfectly fine in attributes table (advanced filter (expressions)) and also for the rule-based symbolizing -> true if value of 'Field X' in 'Layer 2' exists in 'Layer 1' and false if not.
How does the formula need to be changed to work also in query builder?


texts: Establishing a MySQL connection in ArcGIS 10

texts: Is it possible to create a connection to a MySQL database in ArcGIS 10?


texts:
I'm using Geoserver/PostGIS/openLayers

I've succesfully set a vector layer. And now I'm trying to create a search functionality.
But something goes wrong.

First I set up the layer and add some basic style for the features. So far, everything is ok.

Then I create a function, using javascript, which comes after the layer. This function gets triggered after the user completes a HTML form and hits the "Search" button.

Here's what's inside the function:

//check for empty fields in the form-if not pass the value in a var

if (document.getElementById("typeselect").value!=" ") {var filter_a = document.getElementById("typeselect").value;}

if (document.getElementById("eraselect").value!=" ") {var filter_b = document.getElementById("eraselect").value;}

if (document.getElementById("nm").value!="") { var filter_c = document.getElementById("nm").value;}

if (document.getElementById("descr").value!="") { var filter_d = document.getElementById("descr").value;}

if (document.getElementById("chrono").value!="") { var filter_e = document.getElementById("chrono").value;}

//create a style

var vsd = new OpenLayers.Style();

//put together the final filter, using the vars from above

//filter_a and filter_b are numbers

//filter_c to filter_e are texts

var filter_final = new OpenLayers.Filter.Logical({ type: OpenLayers.Filter.Logical.AND,

filters: [ new OpenLayers.Filter.Comparison({

type: OpenLayers.Filter.Comparison.EQUAL_TO, property: "p_type", value: filter_a }),

 new OpenLayers.Filter.Comparison({ type: OpenLayers.Filter.Comparison.EQUAL_TO, property: "p_era", value: filter_b }),

new OpenLayers.Filter.Comparison({ type: OpenLayers.Filter.Comparison.LIKE, property: "p_name", value: filter_c }),

 new OpenLayers.Filter.Comparison({ type: OpenLayers.Filter.Comparison.LIKE, property: "p_descr", value: filter_d }),

new OpenLayers.Filter.Comparison({ type: OpenLayers.Filter.Comparison.LIKE, property: "p_chron", value: filter_e }) ] });

//add the filter as rule to the style

vsd.addRules([filter_final]);

//apply the style as default

var vst = new OpenLayers.StyleMap({'default': vsd});

//apply the new style to the layer (named pins)

pins.styleMap = vst;

This is not working at all. I've tried editing my code, I don't know what to do, I've searched for an answer and nothing really helped me. Please advise...

Thanks.

EDIT

I also added this line at the end and still nothing pins.refresh();
Do I need to use the WFS protocol, combined with the read() ?

EDIT #2
I've replaced pins.refresh(); with pins.refresh({force:true}); and all the features dissappear. I also tried something like this and still nothing. This is a torture...


texts: Tools in the Sextante toolbox in QGIS 1.8 for interpolation?

texts:
What tools in the Sextante toolbox in QGIS 1.8 can be used for interpolation (kriging and spline)?

I have looked at every single one and could not figure it out.


texts: OpenLayers featureclick event does not work with named graphic (graphicName)?

texts:
OpenLayers 2.13. (and 2.13.1) I have been trying to make featureclick event working with vector layers and named graphic features. It works with externalGraphic, but not named graphic (graphicName: ...).

Although it does work in IE8, but not IE9, IE10 or IE11. Does not work in FF and Chrome neither. See this example: http://jsfiddle.net/rwgHY/

The triangle is a named graphic feature and clicking it does not do anything (this is the problem). The marker is a PNG graphic (external graphic) and when clicked, the event is triggered and the alert shown.

The event is hooked on the map level. Features are placed on two different vector layers.

I've been looking for solutions on the net, but did not find anything. Is it something silly I am missing here?


texts:
In various screen grabs of the dbeaver spatial view such as the blog post Geometry viewer added to DBeaver a basemap is shown for the spatial data.

How do I add a basemap to OSM?


texts: Escaping single quotes for an SQL query in ArcPy

texts:
I have some data with strings containing single quotes (') and can't figure out a way to pass that string into an SQL query in arcpy.analysis.Select().
The SQL query running on a feature class in a GDB requires single quotes for values e.g. Query: "Name = 'Bob'". I want to pass in a name such as "exam'ple2" into the query.
The "Definition Query" in ArcGIS Pro resolves this by using '' (2x single quotes). See below for an example of a working query:
NAME IN ('example1','exam''ple2','''examp''le3')

--- each ' in input is replaced with ''
However, this is not a working option for a python script since each pair of quotes (' ') is treated as a separate string, and all my attempts end up with strings encapsulated by double-quotes.
Here's the code I'm trying to run:
import arcpy

nametuple = ('example1',"exam'ple2","'examp'le3")

arcpy.analysis.Select("Areas","test",f"NAME IN {nametuple}")

ERROR 000358: Invalid expression NAME IN ('example1',"exam'ple2","'examp'le3")


texts: Reset the "Auto open form for single feature results"

texts: Is there a way to remove the setting "Auto open form for single feature results" in QGIS 3.10?


texts: geodango form with geometry field and leaflet map

texts:
i use geodjango and leaflet and i want to create a django form where including geometry field with map widget  (the user must to put in leaflet map) like as django admin page . i try to follow this example

but in html page dont show me the form ,any idea ?

here the code :

urls.py

url(r'^$', TemplateView.as_view(template_name='index.html'), name='home'),

models.py

from djgeojson.fields import PolygonField
from django.db import models

class MushroomSpot(models.Model):
    name = models.CharField(max_length=256)
    geom = PolygonField()

    def __unicode__(self):
    return self.name

forms.py

from django import forms
from leaflet.forms.fields import PolygonField
from mushrooms.models import MushroomSpot

class WeatherStationForm(forms.ModelForm):
    geom = PolygonField()
    class Meta:
        model = MushroomSpot
        fields = ('name', 'geom')

html page :

{% load leaflet_tags %}

{% leaflet_js plugins="forms" %}
    {% leaflet_css plugins="forms" %}

Edit {{ object }}

<form action="POST">
{{ WeatherStationForm }}
<input type="submit"/>
</form>


texts: Concave hull error: not enough arguments

texts:
I am trying to run concave hull on a points layer, but I keep running into this error:

  TypeError: translate(str, str, disambiguation: str = None, n: int =
  -1): not enough arguments

What am I doing wrong?

I'm using QGIS 3.0.


texts: Where is uninstalled plugin tab in Plugin window of QGIS?

texts:
I was trying to install mmqgis plugin for QGIS 2.18. However, i cannot see the 'uninstalled' tab in the plugin window. I need help on:

how to fix this issue?
how to install MMQGIS plugin?


texts:
I am designing a spatial database and I am including a variety of already existent datasources (ex: Lakes, rivers, road network, contour lines, cities, etc). Those sources can be updated every once in a while.

Except from these sources, the database is normalized in fifth normal form (5NF) or Boyce-Codd normal form (BCNF).

Should I also normalize those external datasources?

My concern is that it may be too much trouble for the benefits. I am thinking about the moment when I will update the sources. It would be much more trouble than to just "replace" the existing table.


texts:
I have created a shapefile of polylines that represent a path with multiple intersections. Then I have used the networkx python library to create a network. At the first the results were satisfactory.

But when I tried to use the shortest path plugin of QGIS it returned a result like:

What can I do to fix this?


texts: Python ImportError with libgdal, undefined symbol: opj_codec_set_threads

texts:
GDAL and its Python 3 bindings have been installed through apt package manager for weeks, but, suddenly:

from osgeo import GDAL
---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-1-61e0eefce4de> in <module>()
----> 1 from osgeo import gdal

/usr/lib/python3/dist-packages/osgeo/__init__.py in <module>()
     19                 fp.close()
     20             return _mod
---> 21     _gdal = swig_import_helper()
     22     del swig_import_helper
     23 else:

/usr/lib/python3/dist-packages/osgeo/__init__.py in swig_import_helper()
     15         if fp is not None:
     16             try:
---> 17                 _mod = imp.load_module('_gdal', fp, pathname, description)
     18             finally:
     19                 fp.close()

/usr/local/pythonenv/smap_l4c-env/lib/python3.6/imp.py in load_module(name, file, filename, details)
    241                 return load_dynamic(name, filename, opened_file)
    242         else:
--> 243             return load_dynamic(name, filename, file)
    244     elif type_ == PKG_DIRECTORY:
    245         return load_package(name, filename)

/usr/local/pythonenv/smap_l4c-env/lib/python3.6/imp.py in load_dynamic(name, path, file)
    341         spec = importlib.machinery.ModuleSpec(
    342             name=name, loader=loader, origin=path)
--> 343         return _load(spec)
    344
    345 else:

ImportError: /usr/lib/libgdal.so.20: undefined symbol: opj_codec_set_threads

The same error, with slightly different Traceback, is obtained with import gdal.

I didn't recognize this symbol name, opj_codec_set_threads so I looked it up online. It's apparently part of the OpenJPEG library. I recalled that I had recently run some legacy setup script for another project with OpenJPEG as a dependency. That legacy setup script apparently installed an older verison of OpenJPEG (v2.1.0) on top of what apt had installed (v2.3.0) when I first setup the machine.

I tried reinstalling OpenJPEG from apt:

sudo apt install --reinstall libopenjp2-7 libopenjp2-7-dev
sudo ldconfig

This did not solve the problem. I tried re-installing GDAL and its bindings:

sudo apt install --reinstall gdal-bin libgdal-dev python-gdal python3-gdal

Still no good.


texts: Remove polygons within other polygons in a dataframe

texts:
Given a dataframe, is it possible to find polygons that lie within another polygon in the same dataframe and remove them based on a condition?
For example, consider the following dataframe
                                            geometry          category
0  POLYGON ((332294.432 2794730.319, 332294.032 2...            A
1  POLYGON ((332193.232 2794730.719, 332192.832 2...            B
2  POLYGON ((332271.632 2794731.519, 332271.232 2...            C
3  POLYGON ((332271.632 2794731.519, 332271.232 2...            C

Suppose polygon 3 lies within polygon 2. Since they are both of category C, I would like to remove polygon 3. If they were of different categories, I would not do a removal.
Is there a way to do this within geopandas?


texts:
Are there strategies for dealing with TopologyException errors in PostGIS 2? Nearly every function I run against two geometries results in a TopologyException (e.g. ST_Intersection, ST_SymDifference).

I always try ST_Buffer(the_geom, 0) first, then ST_MakeValid(), then (ST_Dump(the_geom).geom and repeat ST_Buffer. Sometimes a combination of one or more fixes it, sometimes not.

Perhaps my data is just not clean enough for PostGIS 2? I feel like I rarely had this issue on Postgres 9.0 with PostGIS 1.5, but now it is almost every dataset.


texts:
I have created image collection in GEE of LANDSAT 7 NDVI (LANDSAT/LE07/C01/T1_8DAY_NDVI). I chose date range and I know I have 71 images in the image collection.

Is there any way to clip all the image collection according to a shapefille I have and see them as layers? Should I build function that calls every image and visualize it?

I used filterBound and the code run but when I try to visualize even just one image it doesn't really filter it according to the shape. My end goal is to download all the NDVI images I have for that period clipped according to the shapefille


texts: Converting XY coordinates in the Iberian peninsulate to lat, lon?

texts:
I have been looking for a way to convert XY coordinates (from a coordinate system based on ETRS89) to latitude and longitude coordinates.

So far I have this function:

from pyproj import Proj, transform
def utm_to_lon_lat(easting, northing):

    """
    ETRS89 to latitude, longitude converter
    :param easting: Easting coordinates array
    :return: latitude, longitude arrays
    """

    n = len(easting)
    latitude = np.zeros(n)
    longitude = np.zeros(n)

    inProj = Proj(init='epsg:3857')
    outProj = Proj(init='epsg:4326')

    for i in range(n):
        if not np.isnan(easting[i] + northing[i]):
            latitude[i], longitude[i] = transform(inProj, outProj, easting[i], northing[i])
            print(easting[i], northing[i], '->', latitude[i], longitude[i])

    return latitude, longitude

However, it is not working. The error is that not all the points fall into what I know is the area of the domain, The Iberian peninsula. My guess is that I got the projections wrong, is no it?


texts: QGIS plugin output layers not displaying?

texts:
I'm new to using qgis and I am trying to run the transectizer plugin and the feature grid creator.  They basically plot bore holes within selected parameters.  I run the tools and make the right selections.  A layer is added to my table of contents, but it does not display in my window.  I open up the attribute table and it hasn't been populated with the bore holes either.

Can anyone help me understand what I'm doing wrong?


texts: Projection of all the polygons contained in a geojson at once

texts:
I would like to project all the features in a geojson from WGS84 (EPSG:4326) to Web Mercator (EPSG:3857) using python.

This geojson contains only polygons.

I found out the following shapely code combined with pyproj will do it. However after the projection has been applied then I need to manually build the output geojson with the new Web Mercator projection.

import pyproj
import json

from functools import partial

# on windows: conda install -c conda-forge shapely=1.5.17
from shapely.geometry import Polygon, Point, shape, mapping
from shapely.ops import transform

# the shapely projection function
datum_wgs84 = pyproj.Proj(init='EPSG:4326')
projection_web_mercator = pyproj.Proj(init='EPSG:3857')
projection_wm_func = partial(pyproj.transform, datum_wgs84, projection_web_mercator)

# load GeoJSON file containing sectors
with open('my_input_file.geojson') as f:
  js = json.load(f)

# project GPS lat/lon coordinates to web mercator for each polygon in the geojson
for feature in js['features']:
  polygon = shape(feature['geometry'])
  print('{} bounding box (in WGS84): {}'.format(feature['properties']['my_polygon_name'], polygon.bounds))
  print('-> A sample pair in the original WGS84 projection: {}'.format(mapping(polygon)['coordinates'][0][0]))
  projected_polygon = transform(projection_wm_func, polygon)
  print('-> A sample pair in the projected "Web Mercator": {}'.format(mapping(projected_polygon)['coordinates'][0][0]))

# TODO put all the projected polygons together in a new output geojson

Is there a more compact way of doing this manual part to avoid most of the boiler plate code?
Ideally some sort of function with inputs:

input geojson filename
output geojson filename
existing projection
final projection

Or at least a combination of functions to avoid the for loop.


texts: Can you rename the Django GeoManager object instance?

texts:
Is it possible to change the name of the GeoManager object from "objects" to something like "geo_objects"?  The way I've done it in the code below is supposed to be the default way to configure your model if you're using GIS.  I would think that you could so long as you always call the GeoManager using the new name.  However, I don't see this GeoManager object instance used elsewhere in my model so I'm not sure if it's possible to change the name if I can't see where it's being referenced.

from django.contrib.gis.db import models
class PersonalProfile(CommonProfileInfo):
    # ...
    # GeoDjango-specific fields
    home_point = models.PointField(srid=4326)
    objects = models.GeoManager()   # <= Change name of "objects"?


texts: I want to calculate a area of NDWI for a region for each month over for 10 years in Google Earth Engine

texts:
I want to calculate a area of NDWI for a region for each month over for 10 years in Google earth engine. I have a starting date and Ending date array for every month for each year.Using Image collection function I need calculate NDWI and its area for the image availability month and add it in a array. If image not available for a month then add zero in the array and so on.

Here is the Link for my code: https://code.earthengine.google.com/322b02e03dffed8d02327eef3474c1d9

function comps(yearArray, SD, ED) {
  var collection = [];
  for (var idx = 0; idx < SD.length; idx++) {
    var year = SD[idx].split('-')[0];
    var month = SD[idx].split('-')[1];
    var Date = year + '-' + month;
    timeseries.Year.push(Date);
    if (year <= 2013) {
      var SD_1 = SD[idx];
      var ED_1 = ED[idx];
      collection = colLT5.filterDate(SD_1, ED_1)
        .filterBounds(ER)
        .filter(ee.Filter.lt('CLOUD_COVER', 10));
      var Count = collection.size();
      Image_count.push(Count);
      if (Count.gt(0)) {
        Map.addLayer(collection.mosaic().clip(ER), vizParams_L5, 'FCC_' + year, false);
        var NDWIcomp_L5 = collection.map(NDWI_L5).mosaic().clip(ER);
        var L5_Aqua = NDWIcomp_L5.updateMask(NDWIcomp_L5.gte(-0.1));
        Map.addLayer(L5_Aqua, {
          bands: 'NDWI',
          min: -1,
          max: 1,
          palette: ['#eff3ff', '#bdd7e7', '#6baed6', '#3182bd', '#08519c']
        }, 'NDWI_' + year, false);

        var Water_Count_L5 = Pixel_Count(L5_Aqua);
        var ndwi_area_L5 = (Water_Count_L5.multiply(900)).divide(1000000);
        timeseries.Area.push(ndwi_area_L5);
      } else {
        timeseries.Area.push(0);
      }

This is the relevant part of code. If the collection.size()>0 then calculate NDWI and area. Push the area in the Area array otherwise zero value to be pushed in the area for that date


texts: I have been trying to import a CSV file into QGIS 3 on Linux. I end up with all my points at 0 lat and 0 long, off the coast of Africa. They should be in Florida.  I cannot find this same issue in the forums. I have specified my X and Y, defined the geometry, and switched these around. I tried importing as degrees, minutes, seconds. I get the message that this layer has no projection specification, but I have one selected? What do I need to do?


texts: Layer.bindTooltip doesn't work in Leaflet

texts:
My objective is to plot polygons and assign a tooltip label to each one. For some reason, my code doesn't seems to work.

poly_geojson = L.geoJson(poly_gon,{
    onEachFeature(feature, layer) {
        if (feature.properties.sourceLang == "en") {
            layer.bindTooltip("
"+feature.properties.name+"
")
        }
        else {
            layer.bindTooltip("
"+feature.properties.name+"
")
        }
    }
}).addTo(map);


texts: Using Select by Attribute within a for loop

texts:
I am trying to iterate through a feature layer of counties, using each row as the selection expression for Select by Attribute. This selection will then be used to clip another layer. The goal is to break the initial layer into smaller pieces to aid in geoprocessing.

`import arcpy
counties = arcpy.MakeFeatureLayer_management(r"R:\Data\Base_Data\Administrative\Admin_Political_Bounds.gdb\HL_Political_Boundaries\HL_Counties_NJGIN", "counties")
workspace = r"R:\Data\Data_Development\HDC_Allocations\HDC_raw3.gdb"
with arcpy.da.SearchCursor(counties, ['SHAPE@AREA','COUNTY']) as cursor:
    for row in cursor:
        i=row[1]
        expression = '"COUNTY" = {0}'.format(i)
        county = arcpy.SelectLayerByAttribute_management(counties,"NEW_SELECTION", expression)
        outname = workspace + "_" + row[1]
        arcpy.Clip_analysis(HDC_raw, county, outname)
        arcpy.Delete_management(county)`

My problem is that I keep getting a syntax/invalid expression error on the Select Layer by Attribute line.

This is part of a much larger script that is run as a standalone script.


texts: How to print only a selection of the map using MapFish Printing Module?

texts:
I am working on a GIS API that uses OpenLayers 2. I want to make snapshots of a selected region. For example. I draw a rectangle, and the MapFish Printing Module would print only the selected region.

I have been trying to use the bbox property, but I get the wrong image, i.e., instead of a rectangle near Lisbon, I get a rectangle with south Europe and North/Center Africa.

Here are my print options:

var printOptions = {
    units: 'm',
    srs: 'EPSG:900913',
    layout: options.layout,
    dpi: '300',
    mapTitle: options.mapTitle,
    comment: options.comment,
    layers: [
        {
            baseURL: 'http://a.tile.openstreetmap.org/',
            opacity: 1,
            singleTile: false,
            type: 'OSM',
            maxExtent: [-20037508.3392, -20037508.3392, 20037508.3392, 20037508.3392],
            tileSize: [256, 256],
            extension: 'png',
            resolutions: [156543.03390625, 78271.516953125, 39135.7584765625, 19567.87923828125, 9783.939619140625, 4891.9698095703125, 2445.9849047851562, 1222.9924523925781, 611.4962261962891, 305.74811309814453, 152.87405654907226, 76.43702827453613, 38.218514137268066, 19.109257068634033, 9.554628534317017, 4.777314267158508, 2.388657133579254, 1.194328566789627, 0.5971642833948135]
        }
    ],
    pages: [
        {
            bbox: [-1050091.8944356, -990623.88643799, 4675576.582525099, 4733515.8499572]
        }
    ]
};

I am assuming that the bbox array is defined like this: [minX, maxX, minY, maxY].

Is this correct?

I managed to make it work using Thomas B's answer.

I draw a box, I get the bounds of it, and then I get the BBOX value in an array, like this:

var printBbox = geometry.getBounds().toArray();

And the, in the printing module definitions, I do this:

pages: [
    {
        bbox: printBbox
    }
]


texts: Iterating model that solves multiple routes with multiple locations?

texts:
I have multiple line feature classes that I would like to make routes and solve. For each of these routes I have specific start/stop points (feature classes) that will act as locations.

How could I build a model that would pull each route with its associated locations, solve, and write out the solved route?


texts: Converting DEM in degree decimal to meter using ArcMap?

texts:
i am trying to make cell size of satellite image and DEM equal for Maxent.

cell size of satellite image is in metric but DEM is in degree decimal.

how can i convert DEM in degree decimal to meter using ArcMap?


texts: arcpy.da.InsertCursor() Fails with SDE But Works with Local FGDB

texts:
I am trying to use arcpy.da.InsertCursor() to insert rows to an SDE table but no rows were added and without any error message. I tried the same code with an exact replica table in a local FGDB and it worked fine. The code is simple:
workspace = r"Database Connections\SDE.sde"
arcpy.env.workspace = workspace
arcpy.env.overwriteOutput = True

edit = arcpy.da.Editor(workspace)
edit.startEditing(False, True)
edit.startOperation()

cursor = arcpy.da.InsertCursor(sde_table, ["Date_"])
cursor.insertRow(["10/21/2020"])

I also tried to use arcpy.InsertCursor() and it works well in both SDE and local FGDB. Is there anything missing for getting arcpy.da.InsertCursor() to work in SDE?
RDBMS: SQL Server 2012
Python: 2.7
SDE: 10.6.1


texts: Count unique features of points in polygon

texts:
I have a series of points with various fields in the attribute table such as car_id, time etc.,

I've created a hexagon vector grid, I want to count how many unique car_id's occur in one of the hexagon polygons.

Essentially a step further than points in polygon, if there are multiple points in a polygon where car_id = 25, the car_id = 25 would be counted once. The car_id count for a polygon using the sample data below would be 4

car_id = 25
car_id = 25
car_id = 25
car_id = 22
car_id = 22
car_id = 23
car_id = 24

My question is what is the most appropriate QGIS tool to achieve this? Would it be Count Unique Points in Polygon, in the processing toolbox?


texts: How to update url parameters (CQL_FILTER) in OpenLayers 3

texts:
I have a vector layer setup with a CQL Filter:

vectorSource = new ol.source.Vector({
                loader: function (extent, resolution, projection) {
                    var url = 'http://server01:8080/geoserver/wfs?service=WFS&' +
                        'version=1.1.0&request=GetFeature&' +
                        'typename=sf:TABLE01&' +
                        'CQL_FILTER=ID=355&' +
                        'outputFormat=text/javascript&' +
                        'format_options=callback:gotFetch';
                // use jsonp: false to prevent jQuery from adding the "callback"
                // parameter to the URL
                $.ajax({
                    url: url,
                    dataType: 'jsonp',
                    jsonp: false
                });
            }
        });

I want to have the user press a button and update the id above to 455 (I want to refresh the layer with latest CQL_FILTER data):

CQL_FILTER=ID=455

And then refresh the layer.  How can I refresh the WFS Layer?

(Also I need a wfs layer because I am setting up markers, I cannot use a WMS layer)

This did not work (gotFetch was not called):

var params = {
                CQL_FILTER: 'ID=455'
            };
vectorSource.updateParams(params);


texts: Error in multi-temporal filtering of images

texts:
I tried multi-temporal filtering of the Sentinel-1 images by adopting the code available at https://code.earthengine.google.com/01e247bc0c59442437a737cfc28f7430

I need to apply the multi-temporal filtering for the images of whole year (2017). However, I got this error:

ImageCollection (Error)
Error in map(ID=S1B_IW_GRDH_1SDV_20170101T011828_20170101T011857_003648_006414_E48C):
Image.multiply: If one image has no bands, the other must also have no bands. Got 1 and 0.

Here is my code:

var aoi = ee.Geometry.Polygon(
        [[[-113.04179687499999, 44.484512810104995],
          [-113.04179687499999, 40.26669922631106],
          [-105.13164062499999, 40.26669922631106],
          [-105.13164062499999, 44.484512810104995]]], null, false);

function multitemporalDespeckle(images, radius, units, opt_timeWindow) {
  var timeWindow = opt_timeWindow
  var bandNames = ee.Image(images.first()).bandNames()
  var bandNamesMean = bandNames.map(function(b) { return ee.String(b).cat('_mean') })
  var bandNamesRatio = bandNames.map(function(b) { return ee.String(b).cat('_ratio') })
  var meanSpace = images.map(function(i) {
    var reducer = ee.Reducer.mean()
    var kernel = ee.Kernel.square(radius, units)
    var mean = i.reduceNeighborhood(reducer, kernel).rename(bandNamesMean)
    var ratio = i.divide(mean).rename(bandNamesRatio)
    return i.addBands(mean).addBands(ratio)
  })

  function multitemporalDespeckleSingle(image) {
    var t = image.date()
    var from = t.advance(ee.Number(timeWindow.before), timeWindow.units)
    var to = t.advance(ee.Number(timeWindow.after), timeWindow.units)
    var meanRatioBefore = ee.ImageCollection(meanSpace).select(bandNamesRatio).filterDate(from, t)
      .filter(ee.Filter.eq('relativeOrbitNumber_start', image.get('relativeOrbitNumber_start')))
    var meanRatioAfter = ee.ImageCollection(meanSpace).select(bandNamesRatio).filterDate(t, to)
      .filter(ee.Filter.eq('relativeOrbitNumber_start', image.get('relativeOrbitNumber_start')))
    var SDmeanRatioBefore = meanRatioBefore.reduce(ee.Reducer.stdDev());
    var SDmeanRatioAfter = meanRatioAfter.reduce(ee.Reducer.stdDev());
    var b = image.select(bandNamesMean)
    var monoLookBack = b.multiply(meanRatioBefore.sum()).divide(meanRatioBefore.count()).rename(bandNames)
    var biLookBackForward = monoLookBack.where(SDmeanRatioBefore.gt(SDmeanRatioAfter), b.multiply(meanRatioAfter.sum()).divide(meanRatioAfter.count()).rename(bandNames));
    return biLookBackForward;
    }
    return meanSpace.map(multitemporalDespeckleSingle).select(bandNames)
}

var s1 = ee.ImageCollection('COPERNICUS/S1_GRD')
  .filterBounds(aoi)
  .filterDate('2017-01-01','2017-12-31')
  .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
  .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
  .filter(ee.Filter.eq('instrumentMode', 'IW'))
  .filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'));

print (s1);

var s1vv = s1.select('VV')
  .sort('system:time_start', true);

print (s1vv);

var s1vh = s1.select('VV')
  .sort('system:time_start', true);

var s1Denoised_vv = multitemporalDespeckle(s1vv, 7, 'pixels', { before: -8, after: 8, units: 'week' });
var s1Denoised_vh = multitemporalDespeckle(s1vh, 7, 'pixels', { before: -8, after: 8, units: 'week' });

print (s1Denoised_vv);


texts: How can I change the line width of the outline for polygons in the Mapfile?

texts: I display Polygons in the UMN Mapserver. With the OUTLINE-keyword I set a color for the outline of the polygons. How can I change the line width of this outline?


texts: Getting image band statistics in Google Earth Engine code

texts: What is the command to get the band statistics (min/max/etc.) for an image? I know I can do it for a region (using image.reduceRegion) but I want to have it for the whole image without defining a region.


texts: ArcGIS Raster Calculator Erase raster from raster

texts:
Two rasters

slope (black and white grid with slope values)
rasterized property (red, values do not matter)

I want to erase the red raster from the slope raster

I found this thread

https://community.esri.com/thread/190381-how-to-perform-reverse-clip

and they suggested something like this (for my case block is the red raster and buff is the slope)

Con(IsNull("block.tif"),"buff.tif")

this produces the output in the multi colors. But it is not the full erase, part of the buff/slope is missing

I do not want to go the vector conversion route.

Answer needs to use Raster Calculator syntax


texts: Scaling layer symbols from ArcMap to QGIS

texts:
I am recreating the symbology of a cartography moving from ArcMAP to QGIS substituting some font symbols in arcmap and reproducing them with new  .svg symbols. The problem is that setting the size of symbols according to ArcMAP (points, millimeters) in QGIS the dimensions of the same symbols at the same displaying scale are very different (larger). The definition (and printing) scale is 1: 10,000 and map unit is meter. I'm confused about a way on how to match the same scale of symbols between ArcMAP and QGIS.
Update:


texts: Showing link for every city at different zoom levels using Leaflet?

texts:
I have a Leaflet map that I am using to render a map of the U.S. The map will have specific zoom levels. For example, when the map first loads it will show major U.S. cities. Then if you zoom in it will show more cities on the point where you zoomed. I already have this set up in Leaflet using Mapbox for tiles.

The problem is that I need to be able to show a link on the map for every city - only on the cities that are visible on the map at the current zoom level.

Is there any way for me to know which cities are shown at the current zoom level? If not, is what I'm trying to accomplish possible with d3 or any other libraries or plugins?


texts: R: aggregate raster with 'mode' function - how does it work?

texts:
I want to aggregate a very fine resolution raster at various coarser resolutions (that represents many 1000s of categorical polygons).

The aggregate{raster} function is what i've traditionally used, the factor arguement reduces the resolution while the 'fun' allows me to sum, max or whatever.

However, when using fun=modal, how does aggregate behave? If I simply reduce the resolution by a factor of 2 (i.e. quarter the total number of cells) how does aggregate work with new multimodal cells?

example code:

ras <- raster(nrows=10, ncols=10,xmn=0, xmx=10, ymn=0, ymx=10)
ras[1:50] <- 1
ras[51:100] <- 2
x1 <- c(9,10,20,34,54,56,67,88,99)
ras[x1] <- NA
plot(ras)
grid(nx = 10, ny = 10, col = "lightgray", lty = "dotted",lwd = 1)

ras2 <- aggregate(ras,fact=2,expand=FALSE,fun=modal,na.rm=T)
x11()
plot(ras2)
grid(nx = 10, ny = 10, col = "lightgray", lty = "dotted",lwd = 1)

Ideally, when reducing by a factor, if there is a multimodal result I'd like aggregate to randomly assign the new cell one of the modal values and not always choose the same (if that's indeed what it does). Larger factors such as 5 etc are less likely to cause this issue but it is still entirely possible.

FURTHERMORE: in the aggregate function above, na.rm is set to TRUE because i dont want every single NA cell to set the new, coarser, cell automatically to NA. However in the top right of the data, that new cell should really be NA as it has 3 smaller NA values. Do i just need to re-assign NA to 999 before running? and how will it work if there are 2 NA cells and 2 cells of another value?


texts:
I get this error every time I try to get the longest flowpath. It says  that the spatial reference is not set for the input(s), yet all layers in my table of contents are set in the same spatial reference which is WGS_84_UTM_zone_51N.

Does anyone know how to solve this problem?


texts: How to Select "Inverse" of intersection of two polygon layers?

texts: I'm wondering if there's a way to intersect two polygon layers, take the intersection, and remove that from one of the polygon layers, making a new layer out of whatever is left over? Basically, the idea I have is one of "punching holes" in a map of Africa where the holes are points around which I construct some sort of buffer. I have QGIS and ArcGIS 10.0. Is there a built-in or user created script for this?


texts: Finding quality flag in Sentinel-5P data to filter by

texts:
I would like to filter S-5P data to show only data with clouds less than 5% using Google Earth Engine. I'm working with Sentinel-5P NRTI NO2: Near Real-Time Nitrogen Dioxide.

A guide from Copernicus states that there are "quality flags related to the presence of clouds within each pixel: 0% pixels means that the sky above this particular area is free of clouds and, conversely, 100% means that the entire sky is covered by clouds." There are many cloud related S-5P parameters on p140 of the documentation

How can I filter by the quality flag mentioned using earth engine?

With Landsat data you can view the cloud property like this

// Get a specific metadata property.
var cloudiness = image.get('CLOUD_COVER');
print('CLOUD_COVER: ', cloudiness); // ee.Number

I've tried viewing the metadata properties of S-5P like this, but the properties returned are only high-level

// Get a list of all metadata properties.
var properties = collection.propertyNames();
print('Metadata properties: ', properties); // ee.List of metadata properties

I've also looked at the earth engine page Image Properties tab however those don't seem to relate to Quality or Cloud cover

Here is the basic code

var collection = ee.ImageCollection('COPERNICUS/S5P/NRTI/L3_NO2')
  .select('tropospheric_NO2_column_number_density')
  .filterBounds(geometry)

var current = collection.filterDate('2019-03-20', '2019-03-27').mean();

var vizParams = { min: 0, max: 0.0002, palette: ['black', 'red']};
Map.addLayer(current, vizParams, 'S5P N02');

UPDATE

I haven't been able to find any of the references to cloud within the NO2 data, however S-5P has a "cloud_fraction" band. I am able to set this to 5% and apply it as a mask to the NO2 dataset. Is this the correct way to get my result?

// Cloud mask for pixels with less than 5% cloud cover
var cloudMask = function(image) {
      return image.updateMask(image.select("cloud_fraction").lt(5));
    };

// Map the cloud masking functions over NO2 data
var collection = ee.ImageCollection('COPERNICUS/S5P/NRTI/L3_NO2')
  .select(['tropospheric_NO2_column_number_density','cloud_fraction'])
  .filterDate('2019-03-01', '2019-03-28')
  .filterBounds(geometry)
  .map(cloudMask);

// Reduce the collection to the mean of each pixel and display.
var meanImage = collection.reduce(ee.Reducer.mean());
var vizParams = {min: 0, max: 0.0002};
Map.addLayer(meanImage, vizParams, 'S5P N02 mean');


texts: Putting coded domains into list to be created in Python script?

texts:
I have a large list of coded domains.

They are numbers like:
01.6530.01, 05.7991.08, 05.7991.01, 05.7991.06

The code will be the exact same as the description.

I've tried to make a list using:

FRClist = [01.6530.**01**, 05.7991.08, 05.7991.01, 05.7991.06]

but it gave me a syntax error on the 01 that is bolded.

So I tried to change all the values to a string, using:

FRCstr = str(01.6530.**01**, 05.7991.08, 05.7991.01, 05.7991.06)

but it still gave me a syntax error on the same bolded 01.

In the end I hope to put the list into this:

arcpy.AddCodedValueToDomain_management(GDBpath, "FRC", FRClist[x], FRClistp[x] )

I have quite a few domains which require a lot of coded values.


texts: How can I be sure that line segments are joined at the same node?

texts:
I have received a single ArcGIS feature class which contains polyline segments. The segments appear to be joined together at their endpoints, but upon closer inspection they are not perfectly docked to the next polyline. I want to merge the segments according to a specific attribute, but need to maintain valid geometries without creating multipart features. How can I test whether the the vertices of 2 lines actually share the same coordinate and if not how can I force them to snap together if within a certain threshold (e.g 1 meter)?
Software: ArcGIS Advanced or SAFE FME


texts: How to get z coordinate of a point using pyqgis

texts: When I edit a vector layer in qgis, I select a lineString using Node Tool and I can see z coordinate of points in "Vertex Editor". In the image,we can see x,y,z coordinates of a list of points. But how can I get z coordinate of a point using pyqgis. In pyqgis, QgsPoint is a two demension point without z coordinate. And QgsPointV2 is a three demension point. When I get features from a VectorLayer and then get geometry from feature, It seems the geometry can only be convert to a two demension geometry (like QgsPoint, QgsPolyline, QgsPolygon). So how to get z coordinate of point.


texts: Angle at intersection point from two Lines

texts:
I would like to find the angle of two lines at their intersection point.
In my opinion it is the azimuth I'm looking for, right?

So I'm using Python with shapely for finding the intersection and QGIS for visualisation, but I'm not able to calculate the angle with Python.... So I decided to look how QGIS is doing this

  double QgsPoint::azimuth( const QgsPoint& other ) const
   {
     double dx = other.x() - m_x;
     double dy = other.y() - m_y;
     return ( atan2( dx, dy ) * 180.0 / M_PI );
   }

QGis

But i'm not able to reproduce this result

>>> x2=47.052246
>>> x1=47.052257
>>> y2=8.296885
>>> y1=8.296906
>>> xAns=x1-x2
>>> yAns=y1-y2
>>> np.arctan2(xAns,yAns)*180/np.pi
27.645975364825276

With QGis i'm able to get the angle graphically from Point2 over the intersection too Point1

Furthermore I tried

α = tan-1(cos(lat1) (long2 − long1) / (lat2 − lat1))

from ->charlespetzold.com/etc/AvenuesOfManhattan/

without success....

#

Ok after the help of @vagvaf and @vince
i wrote the following (not final)

from shapely.geometry import Point
import numpy as np

P1 = Point(47.052257,8.296906)
P2 = Point(47.052246,8.296885)
interP = Point(47.052150,8.296925)

dx = P1.x-interP.x
dy = P1.y-interP.y

dx2 = interP.x-P2.x
dy2 = interP.y-P2.y

azimuth1 = np.arctan2(dx,dy)*180/np.pi
azimuth2 = np.arctan2(dx2,dy2)*180/np.pi

print 'dx ' + str(dx) + '\n' + 'dy ' + str(dy) + '\n' + 'azimuth P1->interP ' + str(azimuth1)
azimuth11 = 360-azimuth1    # dx>0 & dy<0
print '360 - [azimuth1 P1->interP] ' + str(azimuth11)
print 'dx2 ' + str(dx2) + '\n' + 'dy2 ' + str(dy2) + '\n' + 'azimuth interP->P2 ' + str(azimuth2)
azimuth22 = 180-azimuth2    # dx<0 & dy>0
print '180 - [azimuth2 interP->P2] ' + str(azimuth22)
print 'azimuth11 - azimuth22 ' + str(azimuth11-azimuth22)

that results in:

dx 0.000107
dy -1.9e-05
azimuth P1->interP 100.069062699
360 - [azimuth1 P1->interP] 259.930937301
dx2 -9.59999999992e-05
dy2 4.00000000003e-05
azimuth interP->P2 -67.3801350517
180 - [azimuth2 interP->P2] 247.380135052
azimuth11 - azimuth22 = 12.5508022494

that is more or less equal to the qgis result


texts: How to convert UTM (easting,northing) to Geographic (lat,long) in degree decimal using geotools

texts: I want to convert easting, northing values in meters of UTM to lat,long of geographic coordinate system using geotools


texts: How can I find a county fips given a yahoo woeid

texts:
I've used the Yahoo API to pull some woeids for counties.  I need to map those woeids back to fips codes in order to associate them back to other data in my database, but I can't seems to figure out how to get a fips given a woeid.

Yahoo's concordance API handles the other direction.  With that, I can get a woeid for a given fips code, but how do I go the other direction?

Thanks!


texts: How to create "AGREE Stream" for DEM Reconditioning in Arc Hydro tool?

texts:
I want to deliniate streams and watershed using SRTM DEM. for this I am using Arc Hydro tool. While performing  DEM Reconditioning for DEM manipulation it asks AGREE Stream as a input.

How to create this so caled "AGREE Stream" to feed for DEM Reconditioning?


texts:
I have a file with over 27,000 addresses, and I'm trying to geocode it using the OpenStreetMap/Nominatim web service via the MMQGIS plugin on QGIS 2.10.

I keep getting an error after about 30 minutes of this plugin running in which 0 addresses get geocoded. Presumably, it is because my .csv file is not clean.

An error has occured while executing Python code:

Traceback (most recent call last):
  File "/Users/derekkaden/.qgis2/python/plugins/mmqgis/mmqgis_dialogs.py", line 1198, in run
    shapefilename, notfoundfile, fields, service, 1)
  File "/Users/derekkaden/.qgis2/python/plugins/mmqgis/mmqgis_library.py", line 1901, in mmqgis_geocode_web_service
    attributes.append(unicode(row[z], 'utf-8').strip())
UnicodeDecodeError: 'utf8' codec can't decode byte 0xc9 in position 14: invalid continuation byte

Is anyone familiar with this error?


texts: Creating precipitation raster map in GRASS GIS?

texts:
I am trying to create a potential erosion map of my basin for which I require a rainfall distribution map. I have 3 stations' rainfall data for the year 2007- 2015.

How do I create a rainfall map in GRASS GIS?


texts: GeoServer v2.18.0 Upgrade Web UI version of JQuery

texts:
I am looking into how I would update the version of JQuery used in the Web UI part of GeoServer (v2.18.0). The current version (JQuery v1.12.4) has 4 medium vulnerabilities, so I am looking to replace that with v3.5.1
Any GeoServer Developers around? I am a developer myself, but I have very little experience with GeoServer.
Update: This isn't off-topic, it is a library that GeoServer include and it needs updating. It isn't directly affecting GeoServer, but JQuery is backwards compatible. So there must be a way to update it?


texts: How to obtain percentiles of classes in classified image in Google Earth Engine?

texts:
https://code.earthengine.google.com/e430e2a79aecb65e3efbce23779a0abc

I have written this code for supervised classification using the cart classifier.
Now I want to extract the no. of pixels in each class i.e. water, urban, vegetation etc.

var rgb_vis = {min: 0, max: 0.3, bands: ['B4', 'B3', 'B2']};
var image =  ee.ImageCollection('LANDSAT/LT05/C01/T1_TOA')
  .filterBounds(ganga)
  .filterDate('2009-10-01', '2009-10-31')
  .median();
var region = image.clipToCollection(ganga);
// Map.addLayer(region, rgb_vis,'map');
//Data
// print(region);
var fc = water.merge(vegetation).merge(urban).merge(fallow);
var bands = ['B2','B3','B4','B5','B6','B7'];
//Train
var training = region.select(bands).sampleRegions({
  collection:fc,properties: ['landcover'],scale : 30});
var classifier = ee.Classifier.cart().train({features: training, classProperty: 'landcover', inputProperties: bands});
//Classification
var classified = region.select(bands).classify(classifier);
// print(classified);
//Display
Map.centerObject(fc, 11);
// Map.addLayer(region, rgb_vis, 'Landsat image');
// Map.addLayer(classified, {min: 0, max: 3, palette: ['0FD1F8','0E8205', 'EAF10B','EE0BF1']}, 'classification');
var output = classified.visualize({min:0,max:3,palette: ['0FD1F8','0E8205', 'EAF10B','EE0BF1'],forceRgbOutput:true});
Export.image.toDrive({image:output, region:geometry,description:'2009', scale: 30,crs:'EPSG:4326'});
Map.addLayer(output)


texts: What is name of vertices order property?

texts:
I'm trying to find the next vertices on polyline using python script in order to calculate the angle between them.

What is the name of vertices order property ?

Vertices have an attribute for X,Y and for his number- at polyline the number  is also the sequence of the vertices.

How can I use this property ?


texts: Using geometry points to get data in order to plot heatmap?

texts:
Excuse me if the title is confusing, English is not my main language, I'll try to explain myself as best as possible, also, I'm not great graphing, I'm learning.
My code plots points on a shapefile representation using Point from geometry and geopandas. I'd like to create a heatmap with the data I have. Some research has showed me that in order to do so I'd need to have a value that would represent the size of the ''heat'', could I create such data column using the points I already have, so, for example, if there's a lot of points at xy area at the plot, then the value for that point would represent small and bright colored heat?
Here's my code:
import pandas as pd
import matplotlib.pyplot as plt
import descartes
import geopandas as gpd
from shapely.geometry import Point, Polygon

street_map = gpd.read_file('shape_files/Street_Network_Database_(SND).shp')

fig, ax = plt.subplots(figsize=(15,15))

street_map.plot(ax=ax, color="grey", alpha=0.3)

df = pd.read_csv('data.csv')
crs = {"init":"epsg:4326"}

geometry = [Point(xy) for xy in zip(df["longitude"], df["latitude"])]

geo_df = gpd.GeoDataFrame(df, crs=crs, geometry=geometry)

geo_df.plot(ax = ax, markersize = 10, color = "blue", marker = "o", label = "POI")

Here's the results of this code:


texts: Obtaining XY coordinate when clicking on globe using ArcGlobe

texts:
Work in ArcGlobe and I want to calculate the coordinates of the globe latitude and longitude, but the coordinates from the mouse click is read from the screen, not the globe, because when I clicked on the corner of the screen it shows me the coordinates 0,0. How can I make it so that coordinates read with globe and not with the screen?
import arcpy
import pythonaddins

class ToolClass2(object):
    """Implementation for coord_addin.tool (Tool)"""
    def __init__(self):
        self.enabled = True
        self.shape = "NONE" # Can set to "Line", "Circle" or "Rectangle" for interactive shape drawing and to activate the onLine/Polygon/Circle event sinks.
    def onMouseDown(self, x, y, button, shift):
        pythonaddins.MessageBox("Long" + " " + str(x) + '\n'+ "Lat"+ " " + str(y), 'Coordinates', 0)
        pass
    def onMouseDownMap(self, x, y, button, shift):
        #pythonaddins.MessageBox("Long" + " " + str(x) + '\n'+ "Lat"+ " " + str(y), 'Coordinates', 0)
        pass
    def onMouseUp(self, x, y, button, shift):
        pass
    def onMouseUpMap(self, x, y, button, shift):
        pass
    def onMouseMove(self, x, y, button, shift):
        pass
    def onMouseMoveMap(self, x, y, button, shift):
        pass
    def onDblClick(self):
        pass
    def onKeyDown(self, keycode, shift):
        pass
    def onKeyUp(self, keycode, shift):
        pass
    def deactivate(self):
        pass
    def onCircle(self, circle_geometry):
        pass
    def onLine(self, line_geometry):
        pass
    def onRectangle(self, rectangle_geometry):
        pass


texts:
I have a question about GeoServer comportment when sending a GetFeature request with a DWithin filter.

<wfs:GetFeature version="2.0.0" service="WFS" outputFormat="JSON"  count="10"
xmlns="http://www.opengis.net/wfs/2.0"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xmlns:wfs="http://www.opengis.net/wfs/2.0"
xmlns:xlink="http://www.w3.org/1999/xlink"
xsi:schemaLocation="http://www.opengis.net/wfs/2.0 http://schemas.opengis.net/wfs/2.0/wfs.xsd">

<wfs:Query typeNames="myFeature">
    <wfs:PropertyName>label</wfs:PropertyName>
    <fes:Filter xmlns:fes="http://www.opengis.net/fes/2.0" xmlns:gml="http://www.opengis.net/gml/3.2">
        <fes:And>
            <fes:PropertyIsEqualTo matchAction="ANY" >
                <fes:ValueReference>label</fes:ValueReference>
                <fes:Literal>XXXX</fes:Literal>
            </fes:PropertyIsEqualTo>
            <fes:DWithin>
                <fes:ValueReference>geom</fes:ValueReference>
                <gml:Point srsDimension="2">
                    <gml:pos>43.55749 1.525864</gml:pos>
                </gml:Point>
                <fes:Distance>100</fes:Distance>
            </fes:DWithin>
        </fes:And>
    </fes:Filter>
</wfs:Query>

Considering this defect : https://osgeo-org.atlassian.net/browse/GEOS-937, I suppose the unit uses for the distance by GeoServer is the degree (unit defined by the coordinate system of data: EPSG:4326).

But in my PostGIS log, I can see this request:

ST_DWithin("geom",ST_GeomFromText('POINT (1.525864 43.55749)', 4326),0.001035624608029343))

Why does GeoServer apply a transformation on my distance value?


texts: Is it possible to query features based on attributes from the client side of qgis2web?  How can this be achieved?  I have a single data set that I want to serve over the web, while allowing users to search for locations that have specific attribute values.


texts:
When processing a dtm from lidR::grid_terrain() I'm seeing several warnings that say:

  There were 4 degenerated ground points. Some X Y Z coordinates were
  repeated. They were removed.

What is a degenerated ground point?

I am wondering if this this normal, or indicative of some larger problem with my .las tiles.


texts: How to remove "ERROR 1" messages from a standalone QGIS script?

texts:
I am running a standalone script for QGIS which creates a high resolution grid and clips shapefiles onto it. Works perfect if the cell width and height have a high values (in my case, >200m for both). Anything lower than this and I suddenly get the following error:

  ERROR 1: Attempt to write non-polygon (POINT) geometry to POLYGON type shapefile.

I have read this post with a similar problem, and as mentioned in the answer, the results I obtain are correct so I am not worried about that.

Is there a method to remove this error message(s) from appearing? It looks incredibly messy:

I've also included a stripped-down snippet of code I use:

def run():

    Cellsize = 10
    layerPath = "path\to\shapefile"
    extent = QgsVectorLayer( layerPath, '', 'ogr' ).extent()

    centerx = (extent.xMinimum() + extent.xMaximum()) / 2
    centery = (extent.yMinimum() + extent.yMaximum()) / 2
    width = extent.xMaximum() - extent.xMinimum()
    height = extent.yMaximum() - extent.yMinimum()

    outputs_0=general.runalg("qgis:creategrid", Cellsize, Cellsize, width, height, centerx, centery, 1, 'EPSG:7405', None)
    outputs_1=general.runalg("qgis:clip", outputs_0['SAVENAME'], "path\to\shapefile", "create\path\for\result")

try:
    run()
except Exception:
    pass


texts: Trying to convert coordinates from WGS84 Web Mercator Auxiliary Sphere to WGS84

texts:
I've been struggling with this for way too long now.  I've got some data that was georeferenced in the WGS84 Web Mercator Auxiliary Sphere that should have been done just in WGS84.

The Auxiliary Sphere (ESPG: 3857) is the projection commonly used by Google, Bing, and recently ArcGIS Online.  The projected coordinates I have are in meters, as per Aux Sphere specs, but I've got to convert them to degrees and bring them into WGS84.

I know the Aux Sphere is just a projected version of WGS84, but I can't seem to wrap my mind around converting coordinates from a projected coordinate system to a geographic coordinate system.

I have some code that will convert from WGS84 Web Mercator (WKID: 102113) to WGS84, but because that system uses a spherical globe, rather than the ellipsoid used by WGS84 and the Aux Sphere, my results are a bit off.  I've tried modifying it, but I can't seem to get it just right.

Here's my code if that is helpful:

public static void ToGeographic(ref double[] mercator)
    {
        if (Math.Abs(mercator[4]) < 180 && Math.Abs(mercator[5]) < 90)
            return;

        if ((Math.Abs(mercator[4]) > 20037508.3427892) || (Math.Abs(mercator[5]) > 20037508.3427892))
            return;

        double x = mercator[4];
        double y = mercator[5];
        double num3 = x / 6378137.0;
        double num4 = num3 * 57.295779513082323;
        double num5 = Math.Floor((double)((num4 + 180.0) / 360.0));
        double num6 = num4 - (num5 * 360.0);
        double num7 = 1.5707963267948966 - (2.0 * Math.Atan(Math.Exp((-1.0 * y) / 6378137.0)));
        mercator[4] = num6;
        mercator[5] = num7 * 57.295779513082323;

    }


texts: Intersection/Addition of multiple tables in PostGIS

texts:
I want to intersect multiple table together. For each intersection geometry, I want to sum the value of all the layers intersecting at that place.

My current query does it by step, getting the intersection between two tables first and getting the new value by summing for each new geometry, and then finding the intersection of these two tables with a third table, etc. Eventually I can have up to 10 tables to combine and sum the value.

My current query to do so for 3 tables is as follow:

WITH intersect_ab AS (
SELECT ST_INTERSECTION(a.geometry, b.geometry) as geometry, a.val + b.val as val
FROM mytable_a a
INNER JOIN mytable_b b ON ST_INTERSECTS(a.geometry, b.geometry)
),
     intersect_ab_c AS (
SELECT ST_INTERSECTION(ab.geometry, c.geometry) as geometry, ab.val + c.val as val
FROM intersect_ab ab
INNER JOIN mytable_c c ON ST_INTERSECTS(ab.geometry, c.geometry)
)
SELECT geometry, val
FROM intersect_ab_c;

Here is the query plan associated with this query:

              QUERY PLAN
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
 CTE Scan on intersect_ab_c  (cost=2359.00..2371.08 rows=604 width=40) (actual time=0.735..707.197 rows=2201 loops=1)
   CTE intersect_ab
     ->  Nested Loop  (cost=0.14..179.04 rows=4857 width=40) (actual time=0.611..156.574 rows=1251 loops=1)
           ->  Seq Scan on mytable_b b  (cost=0.00..1.32 rows=32 width=211) (actual time=0.003..0.010 rows=32 loops=1)
           ->  Index Scan using mytable_a_gix on mytable_a a  (cost=0.14..4.79 rows=1 width=217) (actual time=0.115..0.400 rows=39 loops=32)
                 Index Cond: (geometry && b.geometry)
                 Filter: _st_intersects(geometry, b.geometry)
                 Rows Removed by Filter: 1
   CTE intersect_ab_c
     ->  Nested Loop  (cost=0.14..2179.96 rows=604 width=40) (actual time=0.734..705.353 rows=2201 loops=1)
           ->  CTE Scan on intersect_ab ab  (cost=0.00..97.14 rows=4857 width=40) (actual time=0.612..157.457 rows=1251 loops=1)
           ->  Index Scan using mytable_c_gix on mytable_c c  (cost=0.14..0.42 rows=1 width=135) (actual time=0.022..0.032 rows=2 loops=1251)
                 Index Cond: (ab.geometry && geometry)
                 Filter: _st_intersects(ab.geometry, geometry)
                 Rows Removed by Filter: 0
 Planning Time: 12.183 ms
 Execution Time: 707.623 ms

I have a geometry index on all tables (mytable_a, mytable_b and mytable_c). mytable_a has 1027 polygons, mytable_b has 32 polygons and mytable_c has 373 polygons.

Is there a way to make this more efficient? It seems like there ought to be for a "simple"/"likely common" spatial join like this one, but I'm not sure what the best approach would be.


texts:
I use ArcGIS 10.1 on a Windows XP. ArcGIS completely shuts down when I attempt to launch the "Raster Calculator" tool. Has this happened to someone else? Any ideas on how to solve this annoying issue?
Otherwise, can I calculate the sum of raster values using command-line GDal?


texts: How to create an interactive map on a web page?

texts: I am working on a website using Django. A central component of this website will be a map page with accompanying functionality. I was planning on using Javascript to render the map and use AJAX to get data from my server for the usage of the map. I am not sure at the best way to implement the map. I know there a number of Web GIS softwares that have nice prebuilt functionality. The problem is that I want to use my own maps (created from QGIS) and data. I am not sure how to configure them to use my shape files for the map. The other option is to use HTML Canvas. I could easily get something to work but I would rather have the functionality of an existing Web GIS software to use. Any suggestions?


texts:
I have a question when the projected EPSG::4326 and EPSG::900913.
OpenLayers uses by default EPSG::4326 using the projection coordinates in decimal format already uses Google Maps with coordinates EPSG::900913 meters.

I saw several examples of creating the map in the projection and display property EPSG::900913 projection in EPSG::4326.

What difference does it make on the map? If I want to record a point in the database will be in the map projection EPSG::900913 in this case.


texts:
We have point data (x,y in m from a local 0,0) from an archaeological campaign, that we would like to display embedded into real-world cartography.

We know either the real world coordinates of 0,0 (e.g. 337861.49,5510278.69) and a rotation angle 48,5° or the coordinates of a second point (x1,y1).

Version A

One obvious idea would be to use an existing projection in PROJ like e.g. oblique mercator projection (like PROJ omerc) for this, however we end up with a large error (~200 m) and no rotation whatsoever.

Using PROJ tpqed we have a rotation but again a non-offsettable shift in E and N and also mismatched scaling.

Version B

Another option would be starting from a fitting cartesian, projected CRS like EPSG:25832 and just perform an affine transformation. This should be doable via SC_DerivedCRS in a WKT file, however I couldn't find any usable instructions how to do this with WKT2.

How do I wrangle the errors in A or the syntax for the WKT file?


texts: Get  features inside a rectangle after drawing it

texts:
Let's start with https://jsfiddle.net/ziyuang/dfekwzhd/

Say, now I want to allow users to specify areas by drawing rectangles. After the map having a rectangular area, I want to have the information, like properties (i.e., properties in line 85 of the snippet) and feature IDs, of the circles inside the area.

The drawing part seems fine, for which I may use L.rectangle or Leaflet.draw. But I am not sure about the querying part.

Any suggestions?


texts: Connect multiple WMS layers in CARTO querying to GetFeatureInfo. How do I get GetFeatureInfo information from the attached WMS layer?


texts:
I would like to dissolve only polygons that are adjacent to one another.  Ultimately I would like to start with a raster representing landscape patches and delineate the separate patches.  Here is my example, drawing from the raster::rasterToPolygons function in R:

r <- raster(nrow=18, ncol=36)
set.seed(999)
r[] <- runif(ncell(r))
r <- r > 0.7
pol <- rasterToPolygons(r, fun=function(x){x==1})
pol.union <- rasterToPolygons(r, fun=function(x){x==1}, dissolve=T)

plot(r)
plot(pol, add=T)

plot(r)
plot(pol.union, add=T)

Without using the dissolve=T argument in rasterToPolygons, it returns a single polygon for each pixel where r==1 (189 polygons).  If I use dissolve=T it dissolves all of these into a single polygon.  What I want is something in between where each non-adjacent patch is it's own polygon (it would look like pol.union when plotted but would contain more than one polygon).  I believe what I want is to take pol and dissolve only the adjacent polygons.

I realize that this may require some input as to what is considered "adjacent" (i.e., whether only polygons sharing a boundary are considered adjacent or if polygons sharing a single point are considered adjacent).


texts: Selecting features "above" or "below" a line using R

texts:
Given a line and a set of points, I can't figure out how to use sf to identify which side of the line each point falls.

A small reproducible example follows, adapted from a different question

# Load Libraries ----------------------------------------------------------

library('sf')

# Test data ---------------------------------------------------------------

points.df <- data.frame(
    'x' = c(-53.50000, -54.15489, -54.48560, -52.00000, -52.57810, -49.22097, -48.00000),
    'y' = c(-38.54859, -41.00000, -38.80000, -38.49485, -38.00000, -40.50000, -37.74859)
)

line.df <- data.frame(
    'x' = c(-54.53557, -52.00000, -50.00000, -48.00000, -46.40190),
    'y' = c(-39.00000, -38.60742, -38.08149, -38.82503, -37.00000)
)

# Create 'sf' objects -----------------------------------------------------

points.sf <- st_as_sf(points.df, coords = c("x", "y"))

st_crs(points.sf) <- st_crs(4326) # assign crs

line.sf <- st_sf(id = 'L1', st_sfc(st_linestring(as.matrix(line.df), dim = "XY")))
st_crs(line.sf) <- st_crs(4326) # assign crs

# Plots -------------------------------------------------------------------

xmin <- min(st_bbox(points.sf)[1], st_bbox(line.sf)[1])
ymin <- min(st_bbox(points.sf)[2], st_bbox(line.sf)[2])
xmax <- max(st_bbox(points.sf)[3], st_bbox(line.sf)[3])
ymax <- max(st_bbox(points.sf)[4], st_bbox(line.sf)[4])

plot(points.sf, pch = 19, xlab = "Longitude", ylab = "Latitude",
     xlim = c(xmin,xmax), ylim = c(ymin,ymax), graticule = st_crs(4326), axes = TRUE)

plot(line.sf, col = "#C72259", add = TRUE)
text(st_coordinates(points.sf), as.character(1:7), pos = 3)

In this example, it's easy to verify that points 2 and 6 fall south of the line, and the rest north. How can I automate the labeling?

Non sf based answers are welcome too.


texts: QGIS - Extract by location, what does Cross do?

texts:
I don't understand what Cross does when extracting by location, and how it differs from Overlap- the definition on the user guide is very confusing - it says:

  Cross: Returns 1 (true) if the supplied geometries have some, but not all, interior points in common and the actual crossing is of a
  lower dimension than the highest supplied geometry. For example, a
  line crossing a polygon will cross as a line (selected). Two lines
  crossing will cross as a point (selected). Two polygons cross as a
  polygon (not selected).

Can someone explain this in a clearer way, maybe with a simple example?


texts: How to create a dataframe from an RDD using Spark Core?

texts:
I am new using Spark Core and I need to read a csv file, delete the header and after that convert to RDD and to dataframe.

I tried this, but I cant get a dataframe:

from pyspark import SparkContext, SparkConf

conf = SparkConf().setAppName("Pruebas").setMaster("local")
sc = SparkContext(conf=conf)

file = sc.textFile("hdfs://localhost:9000/PECE.csv")

def f(idx, iter):
    output = []
    for sublist in iter:
        output.append(sublist)
    if idx > 0:
        return (output)
    else:
        return (output[1:])

def p(x): print(x)

#Utilizamos la funcion para quitar la cabecera
data = file.map(lambda x : x.replace("\"","").split(";")).mapPartitionsWithIndex(f)

print(data)

file2 = sc.textFile("hdfs://localhost:9000/PECM.csv")

def f2(idx2, iter2):
    output2 = []
    for sublist2 in iter2:
        output2.append(sublist2)
    if idx2 > 0:
        return (output2)
    else:
        return (output2[1:])

def p(x2): print(x2)

data2 = file2.map(lambda x : x2.replace("\"","").split(";")).mapPartitionsWithIndex(f2)

print(data2)

val dataFrame = spark.read.csv("PECM.csv")