texts: Display angle of each line feature within a specific field for that layer using UpdateCursor method

texts:
So I want to calculate angles for each line feature layer using UpdateCursor method. So far I have progressed with the code but what it is doing is; it is giving the result of the last feature angle and substituting it for all the fields. I tried to avoid that function calling and now I am getting this error.
I have a polyline layer in a GDB consisting of smaller line features, I want to calculate the angle as in orientation of each line feature and update them in a marked field above in the attribute table.
Since I am new to this, I don't understand where I am going wrong.
Is there some other concrete method to do this?
Here is my code:
count1 = 0
spatial_Reference = arcpy.Describe(output_split_merge_fibre_drop_sheath2).shapeFieldName
with arcpy.da.UpdateCursor(output_split_merge_fibre_drop_sheath2, ["Fibre_SH_D", "SHAPE@", "OBJECTID"]) as cursor:
            for row in cursor:
                shape = row[1]
                count1 = count1+1
                feature = row.append(shape)
                x1 = feature.positionAlongLine(0.51, True).firstPoint.X
                x2 = feature.positionAlongLine(0.49, True).firstPoint.X
                y1 = feature.positionAlongLine(0.51, True).firstPoint.Y
                y2 = feature.positionAlongLine(0.49, True).firstPoint.Y
                radian = math.atan2((y1 - y2), (x1 - x2))
                row[0] = str(radian)
                print(row[0])
                cursor.updateRow(row)

The error message I am getting is:-

AttributeError: 'NoneType' object has no attribute 'positionAlongLine'


texts: How identify Layer name when removing it from legendInterface qgis?

texts:
I am trying to connect the removing layer action with a slot for doing some processing after that.
I already try this:

QObject.connect(self.iface.legendInterface(), SIGNAL("itemRemoved()"),self.RemoveScenario)

but i need to identify the removed layer name.
anyone have an idea how to do that?


texts: How can I convert the Gaussian coordinates( Gauß-Krüger) to the WGS84 coordinates


texts:
I'm trying to create an interpolation of xyz data in QGIS. The data has successfully interpolated, but there is some strange striping or banding on the final result that I cannot seem to get rid of. See the images here.

The first (black and white) was created by saving the data points as a shapefile, then interpolating. Parameters in TIN interpolation were cubic method, with x and y set to 5 meters. I let col and row define themselves- how would that affect the result?

The second was with the xyz data loaded as a deliminated text layer. Interpolated the same as above.

The methods with the basic interpolation package do not change the result. Banding appeared with both linear and cubic options. I also looked at contours which exhibit the same problem (see the third image). Data itself is fine.

The banding might have a geographic cause- it is occuring over an area deeper than the surroundings. However, I do not expect it to create this level of an issue.


texts:
I'm using a device with a chipset "ublox NEO-M8U" and an external antenna. I tried to write a parser for NMEA data that I received from this device, but I cannot find an answer to the question: why I don't receive milliseconds in timestamps (GNGGA/GNRMC) only always the regex looks like that HHMMSS.00.
Example of sentences:
$GNRMC,152418.00,A,6024.02224,N,00518.92204,E,0.008,,181120,1.38,E,D*35
$GNGGA,152417.00,6024.02213,N,00518.92186,E,2,12,0.85,19.2,M,43.8,M,,0000*74

What kind of pieces of information says me in the datasheet that I will receive (true) milliseconds in Timestamp?


texts: What are some good resources related to Geodjango?

texts:
I am trying to learn Geodjango. Has anyone used it in their own project?

The documentation on the official site is sparse. Can someone point out the best

resources helpful in learning it other than official documentation ?

Please also tell me how is your experience with it and other web application frameworks that you have used for mapping.


texts: How to rotate external graphic using mouse in Openlayers?

texts:
I would like to rotate externalGraphic using mouse something like this URL.

In this URL polygon is used to rotate it but I want an image to be rotated using external graphics tag.

Any help is appreciated.


texts: How do I sort a geoJSON feature collection alphabetically by a property value?

texts:
I am trying to sort the objects in a geoJSON FeatureCollection alphabetically by the values of one of their properties. How do I get to the nested object properties with JavaScript? I have tried using .sort to no avail.

Sample geoJSON to sort on PARK_NAME property:

"type" : "FeatureCollection",

"features" : [{
    "type" : "Feature",
    "id" : 1,
    "properties" : {
        "CRAG_NAME" : "Alberton Road Rock",
        "PARK_NAME" : "Patapsco State Park",
        "LAT_DD" : 39.3168,
        "LONG_DD" : -76.8033,
        "GUIDE_URL" : "http://indy-adventures.net/climbing/climb-maryland/guide.html#Alberton",
        "OTHER_URL" : "http://www.dnr.state.md.us/publiclands/central/patapsco.asp"
    },
    "geometry" : {
        "type" : "Point",
        "coordinates" : [-76.8033, 39.3168]
    }
}, {
    "type" : "Feature",
    "id" : 2,
    "properties" : {
        "CRAG_NAME" : "Annapolis Rocks",
        "PARK_NAME" : "Appalachian Trail",
        "LAT_DD" : 39.5583,
        "LONG_DD" : -77.5987,
        "GUIDE_URL" : "http://indy-adventures.net/climbing/guide.html#AnnapolisRocks",
        "OTHER_URL" : "http://www.dnr.state.md.us/publiclands/at.asp"
    },
    "geometry" : {
        "type" : "Point",
        "coordinates" : [-77.5987, 39.5583]
    }
}]
}

.....


texts: Labelling a multipoygon scattered in different areas

texts: I am using mapnik and tilecache to generate tiles and openlayers to view map in the browser. In one of my shapefiles(multipolygon) with administrative boundaries, I found that more than one records have their polygons displayed in two areas in the map. In the map, label is placed in center position of the two polygons. When I opened the shapefile in Qgis, I found that label is displayed for each of the two polygons correctly. How will this type of labelling be done in mapnik. Thanks in advance.


texts: How to calculate distance to a feature with gdal_proximity?

texts:
I'm using gdal_proximity to find the distance to the nearest major river across the USA (lower 48 states). I've projected the NHD+ network flowlines to Conus Albers (epsg:5070), selected rivers with stream order > 5, and rasterized, burning rivers as 255, no river as 0. This much is fine, but now I need to find distance to nearest river for sites within 50km. The input file is at 30m resolution on continental scale, so is very large, but the conversion should be a simple gdal_proximity command:

gdal_proximity.bat -values 255 -distunits GEO -maxdist 50000 -nodata -999 infile.tif outfile.tif -co COMPRESS=DEFLATE -co BIGTIFF=YES -co TILED=YES

This seems to -almost- work, but is producing a weird geometric pattern in the output (see image). The data that are present in the output have processed correctly. Can anyone suggest why so much of the output is missing?

Edit:
To test if this was caused by any of the optional parameters, I ran gdal_proximity again in this configuration:

gdal_proximity.bat H:\data\tmp\NHDplus_network_flowline_SO6plus.tif H:/data/tmp/NHDplus_network_flowline_SO6plus_proximity.tif -values 255 -maxdist 50000 -of GTiff

Which yielded essentially the same result:

My only thought is that it may be related to the size of the raster (~100 gb uncompressed. As far as I'm aware there isn't a limit to the size of a BigTiff, but maybe there's a limit to what gdal can analyse effectively?


texts: I have a vector map (consisting of several polygons) and I would like to create several overlapping tiles. I did the same for a raster map using "split raster". For the vector map, I created an overlapping fishnet (using "buffer" and "feature envelope to polygon") and I tried to use it as split features, but I don't get the correct result (i.e.the overlap is missing).


texts: GDAL - Polygonize multiband raster and save pixel values of each band to corresponding attribute fields

texts:
I'm trying to polygonize 2 band raster(which was previously reclassified) and put values of each band into seperate attribute table column. Both bands overlap each other perfectly.

This piece of code works perfectly for 1 band rasters:

shapeFile = "shape1.shp"
driver = ogr.GetDriverByName("ESRI Shapefile")
outDS = driver.CreateDataSource(shapeFile)
outLayer = outDS.CreateLayer("lyr1", srs)
f1 = ogr.FieldDefn('DN', ogr.OFTInteger)
f2 = ogr.FieldDefn('DIF', ogr.OFTInteger)
outLayer.CreateField(f1)
outLayer.CreateField(f2)
gdal.Polygonize(raster.dataset.GetRasterBand(1), rasterMask, outLayer, 0, [], callback=None )
#gdal.Polygonize(raster.dataset.GetRasterBand(2), rasterMask, outLayer, 1, [], callback=None )

However the commented piece of code where I try to put pixel values of 2nd band into the 'DIF' field doesn't work (and it's field index is 1, I've checked that).

The Polygnize function doesn't seem to have an option for multiple band rasters and I'm not sure if calling it 2 times in a row is even a good idea.

I was thinking of polygonizing each band to a different shapefile and then somehow merging those shapefiles to match each field value but I'm not even sure if that's possible.

How should I proceed?


texts: Opening MODIS Raster MOD13Q1 .hdf file comes in as Float32 instead of 16-bit signed integer

texts:
Trying to open a MODIS 16 day composite NDVI raster from the .hdf file.  After selecting just the NDVI band, QGIS nicely displays the image.  However, on this image the range of values is between -3.57096e+06 to 7.97734e+07.  So outside the expected range of -3000 (fill value) and 10000.  Looking at the Metadata, says the Data Type is Float32 instead of the expected 16-bit signed integer (Int16?).

Is there a way to choose the Data Type of the raster being imported or to change it once imported?


texts: dateline issue: wrapDateLine with MapFile

texts:
I am using wrapDateLine properties for OpenLayers.Layer.WMS for handling dateline issues. The problem with wrapDateLine is that I am getting the same polygon repetitively around the wrapped map as a user selects a desired polygon around dateline. Is there any way to select the polygon around dateline in a single wrapped map and avoiding repetitions? My OpenLayers.Layer.WMS constructor is as:

var I = new OpenLayers.Layer.MapServer("region", "http://localhost/cgi-bin/mapserv?map=/home/public_html/test1.map",{poly='Bermuda', layers: 'myregion', format:'image/png'},
    {wrapDateLine: true}
    );

e.g. Bermuda in my example is in dateline and if I use wrapDateLine:true, I get complete polygon of Bermuda but it is repeated (wrapped). I want complete polygon of Bermuda but no repetition. Is there any other way for it?

I see that there are options like OpenLayers.Bounds() and layer.adjustBounds(). I am not sure if it works in my problem.


texts: How to combine/merge attribute data from 1 layer to another, polygon to point, merge attributes and wrap/carriage return attributes

texts:
I am trying to merge field areas from a polygon layer into a point layer
So the field area from this polygon onto the dot layer.

Where the Point layer, has an String attribute, and set to Multiline, so it would look like this:

I can generated the area measurement, using Field Calculator QGIS 2.14, Expression $area
(Though I understand there can be/is errors with this method)
Point Layer:

Polygon layer:

Join by field location
You can do a spatial join but that just adds another column with the attributes, but I need another step.
Join by field location I believe adds the area location from the polygon layer to the point layer.  Go to Menu, Vector, Data Management Tools, Join Attributes by location
(Hopefully this will work when you have smaller polygons and points closer together as the distances between points/polygons is quite small)

I need to merge the columns generated as well.  So the area column needs to become more like the Label column. (see picture)
I am looking to avoid manually copying and pasting the data from one column into another in the attribute table.  In that case, it will be easier to measure each polygon individually using the identify tool, and then enter the label name and area one by one.

Try 3 - Pnz answer - Unsuccessful
Create 2 new layers
Names 2 - Point layer
Attribute: Label, String, multiline
Polygon 2 - polygon layer
Created area attribute on Polygon 2 using field calculator
Had both layers editable
Ran Vector -> Data Management Tools -> Join by Location

Created a new point layer but with 3 attribute columns:
Label, id, area

Test 2
Ran Vector -> Data Management Tools -> Intersect

This creates a new point layer with 3 attribute columns
Label, id, area

The issue is that the label needs to be as follows

I don't think "intersect" or "join by location" will join, and "wrap" in effect 2 lots of attribute column data.
I have had both layers editable, but since the Intersect and Join by location tool create new layers, I don't think the layers being editable is the issue.


texts:
I need to implement ISO 19115-1 standards to a collection of datasets using ArcCatalog. However, I'm having trouble locating documentation about actually using the updated standard in Arc.

I've found documentation for 10.5-- it states that previous 19115:2003 versions will be compatible with 19115-1, but does not provide documentation for a revised 19115-1 Metadata Style (This document is for 10.5, and I do not have access above 10.2).

Based on previous experience, I thought saving a new .cfg file in my C:/Arc/Metadata folder would be a solid workaround, but I can't locate the file anywhere.

Have any of you used ISO 19115-1, what do you suggest?


texts: Smooth mosaicking of smaller orthomosaics to form a homogenous mosaic

texts:
This image below is generated using Erdas Imagine where I used mosaicpro tool to create a super mosaic of orthomosaics generated using PIX4d

I used following options:

Color Balancing
Color Dodging
Histogram equalisation
Feathering in overlaping areas

Is there a way to smooth these results further such as to remove those strips of light green patches?


texts:
Im looking at making a standalone script, that on execution will run the QGis MapComposer using a previously generated template and output a PDF/PNG of the map. This process of making a map will be the last step in a operational system that runs twice daily.
I considered using the --snapshot function, but because I'd also like to have a number of images (logos etc) on the layout, the snapshot becomes unworkable.
In essense, what I'd like to do is: From a standalone script, save a PDF or PNG that is based on a qpt-template-file. For each run of the script, I'll have a new .tif file, which I can write into a qgs-project-file, in order to load the right file. With this, what I expect that I'm supposed to do is:

(1) initialize Qgis
(2) open up the qgs-file (in order to specify file, stretch etc)
(3) open up the map composer with the qpt-file (in order to have all my static images, labels etc)
(4) save the map to a PDF / PNG (format is of little concern)
(5) close Qgis

What I've been having trouble with, is steps 2 and 3. Steps 1, 4 and 5 I've grabbed from other scripts and various examples. The pyqgis cookbook doesn't really cover the parts that are troubling me, and API Documentation is significantly out of my league.
Is anyone here willing and able to give me a hand with my problems, either by pointing me at the right functions, or by providing a script doing something along the lines of what I wish to do?


texts: I would like to set up my own custom map by changing the images in openlayers 3. Is something like this possible and how would I go about doing it?


texts: How to determine and send alert SMS for next bus stop in a GPS tracking system?

texts:
In my gps tracking system (PHP/MySql) I have got 5 Stops for a vehicle.

Whenever the bus reaches stop 1, the passengers at stop 2 should get an SMS alert. And in the same way for the sops followed by.

I have given stops as points. How can I determine whether the bus reached a particular stop? Should always monitor all 6 stop whether the bus reached once the trip is started?
Can anybody give some idea?

(Map used: Google Map/Openlayers)

See the image.


texts: What is the easiest way to purge all installed ArcGIS Pro AddIns?

texts:
I ran the BuildAll.cmd in arcgis-pro-sdk-community-samples and ended up with:

I need to reverse this temporarily (or permanently and I'll rerun BuildAll.cmd).  What's the easiest way to get back to "factory default"?


texts:
I am trying to build a geocoding app on OSM data, that in regard to arbitrary location will return points within specified radius,  as well as distances to the points.

I have the following table:

id      parent_id   key       value         geometry
4821719 221252430   natural   water         POINT (-8.69802 40.60405)
3483776 3335087     leisure   sports_centre POINT (-8.71858 40.61685)
3483777 3335087     leisure   sports_centre POINT (-8.71608 40.61672)

all the geometries are points

I do it with the help of Python geoalchemy2 package with the following expression:

query = session.query(Table,func.ST_Distance_Sphere(Table.geometry, arbitrary_point))\
               .filter(func.ST_Distance_Sphere(Table.geometry, arbitrary_point) < radius)

which translates to sql as:

SELECT table.id, table.parent_id, table.key, table.value, ST_AsEWKB(table.geometry) AS
geometry,
ST_Distance_Sphere(table.geometry, ST_GeomFromWKB(:ST_GeomFromWKB_1, :ST_GeomFromWKB_2)) AS
"ST_Distance_Sphere_1"
FROM table
WHERE ST_Distance_Sphere(table.geometry, ST_GeomFromWKB(:ST_GeomFromWKB_3, :ST_GeomFromWKB_4)) < :ST_Distance_Sphere_2

the returned result is like:

id  parent_id   key            value          geometry                 ST_Distance_Sphere_1
2   3788923     natural        water          POINT (-0.47833 52.25222)  9906.627366
5   3865117     leisure       sports_centre   POINT (-0.50365 52.16449)  0.333864
9   4082699     addr:country    GB            POINT (-0.49252 52.16019)  896.977083
12  4082701     addr:country    GB            POINT (-0.49220 52.16020)  915.175355
16  4082703     addr:country    GB            POINT (-0.49205 52.16015)  927.182032

the query takes long time to execute, so tried to add index with

CREATE INDEX geo ON table USING GIST ( geometry );

but it did not help at all

so I tried to remove this index, having the drop INDEX never finish.

I need the ways to speed-up the query.

Also, as well quering individual locations, I will need query
large arrays of locations at a time.

I am entirely new to PostgreSQL.

_____edit_____

EXPLAIN in presence of the index returns the following:

[('Gather  (cost=1000.00..3863720.18 rows=2741158 width=68)',),
 ('  Workers Planned: 2',),
 ('  ->  Parallel Seq Scan on osm  (cost=0.00..3588604.38 rows=1142149 width=68)',),
 ("        Filter: (st_distance_sphere(geometry,
'0101000020E6100000BEDEFDF15E6D21C0EF74E789E74E4440'::geometry) <
'2000'::double precision)",)]

'EXPLAIN' having the index dropped:

[('Gather  (cost=1000.00..3863720.18 rows=2741158 width=68)',),
 ('  Workers Planned: 2',),
 ('  ->  Parallel Seq Scan on osm  (cost=0.00..3588604.38 rows=1142149
width=68)',),
 ("        Filter: (st_distance_sphere(geometry,
'0101000020E6100000BEDEFDF15E6D21C0EF74E789E74E4440'::geometry) <
'2000'::double precision)",)]


texts: Can Virtual Fields be used for a rendering style?

texts:
In QGIS 2.6, creating a virtual field for a layer and adding up 2 integer fields work fine. I get the correct results and I can use that field for styling.

However, when doing a Joins with 2 layers and using a virtual field to add up an integer column from each layer, the results displayed are correct. But QGIS does not recognise the values when styling.

When I try to classify this in the Categorized styling option, no data values are recognised.

The virtual field tooltip describes:

  "A virtual field will be recalculated every time it is used. Its definition will be saved in the project file. It will not be saved in the dataprovider and therefore its values not be available in other software."

As it won't be saved in the dataprovider, does that mean QGIS is not even able to read the data? Strangely, I am able to export the layer, including the virtual field, as a .csv which outputs all the correct results for the virtual field.

Is this a possible bug?

EDIT:

The bug report for this problem has been Closed, although I will wait until the next major version of QGIS to test this.

If anyone can confirm this (with screenshots preferably) then I will happily accept your answer.


texts:
Is there any way to autoinstall the packages another package depends on. I just downloaded RStudio 1.3.62 and can't remember having these issues before. They looked like this when I tried to install mapview:

ERROR: configuration failed for package ‘systemfonts’
* removing ‘/home/robin/R/x86_64-pc-linux-gnu-library/3.6/systemfonts’
Warning in install.packages :
  installation of package ‘systemfonts’ had non-zero exit status
ERROR: dependency ‘systemfonts’ is not available for package ‘gdtools’
* removing ‘/home/robin/R/x86_64-pc-linux-gnu-library/3.6/gdtools’
Warning in install.packages :
  installation of package ‘gdtools’ had non-zero exit status
ERROR: dependency ‘gdtools’ is not available for package ‘svglite’
* removing ‘/home/robin/R/x86_64-pc-linux-gnu-library/3.6/svglite’
Warning in install.packages :
  installation of package ‘svglite’ had non-zero exit status
ERROR: dependency ‘svglite’ is not available for package ‘leafpop’
* removing ‘/home/robin/R/x86_64-pc-linux-gnu-library/3.6/leafpop’
Warning in install.packages :
  installation of package ‘leafpop’ had non-zero exit status
ERROR: dependencies ‘leafpop’, ‘svglite’ are not available for package ‘mapview’
* removing ‘/home/robin/R/x86_64-pc-linux-gnu-library/3.6/mapview’
Warning in install.packages :
  installation of package ‘mapview’ had non-zero exit status

I could install them manually. But then they depend on other packages which I had to install again. Any idea how I can stop this recursive process?


texts: How to display QGIS "processing ..." bar from a plugin?

texts:
QGIS seems to have an integrated animation which purpose is to tell the user "Please wait, I am not freezing, I am processing"

What is especially cool about it is that it is perfectly integrated to QGIS GUI (docked in the bottom) and doesn't fill itself like a QProgressBar : it just bounces from one side to another, giving no informations about the progression, just saying "Wait, I am processing".

I would like to know if it is possible to call this bar from a plugin.

I know of alternate options, which are displaying a QProgressBar (in my plugin gui or in a QMessageBar).
I even found this post, which seemed to ask the same question I am asking but ended accepting "Display a standard QProgressBar in a QMessageBar" as an answer. And this isn't what I want.

I am not considering options on how to tell my user to wait. I would like to know if there is a way to call this specific bar.


texts: Is there scripts of QGIS core plugins available anywhere?

texts: Is there anywhere I can see the scripts of the Qgis core plugins.


texts: Merging two polygon layers with one layer having priority over other using ArcGIS Desktop?

texts:
How do I merge two polygon layers?

One is coarse scale over a large area (land use). The other is fine scale but is patchy (agricultural fields) in extent. I want the same extent of the land use layer and have the information of the land use layer but to replace with the fine detail of the agricultural fields where they occur. I am using ArcGIS 10.2.


texts: Why does node Circle to Polygon module returns oval/ellipse instead of a shape like circle?

texts:
I'm using npm package circle-to-polygon and I wrote following code to return a polygon which is approximately like a circle.

const circleToPolygon = require('circle-to-polygon');

let coordinates = [28.612484207825005, 77.23574939044147]; //lat, lng

let radius = 311;
let numberOfEdges = 32;
let polygon = circleToPolygon(coordinates, radius, numberOfEdges);

let coords = polygon.coordinates;

When I used coords and created a polygon shape using Google Maps JS API, it showed a shape like oval (2 ends were very far from each other, two very near, hope you understand).

In fact, the different in distances one can say is roughly 4:1 ratio. So it's never a shape closer to circle.

Is there any issue, or is it supposed to do this intentionally?

See image. The circle you see is actual circle whose center and radius I used. The other shape is the oval I'm talking about. It's way too different.


texts: What are my options for Server side clustering using open source GIS tools?

texts:
I need to make a webapp that will show over 40K Points on a map. These points are all across the globe.  Using WFS Service and OpenLayers, it is possible to do clustering on small set of points, but it is not really useful to do this clustering on the client side for my case.

How can I do server side clustering, in a WMS or tiled WMS service? There are several such options for Google maps. What are the corresponding tools in the OSGeo field?


texts:
I have a postgreSQL database that contains latitude, longitude, and elevation for a region of interest. The same data is available in a geoTIFF, if that would make things easier.

Given two Latitude/Longitude pairs, I need to find the highest elevation along the straight line path connecting the two points. Using Python as a programing language, how would I do this?

Please keep in mind that while I am reasonably fluent in Python and PostgreSQL I am completely unfamiliar with any specialized tools or terms for solving this issue. I've heard of PostGIS, for example, but know nothing about it or how to use it, if it would help here.

Edit:
At the moment, I've managed to figure out a query to find the elevation closest to a given point (since it is unlikely that any given lat/lng pair will match exactly with the data in the database), and of course if there was some query I could do to get all the data points that lie along the path, I could easily get the maximum elevation of those points.

It's finding the point that lie along the path that is baffling me.

EDIT 2: After some additional searching, it would appear that the PostGIS function LINESTRING may be of use, when combined with something else to select all points "along" that line. Still looking for what to combine with it, however.


texts: GeoWebCache, GeoServer and Apache Architecture

texts:
I'm configuring a new distributed GeoServer environment and I have some doubts about GWC Standalone, GeoServer & Apache architecture.

I have the following architecture set in a single developing machine.

Apache HTTPD **HTTP 80**
2 Geoserver Instances (2 TOMCAT) **HTTP 8080** **HTTP 8081**
1 GeoWebCache Standalone Instance (1 TOMCAT) **HTTP 8082**

I have a load balance configuration in Apache HTTPD with mod_jk that points to the 2 GeoServer instances. The 2 GeoServer Tomcat instances are clustered and point to a single DATA DIR.
I've also configured GeoWebCache DATA CACHE DIR.

If I make a request to http://localhost/geoserver it calls one of the 2 GeoServer instances.
If I make a request to http://localhost/geowebcache it calls the GeoWebcache instance (redirects to /geowebcache/home).
I've set the WMS getCapabilites URL from my Geoserver; at geowebcache-core-context.xml

Everything is working fine!

My first doubt is in concern to GeoWebCache and GeoServer WMS endpoints when consuming it in the front-end app (OpenLayers). What is the best approach when calling WMS service?

Make a URL rewrite with Apache so that I call /geoserver/gwc/service/wms/... and it would point to GeoWebCache instance
Call GeoWebCache URL directly: /geowebcache/service/wms/...
Make a URL rewrite to Geoserver WMS URL so that when I call http://localhost/geoserver/<workspace>/wms/... it would point to GeoWebCache instance

Second doubt is about the integrated GWC configuration.

Should I disable GWC at GeoServer Caching Defaults configuration when using GWC standalone?
What really happens when I enable WMS-C? Does it just enables /geoserver/gwc/sevice/wms URL?
And what really happens when I enable WMS direct integration? Does it make /geoserver/wms URL call the integrated GWC?
Do the configuration made in Tile Caching GeoServer panel of each layer affects GWC Standalone (is it passed at getCapabilities)?

The resources I've read on GWC Standalone aren't so clear in the docs. I hope someone have these answers.


texts: spatio temporal anisotropy in R: stAni argument in vgmST and krigeST

texts:
I have some problems in understanding the usage of argument 'stAni' of vgmST and krigeST functions, for the setting of spatio temporal anisotropy.

How have I to interpret this argument, that is the integer value that I saw I have to assign to it, what does it mean?

How can this argument influence the result of prediction?

Are units of spatial and temporal ranges and lags related to this argument?

How can I study the spatio temporal anisotropy of my data?


texts:
QGIS crashes relatively often on my computer with Windows 10.
How can I research what is happening before the crash happens?
Better, it is not crashing with the dialog "ooops, qgis crashes" it lags with "no response in title" ... when then closing the app crashes resp. win says "closing anyway or wait" ...


texts: Error 00358 when using arcpy.SelectLayerByAttributeManagement

texts:
All I'm trying to do is use SelectLayerByAttribute to select all lines less than 5.1 meters. However, no matter what query I try I keep getting
"error 000358 invalid expression" when I run the script as a tool. I have already tried using various combinations of escape characters and string concatenations like ' \"Shape_Length\" < \'5.1\' ' and '\"Shape_Length\" <' + '5.1'. I've also have tried putting Shape_Length in {}'s.  I tried doing the same query in the python window in ArcMap and all I have to do is "Shape_Length < 5.1" and it works just fine. I've never been so stumped. Anyone have anymore ideas on what to try?

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

arcpy.env.overwriteOutput = True
arcpy.SetLogHistory(False)
arcpy.CheckOutExtension("spatial")

og_polyline = arcpy.GetParameterAsText(0)
polyline =arcpy.CopyFeatures_management(og_polyline,r"in_memory/lines")

layer=arcpy.MakeFeatureLayer_management(polyline,"thelayer.lyr")

calc_fields=["SHAPE@","OBJECTID"]

norepeatlist=[]

query="Shape_Length < 5.1"

arcpy.SelectLayerByAttribute_management(layer, "NEW_SELECTION", query)


texts: QGIS Desktop 2.9.0 - ZoomOut doesn't work

texts:
In the QGIS Desktop 2.9.0-Master the Tool "zoomOut" doesn't work. The "zoomOut" have the same operation of the "zoomIn".
This is a bug?
I don't know where I report this bug, because I know the 2.9.0 is the qgis-dev version.


texts: GeoTools Omitting Quotes in its Queries?

texts:
I have a table in my PostgreSQL 8.4 with a column 'Shape' of type 'geometry'.  When I start up Tomcat my GeoServer throws an exception:

24 Sep 16:20:25 WARN [geotools.jdbc] - Error occured determing srid for my_table.Shape
org.postgresql.util.PSQLException: ERROR: column "shape" does not exist

In my PostgreSQL log I can see the corresponding query and error:

ERROR:  column "shape" does not exist at character 77
STATEMENT:  SELECT ST_SRID("Shape") FROM "core"."my_table" WHERE Shape IS NOT NULL LIMIT 1

Just to explore, I tried that same query from psql, and it failed in the same way. I tried it again, except I put double quotes around the 'Shape' in the WHERE clause:

gtda=# SELECT ST_SRID("Shape") FROM "core"."my_table" WHERE "Shape" IS NOT NULL LIMIT 1;
 st_srid
---------
    4326
(1 row)

What is generating that query, why is it omitting the quotes, and is there something I can do about it? I have a number of similar tables, all of which are popping this exception.


texts: How do I compute a python function which calculates/plot minimum distance between a point and a curve?

texts:
I have this code :

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fmin_cobyla
P = (2487,1858)

def f(x):
    return (1.77425666496e-05 * x**3 + -0.125555135823 * x**2 + 296.656015233 * x + -232082.382985)

def objective(X):
    x,y = X
    return np.sqrt((x - P[0])**2 + (y - P[1])**2)

def c1(X):
    x,y = X
    return f(x) - y

X = fmin_cobyla(objective, x0=[2363, 1948], cons=[c1])
print ("The minimum distance is {0:1.2f}".format(objective(X)))

v1 = np.array(P) - np.array(X)

# position vector
v2 = np.array([1, 2.0 * X[0]])
print ("dot(v1, v2) = ",np.dot(v1, v2))
x = np.linspace(1000, 3000, 500000)

plt.plot(x, f(x), 'r-', label='f(x)')
plt.plot(P[0], P[1], 'bo', label='point')
plt.plot([P[0], X[0]], [P[1], X[1]], 'b-', label='shortest distance')
plt.plot([X[0], X[0] + 1], [X[1], X[1] + 2.0 * X[0]], 'g-', label='tangent')
plt.axis('equal')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

The figure that shows computing distance which is obviously wrong.
Can someone tell me what's wrong with my code ?


texts: Mapping geological contacts often needs true contact and inferred contact. I would like to know if it is possible to divide a single polygon outline which marks the outline of a granite body, so that it can be inferred contacts (dashed lines) and true contacts (solid lines)? Basically can different borders of a single polygon be edited separately? I am using ArcMap 9.3


texts: Adding a Dynamic Title in Print Composer

texts: I am using QGIS 3.4.2 for a project with several layers, and a print composer with a generic legend for all. I would like to add a dynamic map title, that states the name of the displayed layer(s) only. Is this possible, and how?


texts: polyline path to follow a draggable modal

texts:
Here is what I am trying to do: on click (right now just click
anywhere on the map)  will draw the line and display a drag-able
modal pop-up, I have static coordinates set for the start point
(which will later be dynamic), but I would like the second point to
find the draggable pop-up modal and stay with it even when it is
dragged.
I started the logic of following the mouse pointer continuously to get the end coordinates. But now thinking not the best method for
finding and following modal.

Below is what I have: arcgis-js

require([
      "esri/Map",
      "esri/views/MapView",
      "esri/Graphic",
      "esri/layers/GraphicsLayer",
      "esri/geometry/support/webMercatorUtils",
      "dojo/dom",
    ],

  function init (Map, MapView, Graphic, GraphicsLayer, webMercatorUtils, dom) {
    var map = new Map({
      basemap: "topo-vector"
    });
    var view = new MapView({
      container: "viewDiv",
      map: map,
      center: [-80.96135253906438, 35.9411934679851],
      zoom: 3
    });

    var graphicsLayer = new GraphicsLayer();
    map.add(graphicsLayer);

    // Create a line geometry
    var simpleLineSymbol = {
       type: "simple-line",
       color: [13,121,190, .9],
       style: "short-dash",
       width: 3
     };

     function showCoordinates(evt) {
        var point = view.toMap({x: evt.x, y: evt.y});
        //the map is in web mercator but display coordinates in geographic (lat, long)
        var mp = webMercatorUtils.webMercatorToGeographic(point);
        //display mouse coordinates
        dom.byId("info").innerHTML = mp.x.toFixed(3) + ", " + mp.y.toFixed(3);

       const coordinatesAx = -80.96135253906438;
       const coordinatesAy = 35.9411934679851;

       const coordinatesBx = mp.x.toFixed(3);
       const coordinatesBy = mp.y.toFixed(3);

       var polyline = {
         type: "polyline",
         paths: [
           [coordinatesAx, coordinatesAy],
           [coordinatesBx, coordinatesBy]
         ]
       };
       var polylineGraphic = new Graphic({
         geometry: polyline,
         symbol: simpleLineSymbol
       })
       //graphicsLayer.add(polylineGraphic);
        document.getElementById("viewDiv").onclick = function fun() {
           graphicsLayer.add(polylineGraphic);
        }
      }
      view.when(function(){
        //after map loads, connect to listen to mouse move & drag events
        view.on("pointer-move", showCoordinates);
        //graphicsLayer.add(polylineGraphic);
      });
});

And for the draggable modal using jQuery Kendo, that I'll call in a separate file (but shouldn't matter how this part is derived just want point B to follow modal/div elem).

  <script>
     $("#dialog").kendoWindow({
        draggable: true
     });
  </script>

here is a visual of what I am trying to accomplish. This would
que on click, and the polyline would follow the draggable modal...*

Mostly just struggling with getting one point of the line to find, and follow my draggable modal. any ideas how I could accomplish this? on pointer-move is clearly not the answer at all.


texts: How to write shapefiles more efficiently in GeoScript?

texts:
I've created a geoscript (in groovy) to create a shapefile from generated geometries with certain attributes.

It's a heatmap consisting roughly 65.000 polygons with an attribute storing the temperature. The final code to write the shapefile looks like that:

dir = new Directory("/my/path")
shp = dir.create('heatmap',
                 [
                     ["geom", "Polygon", "epsg:4326"],
                     ["temp", "float"]
                 ]
                )
// ...
temperatures.eachWithIndex{ temperature, index ->
    shp.add([geometries[index], temperature])
}

The last line is the issue here. Adding each of the 65k attributes to the shapefile is taking ages.

I've found out that (because shapefiles are a proprietary, binary format) adding each geometry causes the shapefile to completely be rewritten from scratch. The result is, that the efford to write records to a shapefile is not linear but exponential to the number of records.

That's a huge problem. Is there any way to write the stuff more efficiently to the shapefile? I got all geometries and all attributes stored in arrays, isn't there any way to write them all at once?


texts: Python script outputs polygon feature class with no shape_area or shape_length and therefore doesn't show up in arcmap

texts:
I am trying to copy features a number of times based on a numeric field value.

The code below works for points when the POINT type is chosen. However, when I switched the feature type to POLYGON it correctly copies the features as many times as I want, but the polygons do not display when brought into arcmap.  For example, the output attribute table (shown below) is populated, but the SHAPE_length and SHAPE_Area fields are shown as zero.  Is there perhaps something in this code maybe within the input or search cursor '*' function that is not reading the shape_area and shape_length fields? The input feature class is a polygon and has a defined shape_area and shape_length for each record, which is not being carried over into the output feature class.

This table is from the output polygon feature class:

def main():
import arcpy
import os
fc_in = r"input.gdb\inputfeatureclass" # this one exists
fld_count = "PHASE_COUNT"
fc_out = r"input.gdb\outputfeatureclass" # this one will be created
arcpy.env.overwriteOutput = True
sr = arcpy.Describe(fc_in).spatialReference

# create the empty output featureclass
path, name = os.path.split(fc_out)
arcpy.CreateFeatureclass_management(path, name, "POLYGON", fc_in, "SAME_AS_TEMPLATE", "SAME_AS_TEMPLATE", sr)

# insert the features into the output fc
with arcpy.da.SearchCursor(fc_in, '*') as curs_in:
    flds_in = curs_in.fields
    idx_cnt = flds_in.index(fld_count)
    with arcpy.da.InsertCursor(fc_out, '*') as curs_out:
        for row in curs_in:
            cnt = row[idx_cnt]
            for i in range(0, cnt):
                curs_out.insertRow(row)

 if __name__ == '__main__':
 main()


texts: Summarize Elevation in ArcGIS Pro not working?

texts:
The tool Summarize Elevation is not producing the typical fields ie Min, Mean, and Max elevation field values, although fields like Dem Res, Product Name, Source, Source URL are being produced.

Has anyone encountered this issue before?


texts: Does CPU with more cores help for lesser rendering time in QGIS 2.18?

texts:
In the settings, there is an option to check to use multicores.
Would this mean the more cores we have, the faster rendering we get?


texts:
Is it possible for two different layers, both based on the same shapefile, to include different features?

I seem to have two layers (using QGIS); one shows some points, and the other shows different points; yet they seem to have the same layer source.


texts: Labeling streets in extremely dense areas using ArcGIS Desktop?

texts:
I'm looking for some suggestions for labeling a a town that includes residential and urban areas that are extremely dense. So far, the only progress I've made has been by using the "reduce font size" feature in the labeling options, but it has still proven to be unsatisfactory and the font reduction engine seems a bit non-intuitive and poorly prioritized. It doesn't make sense to make the size very tiny everywhere, because there are portions of the map that are not dense at all, and lend themselves to a larger font size. I've attached what I have below, and it's a complete mess. Are there any better ways to label a line feature class that is at times incredibly dense?

Below is the old map. Their labeling seems very intelligent, and almost annotation-like (labels look hand-placed with intuitive leaders), and I'd ultimately like to get it looking close to, or exactly like that if I can.

I am using ArcGIS 10.3


texts: QGIS Web Client - restrict search to single schema of Postgres DB

texts:
I have a single PostgreSQL database with multiple schemata, each schema corresponding to a single town and its geodata. In every schema/town, there are the same type of geodata stored, and for every town, I have very similar QGIS project files that display each town in a similar fashion. For every schema, I have a distinct user/role in my PostgresSQL Database whose search path is limited to its corresponding schema.

I want all of those towns to be featured on the same QGIS Web Client landing page (index.php), so that a user can click on any of the towns name and find the corresponding geodata in the QWC GUI.

BUT: Using the search.wsgi and qwc_connect.py modules to allow users to search for addresses in every town, I want to restrict each towns own QWC-GUI to be able to only search within their corresponding schema. That was the reason I implemented the different schemata for different towns, and created distinct roles that can only access their corresponding schema. Otherwise, I feared, a user might query for a very common street name and get the result from an other town that is not the one he clicked on the landing page.
In the '/wsgi/qwc_connect.py' -file, I have to specify a database connection. And as it presents itself, I can only specify a single database role to login to the DB:

#!/usr/bin/python
# -*- coding: utf-8 -*-
# Script to connect to a PostgreSQL database. Connection object to be used in all wasgi scripts

import psycopg2 #PostgreSQL DB Connection

# configure your DB connection here
DB_CONN_STRING="host='localhost' dbname='mydb' port='5432' user='myrole' password='************'"

def getConnection(environ, start_response):
  #SQL database connection
  try:
    conn = psycopg2.connect(DB_CONN_STRING)
    return conn
  except:
    errorText = 'error: database connection failed!'
    # write the error message to the error.log
    print >> environ['wsgi.errors'], "%s" % errorText
    response_headers = [('Content-type', 'text/plain'),
                        ('Content-Length', str(len(errorText)))]
    start_response('500 INTERNAL SERVER ERROR', response_headers)

    return None

How I understood it, the variable 'searchtables' in the file '/wsgi/search.wsgi' will access a indexed table that is stored within the schema that was specified by the QGIS project file. So there is no need for editing code there, if I give the searchtables in every towns schema the same name, and restrict the access to each schema to the corresponding role (which is also the role used to establish the DB-connection the QGIS project file). But how do I specify multiple roles in qwc_connect.py, so that only the role that is currently used by the QGIS project file to access its own schema is also used for querying the database?

I was thinking about distinct qwc_connect.py files for every QGIS project file, but don't know how to link them. Nor do I know how to rewrite the code above to include multiple db-connections that become active depending on the QGIS project file from where the query originated. Can you show me how it could be done?

What I don't want is to create separate instances of QWC, or separate databases for every town. That would be an easy solution right now, but uncomfortable to manage in the future and difficult to maintain when more and more towns are added to the mix.


texts: Tabulate intersection, but I need it with null values too, not just with matching intersect

texts: I have problem with tabulate intersection tool (In ArcGIS - ArcMap 10.7.1). The number of data is 1298 with 500 m buffer. Then i have polygons, that I need to process. I run the tool,  but get the table with only 170 records, that is matching, but it is not all what i need. I need the complete data even with null fields, where the tool doesn't find the intersect. Is here any way, how to get table with the null values and intersect values together in tabulate intersection?


texts: Couldn't add a column to attribute table (GPX layer) | QGIS

texts:
Running QGIS 2.10.1 on Ubuntu 14.04 x64.

Going to add a column to the attribute table (GPX layer) but couldn't select any type. Then, failed - "Failed to add field 'NameTest' of type. Is the field name unique?".

Need to provide types firstly before adding a column? How to set it?

Reference
https://docs.qgis.org/2.2/en/docs/training_manual/create_vector_data/actions.html


texts: Union on polygons with holes using JSTS?

texts:
I'm trying to union two polygons (see above) with JSTS. My first polygon contains an hole, and my second polygon fill partially this hole.

The expected result is a unique polygon with two hole (or one hole with two points at the same location.

I've tried with PostGIS and the result is two holes). The actual result is what you can see on C: my first polygon has no hole any more.

Is it a defect of the JSTS tool or is it just some configuration to add somewhere?

I'm using the JSTS union function like that:
firstGeometry.union(secondGeometry);

I don't think I have topological errors. I'm easily able to reproduce this case by drawing polygons on OpenLayers, and my polygons seems to be correct. I'm using JSTS and not JTS. I've tried with the buffer(0,0) function, but I have unfortunately the same result.

Here is an example of a union that fails:

POLYGON((335179.0532258 328356.69669355,335178.1782258 328437.63419355,335285.3657258 328437.19669355,335284.9282258 328355.82169355,335179.0532258 328356.69669355),(335213.6157258 328373.32169355,335198.3032258 328400.88419355,335225.4282258 328426.69669355,335267.4282258 328415.75919355,335262.1782258 328371.13419355,335213.6157258 328373.32169355))

POLYGON((335213.6157258 328373.32169355,335225.4282258 328426.69669355,335267.4282258 328415.75919355,335213.6157258 328373.32169355))

And the result of a JSTS union() between these two geometries or a buffer(0,0) on a new geometry collection containing these geometries is the following polygon :

POLYGON((335179.0532258 328356.69669355,335178.1782258 328437.63419355,335285.3657258 328437.19669355,335284.9282258 328355.82169355,335179.0532258 328356.69669355))

It happens in each case where I have

1 to many hole(s)
1 to many polygon(s) within this/these hole(s), creating "sub-holes"
Sub-holes are sharing a vertex

Doesn't seems to happen in other cases.
Geometries are topologically correct.

I've found a defect about that on GitHub: https://github.com/bjornharrtell/jsts/issues/135
I'm still looking for a way to fix it!


texts: Dynamically change Mapnik XML without restarting renderd

texts:
I have multiple tables of data in my PostgreSQL database that have wkb_geometry columns used to generate layers. I would like to allow users to select which layers on a map they will see by accessing UI components such as checkboxes. When doing so, based on the selected UI components, it should send a request to my tile server and update the Mapnik XML file accordingly before rendering the tiles. How can I do this without having to restart renderd each time and how can I make sure that concurrent accesses won't overwrite changes to my XML file?

I do not want to cache these tiles as the same URL could have different parameters (layers turned on) so how would I avoid caching these requests?


texts: Getting DEM Values as NumPY Array in Earth Engine

texts:
I'm trying to convert Earth Engine data to a local GeoTIFF manually, using more or less the method found here. One crucial step is converting the DEM raster values to a NumPy array. I've been trying to do so as follows:

ee.Initialize()

# define ROI
area = ee.Geometry.Rectangle(bounds)

# set image
img = ee.Image('USGS/SRTMGL1_003')

# get the lat lon and add the elevation data back in
latlng = ee.Image.pixelLonLat().addBands(img)

# reduce to list
latlng = latlng.reduceRegion(reducer=ee.Reducer.toList(), geometry=area, maxPixels=1e8, scale=20)

# put the values into lists
elev_values = np.array((ee.Array(latlng.get("elevation")).getInfo()))
lats = np.array((ee.Array(latlng.get("latitude")).getInfo()))
lngs = np.array((ee.Array(latlng.get("longitude")).getInfo()))

While I have had success converting unique latitude and longitude values to a NumPy array, I keep ending up with an empty array of elevation values. The current code throws the following error:

Traceback (most recent call last):
  File "test_one.py", line 392, in <module>
    gDem = GoogleDEM([40.01, -21.01, 40.02, -21.00])
  File "test_one.py", line 41, in __init__
    elev_values = np.array((ee.Array(latlng.get("elevation")).getInfo()))
  File "C:\Users\lmonn\AppData\Local\Programs\Python\Python37-32\lib\site-packages\ee\computedobject.py", line 95, in getInfo
    return data.computeValue(self)
  File "C:\Users\lmonn\AppData\Local\Programs\Python\Python37-32\lib\site-packages\ee\data.py", line 490, in computeValue
    return send_('/value', ({'json': obj.serialize(), 'json_format': 'v2'}))
  File "C:\Users\lmonn\AppData\Local\Programs\Python\Python37-32\lib\site-packages\ee\data.py", line 1186, in send_
    raise ee_exception.EEException(json_content['error']['message'])
ee.ee_exception.EEException: Array: No numbers in 'values', must provide a type.

My guess would be that there is some nesting of objects taking place that I don't quite have a grasp on.

How can I fix this issue? Is there a simpler way to get a NumPy Array of the elevation values for a DEM on Earth Engine?


texts: Conversion of EPSG:3857 coordinate to UTM

texts:
I'm trying to convert coordinates from Google Maps using PostGIS.

I'm using this query:

SELECT ST_AsText(ST_Transform(ST_GeomFromText('POINT(-29.16297 -51.1795034)', 900913), 29182));

The problem is that when I put a point on the map using the UTM coordinates, it is in a different place.

Is there a chance of one of these SRIDs be wrong?


texts: PostGis stopped working: Invalid access to memory location

texts:
I have a PostgreSQL v9.3 database with PostGIS 2.1.8 running locally on Win10 Creative Update 64-Bit.
This database is running for months without any problem. I've used it yesterday afternoon.
This morning when I want to query the table with a geometry column I get this error:
ERROR:  could not load library "D:/Program Files/PostgreSQL/9.3/lib/postgis-2.1.dll": Invalid access to memory location.

When I use pgAdmin to connect to a remote test-server with the same database I don't have this error. So it must be something with my local instance.

Of course, I already restarted my laptop and restarted the PostgreSQL service, without any luck. I checked my Windows Updates and nothing has been installed in the last two weeks.

The dll-file does exist, I have verified it.
I have no idea what is happening and how to solve this.


texts: QGIS 3.12 - Setting Project default CRS causes "Ballpark Transform"

texts:
I'm coming across an issue that I am unsure of whether it is a bug or the result of my error. My work data is all projected in NAD83 (CSRS) / MTM Zone 10. EPSG Code: 2952.

I set QGIS to default to this projection. This was done on a blank/empty project ->Settings->Options->CRS->Set Default CRS. When loading QGIS (Empty project or not), I always get the following error:

At work where I use ArcMap, I haven't experienced this issue at all so am confused on how to address it nor what exactly is the cause.


texts: Get pixel geo info according projection and map extent

texts:
I would like to know if it is possible to calculate the geographic coordinates (lat, lon) for a pixel specified by its x/y (top-left) cartesian coordinates in a map image/raster?

The datum and projection are known: WSG84/EPSG4326 , as the envelope/extent (4 "corners" + the center point).

It seems the goal is somewhat similar to the one of gdallocationinfo, but I don't find the way to indicate the map extent/bbox.

If possible, tools from GDAL or GeographicLib can be used.

I can't find the related formula.


texts: Generate a geometry column from latitude and longitude columns with SQL

texts:
My table has the columns latitude and longitude in float8 type.
I want to generate a geometry column from these two columns.
I make the following query in the SQL Window of the DB Manager in QGIS.
ALTER TABLE populated_places
ADD geographyColumn AS geography::STGeomFromText('POINT('+convert(varchar(20),latitude)+' 'convert(varchar(20),longitude)+')',4326)

I get the following error: syntax error at or near "AS"
How can I solve this problem?


texts: Selecting layer in Road Graph plugin settings in QGIS?

texts: I am currently using QGIS 2.10.1. and I am not able to select any layer when trying to configure the Road Graph plugin. What should I do?


texts:
I'm not a programmer and a fairly outdated mac-user. However, got a new mac again and tried to install QGIS on High Sierra. It opens, but there is an error message about Metasearch:

  2018-01-16T16:59:18   1   Traceback (most recent call last):
              File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/qgis/utils.py",
  line 333, in startPlugin
                plugins[packageName] = package.classFactory(iface)
              File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/MetaSearch/init.py",
  line 29, in classFactory
                from MetaSearch.plugin import MetaSearchPlugin
              File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/qgis/utils.py",
  line 607, in _import
                mod = _builtin_import(name, globals, locals, fromlist, level)
              File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/MetaSearch/plugin.py",
  line 33, in
                from MetaSearch.dialogs.maindialog import MetaSearchDialog
              File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/qgis/utils.py",
  line 607, in _import
                mod = _builtin_import(name, globals, locals, fromlist, level)
              File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/MetaSearch/dialogs/maindialog.py",
  line 43, in
                from owslib.csw import CatalogueServiceWeb
              File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/qgis/utils.py",
  line 607, in _import
                mod = _builtin_import(name, globals, locals, fromlist, level)
              File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/owslib/csw.py",
  line 27, in
                from owslib.util import OrderedDict
              File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/qgis/utils.py",
  line 607, in _import
                mod = _builtin_import(name, globals, locals, fromlist, level)
              File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/owslib/util.py",
  line 35, in
                import requests
              File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/qgis/utils.py",
  line 607, in _import
                mod = _builtin_import(name, globals, locals, fromlist, level)           ImportError: No module named requests

How (and where) can I fix that problem?


texts: Getting value from another table with dbquery

texts:
I am trying to use dbquery to get a value from a table when a record in it is the same as a value in the source table.

In the table where I need the value-added (Sewerage_Manholes...) I am trying

dbquery('asset','commision_date','GID_New=GID')

What am I not doing correctly?

I don't want to run a join by attribute in this case as if this works it opens up the opportunity to use complex SQL.


texts: Query values from list using Select_analysis() in ArcPy?

texts:
I have a list like this...

list = ['1','2','3','4']

how can I use arcpy.Select_analysis to select features that have the attributes in the list.

qry = '"SECDIVNO" IN ' + str(tuple(list))
arcpy.Select_analysis(tempFeat_1, tempFeat_2, qry)

Is the "tuple" incorrect?

I've tried this also and it doesn't work either...

qry = '"SECDIVNO" IN ' + list
arcpy.Select_analysis(tempFeat_1, tempFeat_2, qry)

The error I get is always the "invalid expression error"


texts: Rotate and edit Legend in layout (Python-QGis)

texts:
I'm trying to set up a legend in the lower part of the map and I have only achieved this:
legend = QgsLayoutItemLegend(layout)
legend.model().setRootGroup(root)
layout.addLayoutItem(legend)
legend.attemptMove(QgsLayoutPoint(210, 5, QgsUnitTypes.LayoutMillimeters))

This creates the legend in the layout but I want to rotate it 270°, delete the margins around the colored rectangles, hide the subgroup's title and hide the white color of the background (so that the legend's background appears transparent).
I would also like to make the reduce the Spacing between legend items to 0.05 mm and change the size of the legend rectangles.
The pictures show the code result and the manually edited image.
Thanks in advance!


texts: How can I set a layer visible dynamically using Flex Viewer?

texts:
I have a set of layers in operationallayers tag in my Flex Code connected to ArcGIS Server 10, and I want to make them visible dynamically in a ComboBox Change Handler.

Is it possible?


texts:
I am writing a Python script to reproject a KML file to WGS 84 and then to convert it to a GeoJSON. I am aware of and use the GDAL command line tool ogr2ogr (which was the answer to this question in this post from 2012).

I can use subprocess to use the ogr2ogr from my script, but I was wondering if there are any libraries out there that might handle this better/make it easier for me. Is there anything else out there? Or is ogr2ogr my best option?


texts: How to define a CRS using a projection rotated about a given point

texts:
I have a data set that uses a common projection, but has its own false easting and northing, and is rotated 315 degrees. I need to be able to work in this CRS and easily translate from one CRS to another, so using hacks such as CAD tools is not practical. I have a definition of it for MapInfo, which works fine. However, I need it to also work in QGIS and other software, which require .proj4 format. I have defined a CRS which gives the results I want, but when changing from one projection to another, there is an offset generated which is different depending on the projection being used.

The MapInfo definition is:
"GRDMF", 3018, 31, 7, 173, -41, 2510000, 6023150, 7, 0.686785, 0.726871, -5539726.41, -0.726876, 0.686796, -2108424.34, -1000000, -1000000, 1000000, 1000000

The .proj4 definition I have is:
+proj=nzmg +lat_0=-41 +lon_0=173 +x_0=2510000 +y_0=6023150 +no_defs +a=6378388 +rf=297 +ellps=GRS80 +datum=nzgd49 +nadgrids=nzgd2kgrid0005.gsb +towgs84=84,-22,209,-0,-0,-0,0
+to_meter=1

The projection is New Zealand Map Grid with NZGD1949. The CRS I need is based on NZMG projection, but rotated 315 about the coordinate 2310655.80 5535100.29, which translates to 70501.32 13497.54 in the custom CRS.

Some test coordinates are:

NZMG coords            Custom CRS
2292000 5539000        60522.871 29736.79
2296000 5552000        72720.155 35757.861
2330000 5527000        77898.623 -6126.353
2339000 5521000        79718.234 -16788.504
2333000 5511000        68329.450 -19294.963
2321000 5520000        66629.725 -4392.060

Any help would be appreciated.
Thanks


texts:
I have percent sand/silt/clay values for 30 sample points on grid. I used these values to create 3 interpolated IDW rasters. I want to somehow (I think in raster calculator) to make a code using if/then statements to take values from each IDW and output a texture.

Example: if raster "sand_IDW" value > 50 and < 70 and raster "clay_IDW" value > 20 and <40 then, "Texture" = Sandy clay.


texts: Q-Chainage plugin. Chain points layer does not save properly

texts:
Using Windows 10, QGIS latest version.
Use Q chainage (vector plugin) to make a chain of points on a boundary (line vector layer).
I save the points layer (as layer definition file), and also save the whole project.
When I close and then return to the Project, the points are not there.
The name of the layer is listed in the layers, but it is empty.


texts: Quarterly composites for calculating NDVI and EVI

texts:
I am trying to create a code that could be helpful for making quarterly composites from 2013 to 2016 in order to calculate NDVI and EVI.
I have figured out this code:
var start = ee.Date('2013-01-01');
var end = ee.Date('2019-01-01');
var numbQuarters = end.difference(start, 'month').divide(3).ceil();
//NDVI
var ndvi = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
  .filterBounds(table)
  .filterMetadata('CLOUD_COVER','less_than',30)
  .map(function(image) {
    return image.select().addBands(image.normalizedDifference(['B5', 'B4']));
  });
//EVI?
// make a composite image for every quarter
var quarterlyImages = ee.ImageCollection.fromImages(
            ee.List.sequence(1, numbQuarters).map(function(quarter){
              var startTemp = start.advance(ee.Number(quarter).subtract(1).multiply(3), 'month');
              var endTemp = start.advance(ee.Number(quarter).multiply(3), 'month');
              //var image = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
                                //.filterDate(startTemp, endTemp)
                                //.filterBounds(table)
                                //.filterMetadata('CLOUD_COVER','less_than',30)
                                // define a way to composite an image
                                //.mosaic(); // .mean(); .mossaic(); .max(); etc
              var subset = ndvi.filterDate(startTemp,endTemp);

              //return image.set('system:time_start', startTemp.millis(),
                               //'system:time_end', endTemp.millis());
              return subset.mean().set('system:time_start',startTemp);
            }));

//var median = quarterlyImages.median();
//print(quarterlyImages,median);
//Map.addLayer(median,{ min:-1, max: 1, palette: ['blue', 'white', 'green']},'median');
//Lay anh trong imageCollection
var sort = quarterlyImages.toList(100);
var clip1 = ee.Image(sort.get(1)).clip(table);//1st image 2nd quarter of 2013
var clip2 = ee.Image(sort.get(2)).clip(table);//2nd image 3rd quarter of 2013
var clip3 = ee.Image(sort.get(3)).clip(table);//3rd image 4th quarter of 2013

// print some stuff
///print('image collection', quarterlyImages)
//print('startDates', quarterlyImages.aggregate_array('system:time_start'))
//print('sendDates', quarterlyImages.aggregate_array('system:time_end'))
Map.addLayer(clip1,{ min:-1, max: 1, palette: ['blue', 'white', 'green']}, 'Phu Quoc 01');
Map.addLayer(clip2,{ min:-1, max: 1, palette: ['blue', 'white', 'green']}, 'Phu Quoc 02');
Map.addLayer(clip3,{ min:-1, max: 1, palette: ['blue', 'white', 'green']}, 'Phu Quoc 03');
Map.centerObject(table, 7);

however when I apply this code above to extract images, some of them have too many clouds so that I could not see the image clearly.
I am thinking of using a cloud masking function to solve this problem but I assume that I did not make it correctly and therefore it did not work well. Here is the code:
var cloud_pixel_qa =  image.select('pixel_qa').bitwiseAnd(32).neq(0);//2^5
var cloud_shadow_pixel_qa =  image.select('pixel_qa').bitwiseAnd(8).neq(0);//2^3
var maskedClouds = ((cloud_pixel_qa).or(cloud_shadow_pixel_qa));
return image.updateMask((maskedClouds.add(1).unmask(0)).eq(1));
var L8_coll_mask=L8_coll.map(masking_cloud)
Map.addLayer(L8_coll_mask.median().clip(LĐ_T),vizParams,'false color composite cloud mask');

Does anyone have any suggestions?


texts: PostGIS 2.0 ST_Split rounding errors when splitting line with point

texts:
My apologies for the bad title; I couldn't think of anything more descriptive.

I have been testing ST_Split with a dataset, and found that even with points that are supposed to be exactly on a line (found by ST_Line_Interpolate_Point) often the line doesn't get split.

select ST_AsText(ST_Split(ST_GeomFromText('LINESTRING(0 0, 100 100, 200 200, 400 400)'), ST_Line_Interpolate_Point(ST_GeomFromText('LINESTRING(0 0, 100 100, 200 200, 400 400)'), 0.2)))

Result:

GEOMETRYCOLLECTION(LINESTRING(0 0,80 80),LINESTRING(80 80,100 100,200 200,400 400))

Which is as  expected. Now I introduce some decimals:

select ST_AsText(ST_Split(ST_GeomFromText('LINESTRING(0 0, 100.123456789 100.987654321, 200 200, 400 400)'), ST_Line_Interpolate_Point(ST_GeomFromText('LINESTRING(0 0, 100.123456789 100.987654321, 200 200, 400 400)'), 0.2)))

Result:

GEOMETRYCOLLECTION(LINESTRING(0 0,100.123456789 100.987654321,200 200,400 400))

Stripping some decimals:

select ST_AsText(ST_Split(ST_GeomFromText('LINESTRING(0 0, 100.123456 100.987654, 200 200, 400 400)'), ST_Line_Interpolate_Point(ST_GeomFromText('LINESTRING(0 0, 100.123456 100.987654, 200 200, 400 400)'), 0.2)))

Result:

GEOMETRYCOLLECTION(LINESTRING(0 0,79.6558670241487 80.3434026298966),LINESTRING(79.6558670241487 80.3434026298966,100.123456 100.987654,200 200,400 400))

Nearly all records in my dataset have a lot of decimals (which I think the problem is), meaning NONE of my lines get split, even though the point is supposed to be exactly on the line. Is there any way to work around this (quickly; I have to split 300k lines like this)? Preferably without modifying my dataset by using ST_SnapToGrid?

I'm Using POSTGIS="2.0.3 r11132" GEOS="3.3.8-CAPI-1.7.8" PROJ="Rel. 4.8.0, 6 March 2012"


texts:
I'm trying to turn a QGIS 2.4 project into a TMS tileset.  What method have folks used, with what success and what level of effort?

The project in question has about 50 layers and uses OTF projection since different layers are in different CRSs.

the QTiles plugin hangs (and takes down the rest of QGIS with it) at 'Rendering: 0 from ... (0%)' even for a fairly small area and just a couple of zoom levels; I think this is just because the render function hangs in the same manner if you call it directly?  Is there any way to avoid this problem with QTiles?
gdal2tiles, using geotiff rendered from QGIS as input - I'd have to assemble a VRT and use that as input since QGIS errors out if I try to render one geotiff of even a small fraction of the size needed, so, this would be a question of writing a script to render many smaller geotiffs and building a corresponding VRT, right?  A current problem I'm having with this is trying to effectively 'apply OTF projection for the render function' since the render output only has the data of whichever layers' native CRSs result in data in the area defined by the QgsRectangle; is there an easy way to do that?
a slight variant of the above - skip the gdal2tiles step and write a script to render 256x256 .png tiles directly from QGIS
QGISServer and tilecache - I haven't tried this option but the overhead of setting up the local web server seems daunting; has anyone used this option with good success?
any other plugins or (simpler) methods to accomplish this task?


texts: Smoothing or resampling DEM to find steep roads?

texts:
I am trying to find sections of road with steep gradients.

My process is:

pull a linestring from Openstreetmap data
resample it to create a point every 10 metres
look up the height data for each point from ASTER Global Digital
Elevation raster
calculate the gradient at 10m intervals and look for any in excess
of 14%

This works but throws up a LOT of anomalies which must be due to the height data as I have ridden my bike along these sections. One 10m section on a road segment local to me shows a height difference of 14m between two points 10m apart on the line.

Now clearly this is due to the height data. It is stated to have an accuracy of:

  20 meters at 95 % confidence for vertical data and 30 meters at 95 %
  confidence for horizontal data.

How should I go about smoothing out these anomalies and having a best guess at where the steepest road gradients really are?

Everything has been written in PGPLSQL/POSTGIS so far. I have a set of line strings with points every 10m and a height for each point.


texts: Duplicating layer-table in GeoPackage layer in QGIS3

texts: How can I create new Geopackage layer with the same styling and fields, but empty of features, within same dpkg file. I have yet, houses.dpkg and within this dpk I have created houses and fences tables (layers), but now I want to create houses2 and fences2 within the same dpkg, but with same styling and table structure (fields), as the previous two tables=layers, but to be independet (not linked in any way) with their originating layers (houses and fences), except they are stored in the same dpkg file.


texts: Extracting color values from layer to build index using ArcPy/ArcMap?

texts:
My university wants the maps we create to become more standardized and so we're changing all the colors on our base vector map.

The new standard needs to be indexed and I have the task of creating the spreadsheet/document that will help us in the future.

There are close to 100 layers that I need to index with, for example, the polygon layer name, fill color (rgb values), border color, and border width.

I have only recently started learning ArcPy so I was hoping there would be a GUI method I could follow and if not, what I should be looking for in an ArcPy script.

How does one extract this information?


texts:
I have 4 Lat Lon positions. They represent a (rotated) rectangular region. But due to some calculation errors, sometimes they are not proper rectangular.

I want to convert those into (rotated) rectangular regions.
See Image for more clarification

In the above image I want the lat lon of the corners of the rectangle drawn given 4 lat lon as shown on the left side.

I tried searching but couldn't find any solution. Most of the solutions did not take care of rotated rectangle case.
How do I solve this problem ?

Note : Are there any libraries that does this, I am using JavaScript for implementation.
Note : I don't know if my question is more related to math.


texts:
I am labeling a 2-D area with several intersecting dimension labels, and I would like to move the dimension text ('512.64 ft') over to the sides of the polygon to make it easier to interpret.

However, I can't find a way to move these labels, nor can I convert the labels to annotations like I would with a normal feature class.

Using ArcMap 10.3.


texts: Looking for free web service that calculates driving distance between 2 addresses

texts:
I'm looking for a free web service that gets an input of 2 addresses and returns an output
of driving distance between the 2 points. In addition I'm looking for a web service that doesn't require submitting a domain in order to use it (like google's distance matrix web service).

Thanks.


texts: How to load and view .gri and .grd files in QGIS

texts:
An online GIS Repository called DIVA-GIS (http://www.diva-gis.org/gdata) offers data in three formats - ERSI files, "DIVA gridfiles", and DBF files

The DIVA gridfiles have extensions .gri and .grd

How can these file types be opened in qGIS?  The "Add Raster" menu option does not seem to work


texts: Is there an industry-standard official mapping of Galileo satellites to global "PRN”/ID values?

texts:
In the U.S. GPS system, satellites are identified by the PRN number broadcast by each satellite.  This ID is usually visible on most GPS receiver platforms that provide access to individual GPS satellite signal information to identify the broadcasting satellite (e.g., accessible on Android via the GPSStatus.Listener and the GpsSatellite.getPrn() method (http://goo.gl/QLse3)).
For other GNSS, there have been an adoption of global PRN/ID values for satellites so apps end up getting globally unique IDs per satellite that don't collide with other GNSS.
For example, according to an unofficial version of the NMEA 0183 spec (v2.19, Apr 2015):

Russian GLONASS - IDs 65-96 reserved (See also https://gis.stackexchange.com/a/91629/12198)
Japanese QZSS - IDs 193-200 reserved
Chinese BeiDou - IDs 201-235 reserved

As of April 2015 version of the NMEA 0183 spec above, Galileo isn't mentioned as far as ID reservations.
Is there any industry-standard official mapping between Galileo satellites and PRN/ID values?
EDIT Aug. 9, 2016
According to GPS World:

According to two Notice Advisories to Galileo Users (NAGUs), the two Galileo satellites launched into elliptical orbits in August 2014, GSAT0201 using PRN code E18 and GSAT0202 using PRN code E14, will start transmitting navigation messages for test purposes this Friday, Aug. 5.
The Signal Health Status (SHS) flags will be set to “Test” and the Data Validity Status (DVS) flags will be set to WWG (working without guarantee).
The satellites will not be included in the broadcast almanacs.

These PRNs are alphanumeric, while the existing agreed-upon ranges for other GNSS are numeric.  I'm looking for the numeric PRNs for Galileo that would complement the other existing known GNSS ranges and that would be available on Android.
If you have an Android device that supports Galileo and you're willing to help identify the ID range, please comment on this Github issue.


texts:
Is there a way to specify a feature of a shapefile or all features within a shapefile to serve as masking features for all other layers within a qgis project?
I would like to crop all other features of other layers to the specified extent (hide everything that falls outside of the mask) without having to use geoprocessing tools. I know there is a mask plugin, but that does not allow to achieve what I try to achieve. This should work not just in print composer through the atlas.

It would be great to have that option within the layeer options or even the projects settings.


texts: Removing Frame Border from USGS scanned images in ArcGIS with Python?

texts:
I have a couple of USGS scanned maps like

I want to remove the frame border from these images so that I can use them in ArcGIS.

Is there algorithm or python script that I can run to remove the borders or get the extents of the actual data?


texts: ArcGIS clip tool between a polyline and a polygon does not clip all features

texts:
I am using ArcGIS 10.5 Clip tool to locate flooding on the road with a polygon file for the flooding (Clip Features) and a polyline file for the road network (Input Features). The polygon file consists of many triangularly shaped ponds. To filter it a bit, I have selected only the polygons, intersecting the road and have copied them to a new layer file. When I run the Clip tool, most of the flooded roads sections are extracted correctly, but others are not clipped at all as if some of the polygons are disregarded. See the image for more information.

The blue polygons are ponds of water on the road and the green lines are the result of the Clip tool. The road network is not interrupted.

I observed it visually and could not draw any conclusions why it might be happening and how to fix that. Performed the same thing with Intersect and it selected whole roads instead of the flooded segments of the road.


texts: Erasing line feature (erase feature) from polygon feature in ArcGIS Desktop?

texts:
How can I erase a "line feature" (erase feature) from a "polygone feature" (input) in ArcGIS 10?

With the tool "erase" I can only use an other "polygon feature" as erase feature for a polygon (input).But I can not use an "line feature" as erase feature for a polygon.

I there an other tool in ArcGIS or an other way to solve this problem?


texts: Using integer variable in where clause of ArcPy script?

texts:
I'm trying to include an integer variable in the SQL where-clause of a tool like this (ArcGIS 10 Python script):

newR = ExtractByAttributes(inR, '"IntField">=2508')

where inR is an integer raster, IntField is an integer field in its attribute table.

The above works, but how can I substitute an integer variable instead of 2508?

I've tried all kinds of different quotes and concatenation, but have not been able to.


texts: Display attributes of feature with clustering

texts:
I have one KML layer that has roughly 900 features and I added OpenLayer.Cluster.Strategy to help the layer display a little quicker.

My KML displays the attributes when clicked on before clustering was enabled, but now the bubble displays 'undefined, undefined'.

How did the clustering change the extractAttributes property?

            var state = new OpenLayers.Layer.Vector("State Probation", {
            projection: map.displayProjection,
            strategies: [new OpenLayers.Strategy.Fixed(), new OpenLayers.Strategy.Cluster()],
            protocol: new OpenLayers.Protocol.HTTP({
                url: "state.kml",
                format: new OpenLayers.Format.KML({
                    extractStyles: true,
                    extractAttributes: true
                })
            })
        });

        select = new OpenLayers.Control.SelectFeature([sample, state],
            {
                clickout: true, toggle: false,
                multiple: true, hover: false,
                toggleKey: "ctrlKey", // ctrl key removes from selection
                multipleKey: "shiftKey" // shift key adds to selection
            });

        sample.events.on({
            "featureselected": onFeatureSelect,
            "featureunselected": onFeatureUnselect
        });

        state.events.on({
            "featureselected": onFeatureSelect,
            "featureunselected": onFeatureUnselect
        });

        map.addControl(select);
        select.activate();

        map.addControl(new OpenLayers.Control.LayerSwitcher());

    function onPopupClose(evt) {
        select.unselectAll();
    }

    function onFeatureSelect(event) {
        var feature = event.feature;
        var selectedFeature = feature;
        var popup = new OpenLayers.Popup.FramedCloud("chicken",
            feature.geometry.getBounds().getCenterLonLat(),
            new OpenLayers.Size(100,100),
            "<h2>"+feature.attributes.name + "</h2>" + feature.attributes.description,
            null, true, onPopupClose
        );
        feature.popup = popup;
        map.addPopup(popup);
    }
    function onFeatureUnselect(event) {
        var feature = event.feature;
        if(feature.popup) {
            map.removePopup(feature.popup);
            feature.popup.destroy();
            delete feature.popup;
        }
    }


texts: ArcGIS Pro/Online: deleting a field "General function failure[Incorrect syntax near ','. Incorrect syntax near','.]"

texts: In ArcGIS Pro I accidentally deleted two fields by clicking the "Delete Field" button at the top of the attribute table. The two fields were "GlobalID" and "Appointment".  I was trying to delete selected features at the time. Now I can't edit or delete any features in this layer. I get an error message of: "Update table row. Failed to "Update table row". General function failure[Incorrect syntax near ','. Incorrect syntax near','.]".  If I click on one of the features in the Map window the deleted fields still appear, but instead of having information about them they just say the field name again. Preferably I could restore these fields and go back to where I was previously. Otherwise, how do I push this delete through to the individual features. "GlobalID" has no bearing on my workflow, and "Appointment" will only affect a few features at this point.


texts:
Does any one know where can I download the 2.18 nightly version of qgis?
and would you give me some information about it?


texts: What is the wire format for Google's new vector map tiles?

texts: Google recently released new vector cartography for their mobile maps. I'm curious if anyone here has any insight into the wire format they use for transmission? I've experimented with a few formats over the years, including bencode, bson, variants of WKT and WKB, and gzipped, limited-precision geojson which I've found to be the best trade-off in comprehensibility and file size. I'd love to know what Google decided on, but I don't have an Android device to test with.


texts: How to improve poor performance times for PostGIS 2.1 Tiger Geocoder?

texts:
I've recently upgraded my PostGIS DB from 2.0 to 2.1.  I initially followed the instructions from Waiting for PostGIS 2.1 - Install PostGIS Tiger Geocoder as an Extension and  2.7. Installing, Upgrading Tiger Geocoder and loading data to load (for the first time) the 2012 Tiger Geocoder into the DB by grabbing the updated geocoder files from the postgis-pg92-binaries-2.1.1devw32.zip file.  Then I downloaded and installed the data for 5 western states. (OR, WA, ID, CA, MT).  Finally, I ran the SELECT install_missing_indexes(); command after loading all of the data. I also ran vacuum analyze on the whole DB afterwards.  That all went fine, and the geocoding process is working.

However, while I am admittedly running this DB on a not-so-speedy laptop, my results of 2-8 sec results versus the 61ms results noted here seem a way out of whack.

Any ideas on what I can do to improve performance?

Guides I worked through to try to improve performance:

OpenGeo "Section 29: Tuning PostgreSQL for Spatial"
OpenStreetMap "Tune the database "
Part 1: Getting Started With PostGIS: An almost Idiot's Guide (PostGIS 2.0)
Best practices for PostGIS production environments?

and I've updated all of my postgresql.conf settings to match the highest settings from any of the documents.

System specs: Win 7 64 bit - Intel 2.3Ghz, 4GB RAM, Paging space 2000MB to 8000MB (fixed), PostgreSQL 9.2 (32 bit), Postgis 2.1

Additional Details:  When I run the following query using "explain analyze" I get the following results:

SELECT g.rating, ST_X(g.geomout) As lon, ST_Y(g.geomout) As lat,
    (addy).address As stno, (addy).streetname As street,
    (addy).streettypeabbrev As styp, (addy).location As city, (addy).stateabbrev As st,(addy).zip
    FROM geocode('<house #> <street name>, <city> <state> <zip>') As g;

Function Scan on geocode g  (cost=0.25..15.25 rows=1000 width=68) (actual time=185.606..185.606 rows=1 loops=1)


texts: Shapefile and DEM: check rivers behaviour

texts:
I'm working with shapefiles representing rivers (rivers.shp), and DEM for the terrain elevation (terrain.tif).

Basically, I want to check that my rivers (stored as polylines in my shapefile) behave like "natural" water flow, meaning that it mainly go downstream, can not go uphill first and downhill then, etc...

I'm using GDAL/OGR in Python, and don't really see how to proceed.
Should I extract the points of the polylines in the shp? Look for the elevation at those points in the DEM, and then compute my decision?

Thanks,
eo


texts: SDE GDB Versioning - reconcile and post from DEFAULT to Child

texts:
I am relatively new to versioning and think I have made a mistake in my workflow.

I am using a SDE geodatabase, with a DEFAULT and 'Field Editing' child version in ArcGIS 10.3.

I received a lot of new data, so populated the DEFAULT version using the Object Loader tool - this works fine.  However, I am now setting up a MXD which will be used to publish a map service (surfaced on Collector), and want to use the Field Editing version - but the changes I made are obviously not shown in this version.  I don't seem to be able to reconcile and post the changes to the child version.

Is there a way I can do this automatically? Or is the best method to copy the features to the child version, delete those in the Default, and then post and reconcile the changes?


texts: Add two bands to image gives the bands the same name

texts:
I have a sentinel -1 image collection and I have a function that suppose to use focal median in order to apply speckle correction. Each image has 2 bands: VV and VH and the function I have suppose to run this filter on each one and then to add those two bands to the image.

The problem: The added bands don't get the names I chose, as you can see in the following code:

var Sentinel1 = ee.ImageCollection('COPERNICUS/S1_GRD')
        .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
        .filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'))
        .filter(ee.Filter.eq('instrumentMode', 'IW'))
        .select('VH','VV')
        .filterDate('2020-01-01','2020-01-31')
        .filterBounds(geometry);

var clippedSAR= Sentinel1.map(function(im){
  return im.clip(geometry);
});

//print(clippedSAR);

var filterSpeckles = function(img) {
  var VH = img.select('VH') ;//select the VH polarization band
  var VV=img.select('VV');
  var VH_smoothed = VH.focal_median(5,'circle','meters').clip(geometry).rename('VH_Filtered'); //Apply a focal median filter
  var VV_smoothed = VV.focal_median(5,'circle','meters').clip(geometry).rename('VH_Filtered'); //Apply a focal median filter
  return img.addBands(VH_smoothed).addBands(VV_smoothed);
};

var Clipped_Speckles = clippedSAR.map(filterSpeckles);
print(Clipped_Speckles);

When I print this I get those bands:

Even though in the code I gave different names to the bands:

  return img.addBands(VH_smoothed).addBands(VV_smoothed);

why is this happening? how can I be sure it adds the correct bands+ give them the right names?


texts: Show popupson hover doesn't work with WMS layer

texts:
I am using qgis2web on QGIS for making a simple webmap. I imported a WMS layer (road) and i exported the webmap with "show popups hover" checked.But popups it doesn't work.

When i put the layer from my pc (shapefile not WMS) the popups works.

I am something wrong or popups combined with WMS layers It doesn't supported;;


texts: Adding field in vector layer from attribute table in QGIS?

texts:
I have an Attribute Table named 'data.shp' and a Vector Layer 'pop.shp'.

The 'data.shp' has the following fields "ID", "value1", "value2"

The 'pop.shp' has the "ID" field (same with the 'data.shp')

I want to update 'pop.shp' and add the fields "value1", "value2" from 'data.shp' where "ID" (from 'data.shp') is equal to "ID" (from 'pop.shp').

How could I do that in QGIS?


texts: Reproject a NumPy array with affine transform

texts:
I have NumPy array with one band and no CRS. I do have the shape of the image and the coordinates of the bbox, in WGS84. I would like to reproject it based on bounding box and the shape.
I have managed to create an affine matrix as tuple using rasterio but not sure how can I proceed with that.
This is what I did:
rasterio.transform.from_bounds(
    -101.7359960059834,
      20.21904081937658,
    -100.5717967351885,
      20.8312118894487,
    1103,
    2039
)

>>>Affine(0.0010554843796871222, 0.0, -101.7359960059834,
       0.0, -0.0003002310299519955, 20.8312118894487)

So now I have the affine but I don't know how to reproject my array.


texts: QGIS: How to compute mean of raster values without no data?

texts:
I was able to clip a raster against a vector boundary using Clip Raster by Mask layer. I want to compute the mean of the values inside a single raster. I am using Raster Layer Analysis. However it seems to add in the computation the no data values set to 0. How to work around this?

I am using QGIS 3.10 on Windows 10.

The black "outer" values are 0.


texts: Trigger - Calculate area of a layer inside feature by changing the inner layer

texts:
Working with PostGIS database and QGIS, I have two polygons layers, an external with land lots (lotes) and an internal with building area (areaconstruida).
In the "lotes" layer, I need a field in attribute table that calculate the building area inside the land lot.

I'm tying to create a trigger function that re-calculate the building area inside affected land lots every time the buildings layer is changed.

Before run the trigger, I tried some select in SQL to ensure-me that my code was correct.

     select *
            from
            (select lotes.cod_lote, sum(st_area) from public.lotes,

(select lotes.cod_lote, ST_area(st_intersection(lotes.geom,areaconstruida.geom))
            from public.areaconstruida, public.lotes
            WHERE ST_intersects(lotes.geom,
    (select geom from public.areaconstruida where public.areaconstruida.id=3789)))
    AS TAB1
    -- TAB 1 returns, for all features in lotes that intersects an especific building
feature, the area of intersection with all buildings in project

                WHERE lotes.cod_lote=TAB1.cod_lote
                group by lotes.cod_lote) as TAB2
    -- TAB2 returns, for all features in lotes that intersects an especific
building feature, the sum of interception areas with the buildings layer
                WHERE TAB2.cod_lote='00309SE00006'
-- Trying the result of an specific parcel lot

from this code I got expected result, with the total building area inside the land lot tested.

But I'm
unsuccessful trying to activate the trigger

For the Trigger, I replaced the "public.areaconstruida.id=3789" for "NEW.geom", once I want to change the values of land lots affected by new or modified features in the buildings layers
and replaced the WHERE TAB2.cod_lote='00309SE00006' for WHERE TAB2.cod_lote=lotes.cod_lote

Resulting in this code:

CREATE OR REPLACE FUNCTION area_constr2()
RETURNS trigger AS
$BODY$
BEGIN
UPDATE public.lotes
SET lotes.areaconstr = TAB2.sum
from
(select lotes.cod_lote, sum(st_area) from public.lotes, (select lotes.cod_lote, ST_area(st_intersection(lotes.geom,areaconstruida.geom))
from public.areaconstruida, public.lotes, NEW
WHERE ST_intersects(lotes.geom, (select geom from NEW.geom))) AS TAB1
WHERE lotes.cod_lote=TAB1.cod_lote
group by lotes.cod_lote) as TAB2
where TAB2.cod_lote=lotes.cod_lote;
END;

$BODY$
LANGUAGE plpgsql;

CREATE TRIGGER area_constr2 BEFORE INSERT OR UPDATE ON public."areaconstruida"
 FOR EACH ROW EXECUTE PROCEDURE area_constr2();

The query runs normally, but when I try to change something in buildings geometry on QGIS, results in errors.

This question is similar to this topic: (Trigger - Calculate area of a layer inside feature) from yesterday, but now I'm changing the inner layer, and not the outer layer.

Edit: I tried another code to re-calculate the building area of all land lots, instead of only affected ones.
Of course, it is expected that it demands long time to run.

CREATE OR REPLACE FUNCTION area_constr2()
RETURNS trigger AS
$BODY$
BEGIN
UPDATE public.lotes
SET areaconstr = (SELECT T2.SUM
FROM (select TAB1.cod_lote, sum(st_area)
from public.lotes,
(select public.lotes.cod_lote, ST_area(st_intersection(public.lotes.geom,public.areaconstruida.geom))
            from public.areaconstruida, public.lotes
            WHERE ST_intersects(public.lotes.geom, public.areaconstruida.geom))
    AS TAB1
                WHERE lotes.cod_lote=TAB1.cod_lote
group by TAB1.cod_lote) AS T2
WHERE T2.cod_lote=lotes.cod_lote);

RETURN public.lotes;
END;

$BODY$
LANGUAGE plpgsql;

CREATE TRIGGER area_constr2 BEFORE INSERT OR UPDATE ON public."areaconstruida"
 FOR EACH ROW EXECUTE PROCEDURE area_constr2();

It runs in QGIS and takes about 15 seconds, but no changes at all occours.


texts: Plugin/etc. for users adding points to a map on your homepage

texts:
Is there a plugin for if you want to create a map where users (on your homepage) can enter point data?
It is a pretty simple thing, I assume. Just want to have a map on my homepage where any person can add points (maybe even lines) to my map, which is then gathered in a spreadsheet of some sort for me (readable by QGIS). They should be able to choose a category for their point. It's for a map where cyclists can report where they feel unsafe, lack of bike lane etc.


texts: Point z value Return 1.#QNAN while set before

texts:
I am makin a python code for ArcGis 9.3. The code is as following :

rowsLine = gp.updatecursor(line_Feature_Class)
count = 0
rowLine = rowsLine.next()

rowsPoint = gp.InsertCursor(point_Feature_Class)
feat = rowsPoint.NewRow()

while rowLine:

#...

ptObj = gp.CreateObject('Point')
    ptObj.x = x1
    ptObj.y = y1
    ptObj.z = z1
    feat.Shape = ptObj
    rowsPoint.InsertRow(feat)

del rowsPoint

desc = gp.Describe(point_Feature_Class)
shapefieldname = desc.ShapeFieldName

rowsAlti = gp.UpdateCursor(point_Feature_Class)
rowAlti = rowsAlti.Next()

while rowAlti :

    featAlti = rowAlti.GetValue(shapefieldname)
    pntAlti = featAlti.GetPart()
    Alti = pntAlti.z
    rowAlti.SetValue("Altitude", Alti)

    rowsAlti.updateRow(rowAlti)
    rowAlti = rowsAlti.Next()

The problem is that the return value of pntAlti.z is : 1.#QNAN. But if I replace pntAlti.z by pntAlti.x, then it works !

I saw that 1.QNAN means the value is not set. However I set the z value in my code : ptObj.z = z1

Thanks for your help


texts: OpenLayers Scale Line not showing on the map

texts:
I am trying to generate a scale line on my map using ScaleLine but it is not showing on the map.

//Scale
import {defaults as defaultControls, ScaleLine} from 'ol/control';

map = new Map({
      target: 'map-container',
      controls: defaultControls().extend([
        new ScaleLine()
      ]),
      layers: [
        new TileLayer({
          source: new OSM()
        })
      ],
      view: new View({
        projection: 'EPSG:4326',
        center: [114.190750, 22.339914],
        zoom: 12,
        minZoom: 12,
      }),
    });


texts: Generate different number of points for multiple features in shapefile using R

texts:
I have a shapefile that contains roughly 200,000 features, each with a unique Zone Code. Within each of these features I would like to randomly place n points. A separate dataset contains a list of the unique Zone Codes and the number of random points required. The end result I am after is a list of Zone Code's with a set of unique x,y coordinates for however many points are required for each feature.

I need to do this in R and have some code that does the job. However, it is not very quick and takes almost a day to complete. I am looking for any help in speeding up the process.

Here is an example of the code I am using (with a shapefile and dataset that contains 196 features):

require("rgdal")

download.file("https://dl.dropboxusercontent.com/u/27869346/Zone_ZIP.zip", "Zone_ZIP.zip")
unzip("Zone_ZIP.zip")
Zone_Shape <- readOGR(".", "ZONE_CD")

Zone_Data <- read.csv("https://dl.dropboxusercontent.com/u/27869346/ZONE_DATA.csv",header=TRUE)

Points_All <- NULL
for(i in 1:nrow(Zone_Data)){
if(Zone_Data[i,2]>0){
Zone_Shape_Single <- Zone_Shape[Zone_Shape@data$ZONE_CD %in% Zone_Data[i,1], ]
Points <- spsample(Zone_Shape_Single, Zone_Data[i,2], type="random", iter=10)
Points_DF <- cbind(Zone_Data[i,1],as.data.frame(Points))
Points_All <- rbind(Points_All, Points_DF)
}
}


texts:
I want to display and work with the data from http://dds.cr.usgs.gov/srtm/version1/Eurasia/
in SAGA-GIS but am having difficulties knowing which module I should use to load the data.

As far as I know hgt files are raw binary data, so this should work with the Import Binary Raw/data but by simply loading the data into the import window, nothing happens meaning I have to define the correct parameters, but how?


texts: Organizing Arcmap folders for portability

texts:
I am doing some basic constraint mapping in ArcMap 10.1 and need some pointers on organizing my project for portability. My project TOC in ArcMap contains shapefiles in the following structure:

Project Details
-Project Lands

Features
-POI's, roads, railways etc

Constraints
-Setbacks from features

Typically the folder structure for our Arcmap and data looks like this:

Arcmap file path:
Client -> Project -> ArcMap -> Date -> mxd file

Shape file paths:
Client -> Project -> Data -> Data Type i.e. Project detail, feature, constraint -> Date -> shp file

Often we will package up the Arcmap and shp files to send to the client. What is the easiest way to do this without breaking links in Arcmap? With the folder structure the way that it is, I am unable to only pull the Arcmap -> Date -> mxd file and Data -> Date -> shp files.

I have played with exporting/importing to/from file geodatabases in ArcCatalog but don't think this is the solution other than it being just one file rather than a bunch of different shape files to send to the client.

The other solution that I think should work would be making the paths relative but it seems whenever I change a folder name it breaks the data source and I end up having to reconnect everything which is error prone.


texts: ArcGIS JSAPI draw toolbar tooltips not displayed in IE10

texts:
I am experiencing an issue where the draw toolbar tooltips are not displayed in Internet Explorer 10 when using the ArcGIS Javascript API.

Try the following example in Internet Explorer 10 and Google Chrome (draw a polygon and see what I mean)
http://developers.arcgis.com/en/javascript/samples/graphics_add/

Has anyone else experienced this issue and is there any way to get the tooltips to display in Internet Explorer 10?


texts: QGIS2web setting condition for the attachment widget

texts:
I have a layer, where I need to set an attachment as a widget.
Unfortunately only 4 elements have the attachment. Other don't. In turn, when exporting the map I have the path icon, which I don't want.

If I choose the widget "Hidden", then for all elements the path won't be visible. I don't want that, because I have 4 photos of 4 individual objects.
I tried to set an expression for it:
where at the "Default Path" label I provided:
       if("Path"='NA',"","Path")

but this icon is still visible when extracting the map.
Does anyone know how to solve it?