texts: Counting % overlays of shapefile polygons with buffers

texts:
I have previously asked question here but it was closed because I did not react fast enough twice already and I have put in the comments why those questions do not answer my problem. Nevertheless, my problem is still the same: I have a set of points (meteorological stations) with 5km buffers around them. Additionally, I have a set of polygons (municipalities) in which those stations are located. Each station buffer can be over 1,2,3... polygons and those can be completely different polygons.

I would like to have data on how much of the buffer is filled by each polygon but I was not able to find a function or a way to do it. So for example for point A I would have that in the buffer of 5km there is 80% of polygon B and 20% of polygon C as the output, so 4 values, two for IDs of the polygons and 2 for the percentage overlap (it could be only 2 values for other buffer or 6 etc. etc.).

Main difference to the other similar question provided is that I have 2 vector layers (Buffer with 500 buffers and Municipalities a shape file with 500 different polygons). Previous questions were analyzing this question in situations in which I know precisely the polygons and buffers analyzed, at least that is how I understood them, so after trying to use those approaches on my dataset I was not able to replicate them.

I have also commented on the approaches posted in the previous question but did not get any response to those sadly. Bottom line being I understand those approaches and was able to replicate them in a simple setting but when I do not have specific buffers/polygons I do not know how to execute them. Another way of approaching this problem was to use join attributes by location -one-to-many and then getting averages of the duplicate results, only problem is that then I take average of polygons instead of their weights corresponding to % of intersection with the buffer. In my case I do not know what polygons of the shapefile are going to intersect what buffers, or even if the buffer is going to be intersected by 1,2 or more polygons


texts: SVG iconed features not rendered correctly during zoom

texts:
Im using Geoserver with PostGIS database to put points to map using GeoExplorer. The points are styled to use SVG icon with fixed size (uom metre, size) so they scale together with the map as user zooms in and out. For now, there are no filters or any other issues. Style code :

<sld:PointSymbolizer uom="http://www.opengeospatial.org/se/units/metre">
        <sld:Graphic>
          <sld:ExternalGraphic>
            <sld:OnlineResource xmlns:xlink="http://www.w3.org/1999/xlink" xlink:type="simple" xlink:href="load_icon.php?id=${param}"/>
            <sld:Format>image/svg+xml</sld:Format>
          </sld:ExternalGraphic>
          <sld:Size>80000</sld:Size>
        </sld:Graphic>
      </sld:PointSymbolizer>
      <sld:TextSymbolizer/>

load_icon.php is returning the SVG code for given param :

<svg xmlns="http://www.w3.org/2000/svg" version="1.1" width="220" height="110">
<g transform="scale(0.5)">
<g transform="translate(90,50)">
<rect x="1.5" y="1.5" height="103" width="167" fill="#fff" stroke="red" stroke-width="3"/>
<path d="M3,3 l164,100 m0,-100 l-164,100" stroke="red" stroke-width="3" fill="none"/>
</g>
<text x="173" y="49" text-anchor="middle" font-family="Verdana" font-size="25" fill="red">XXXXX< /text>
<text x="270" y="153" font-family="Verdana" font-size="25" fill="red">G< /text>
</g>
</svg>

I use the plain GeoExplorer to view the icon on the map, but same issue is when connecting the layer as WMS to Leaflet and other libraries.

The icon is showing ok in some resolution :

but gets cropped as user is zooming in :

and even more cropped with further zoom :

with tiles that do not contain the part of the icon. There is caching turned off for this layer, no help.

What I am doing wrong here ?


texts:
I want a button "export" as I clicked on it the shapefile should be downloaded along with it's all files shp, shx, prj, dbf as zip folder. My project is GIS based in which I have two main tasks:

Import shapefile: shapefile to geojson, that is i am doing with the help of leafLet.
Export shapefile: in this after draw/editing on map/layer/shapefile the duplicate shapefile should be created.

How I can export the copy of that particular shapefile?


texts: Choose which MULTIPLE attributes should be joined -QGIS Graphical Modeler

texts:
I need to select which MULTIPLE attributes to join in a QGIS model. It allows me to do this in the interface (Properties-->Join) but not in the QGIS Model version of Join by Attribute.

So I want the match to be by LOTPLAN and the then some attributes joined and others not.

Join tool in modeler


texts: Problem showing WMS with EPSG:3301 projection using Leaflet

texts:
I am trying to show wms layer EPSG:3301 using leaflet, but have no succes for a long time.
http://kaart.maaamet.ee/wms/alus?service=WMS&version=1.3.0&request=GetCapabilities
Here is my code:

function initialize() {

var crs = L.CRS.proj4js('EPSG:3301',
     '+proj=lcc +lat_1=59.33333333333334 +lat_2=58 +lat_0=57.51755393055556 +lon_0=24 +x_0=500000 +y_0=6375000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs',
    new L.Transformation(1, -40500, -1, 7017000)
);

var wmsMap= L.tileLayer.wms('http://kaart.maaamet.ee/wms/alus', {
    layers: 'MA-ALUS',
    format: 'image/png',
    maxZoom: 14,
    minZoom: 3,
    continuousWorld: true,
    transparent:true
});

var map = L.map('map', {
        layers:[wmsMap],
        continuousWorld: true,
        center: new L.LatLng(58.66, 25.05),
        zoom: 2,
        maxZoom: 13,
        minZoom: 1,
        maxBounds: new L.LatLngBounds([
            [56.42, 20.87],
            [60.9, 29.23]
        ]),
        crs: crs
    }
);

Please anybody check my code to see what's wrong with it?? I see map area, see attribution and layer switcher if I add them, but not a map.


texts: How to get Set-Cookie header from ArcGIS REST services?

texts:
When I call an ArcGIS REST service, I get Cache-Control, Etag, Content-Encoding, Content-Type, Transfer-Encoding and Date headers.

But I don't get Set-Cookies header.

Is there a way to get this header? Is there any configuration to make on server side?


texts: How to speed up this simple MySQL points in the box query

texts:
MySQL doesn't have nearest neighbor search so. So I try to come up with a "distance" that's reasonable so the number of returned query is around 20-200.

That way innodb doesn't have to compute distance for like 30k rows if we only need 20. The latitude and longitude is also stored in a myisam table called tableauxiliary

Well, even after that, queries still take 13 seconds and sometimes 30 seconds.

SELECT TB.ID ,
  TB.Latitude,
  TB.Longitude,
  111151.29341326*SQRT(pow(-6.185-TB.Latitude,2)+pow(106.773-TB.Longitude,2)*0.98839228980165) AS Distance
FROM
  `tablebusiness` AS TB
  join tableauxiliary as TA on TA.BusinessID=TB.ID
WHERE
  MBRContains(
        GeomFromText (
            'MULTIPOINT(-6.1928749092968 106.7651250907,-6.1771250907032 106.7808749093)'
            ),
            TA.Latlong)
  AND MATCH (TA.FullTextSearch) AGAINST ('re*' IN BOOLEAN MODE)
AND TA.Prominent >15
ORDER BY
  Distance

Number of matches are 32 by the way.

Showing rows 0 - 29 ( 32 total, Query took 16.1845 sec)

Also an even simpler query where we do not compute distance also take 16 seconds

SELECT TB.ID
FROM
  `tablebusiness` AS TB
  join tableauxiliary as TA on TA.BusinessID=TB.ID
WHERE
  MBRContains(
        GeomFromText (
            'MULTIPOINT(-6.1928749092968 106.7651250907,-6.1771250907032 106.7808749093)'
            ),
            TA.Latlong)
  AND MATCH (TA.FullTextSearch) AGAINST ('res*' IN BOOLEAN MODE)
AND TA.Prominent >15

I can understand if mysql cannot do nearest neighbor search. What I am doing now is a simple search points in a box here. What am I missing?

The explanation for the query is:

id  select_type     table   type    possible_keys   key     key_len     ref     rows    Extra
1   SIMPLE  TA  fulltext    PRIMARY,Prominent,LatLong,LatLong_2,FullTextSearch  FullTextSearch  0       1   Using where
1   SIMPLE  TB  eq_ref  PRIMARY     PRIMARY     752     gensurv_isikota.TA.BusinessID   1   Using index

It basically says only one row is computed for each of the process.

number of data is 1.6 million, but number of data in a box is only 30


texts: It is necessary that the labels are evenly spaced relative to each other and do not clump together. Where can I specify the parameter in the settings: the distance between the labels?


texts: Increment a counter in the catalog_apply function (lidR)

texts:
following up on this question: Capture and parse output of lascheck from lidR... (while lidR 3 is in development, I thought I'd continue with a workaround.)

I'm using catalog_apply to perform lascheck on individual las files in the catalog. I'm getting the right output for warnings/errors, but I can't identify the name of the las file (chunk) because although the warnings/errors are there in the output of lascheck, the name isn't.

However, a list of names is available in the LASCatalog. I would like to extract the corresponding name, but I can't increment a variable to count (for reasons of scope I think). I.e. capture.output(ctg@data$filename[x]) where x is position of the name in the list.

E.g. the following code does return the warnings/errors but x gets stuck at 1.

library(lidR)

ctg <- readLAScatalog("D:/Projects2020/Sherbrooke/") #, progress = TRUE, list.files())

x <- 0

errwrn <- function(chunk)
{
  las = readLAS(chunk)
  if (is.empty(las)) return(NULL)

  las@header@PHB$`Point Data Format ID` <- 12L # Introduce errors.

  ew <- capture.output(lascheck(las))
  err_pos  = grep("+2717", ew)
  errs = ew[err_pos]
  warn_pos = grep("+26A0", ew)
  warns = ew[warn_pos]

  x <- x + 1
  lasName <- capture.output(ctg@data$filename[x])
  ews = c(x, lasName, errs, warns)

  return (ews)
}

lascheckList <- catalog_apply(ctg, errwrn)

Can anyone see how I could increment x?  Or maybe I should use a loop to go through each file in the directory, bypassing the use of a catalog?


texts: How to get line and polygon points from OpenLayers 3 drawing

texts:
I am working on drawing points/segments/polygons on a map that uses OpenLayers 3. The first part I have done draws points, lines and polygons.  Now I am trying to retrieve the drawn coordinates in order to build them for other sessions.

I have been able to getting just coordinate for the point/image, but can't get coordinates for segment/polygon.  Any ideas or help?  I've included a JSfiddle below

https://jsfiddle.net/Lg8w71kv


texts: How to output multiple raster in QGIS processing scripts

texts:
I am writing a script in R for running in QGIS processing. The code is a oil spill classification algorithm and depending on the number of polygons which user inputs in the shapefile, output raster layers are to be produced, each containing the oil/water classification plot corresponding to each polygon in the shapefile. Since I do not know in advance how many polygons user is going to input, I cannot specify in the processing header like:

##out_plot_raster_1 = output raster
##out_plot_raster_2 = output raster
##out_plot_raster_3 = output raster
.
.
.

QGIS processing scripts support multiple raster input like:

## input_many_rasters = multiple raster

so it is logically inferred that same must be true for output also, but on trying:

## out_many_rasters = output multiple raster

I get an error and the script won't run at all.

Now, it is really important for me to be able to output multiple rasters at once so that user can get all the friendly options like "saving them to desired location" and "opening them in QGIS after the script run" which are automatically provided by the processing framework of QGIS. How do I do this? i.e. I want user to be able to enter, before the script run, the location to store all generated raster outputs and if he wants, the rasters must be opened in QGIS automatically after the execution is complete.
Please suggest me a way to do this.


texts: Output field name when using gdal.Polygonize() in python

texts:
I am trying to Polygonize one raster image that previously I convert to 1s and 0s.

If I do it from QGIS, there is the option to select the output field name, which stores as attributes the value of the pixel, and let me iterate over all the features afterwards and store only the ones with value 1. (Since apparently it's impossible to polygonize only pixels with a specific value)

The problem is, that I can't see in the documentation of gdal.Polygonize() (http://www.gdal.org/gdal__alg_8h.html#a3f522a9035d3512b5d414fb4752671b1) any option to select the output file, so I get a bunch of polygons numerated from 0 to end.

Is there any way of keep the values of the pixels as attributes?


texts:
I hope this is the right place to ask, I was referred here from a topic on SO. The same topic pointed me towards using GeoServer to create maps with overlaid tracks on. I have unpacked the .war file onto my tomcat server but am getting a number of errors as listed in my catalina.out file. Having never used GeoServer before it doesn't really mean much to me at all.

Does it make sense to anyone else:

http://pastebin.com/JNkWdVER

Thanks in advance for any light you can shed and help me get this working.


texts: How to Calculate the average value near a point

texts: I am new to ArcGIS ArcMap. I want to calculate the average value of one point ( layer flowline ). This could be done using the value of the point and the values of the four neighbor points ( raster ), eventually to calculate the mean value of them. Is any way to calculate this mean value for each point of a flowline-layer that have totally 10 points,


texts:
Using ArcGLobe I'm attempting to display elevation data.

I'm using Digital Terrain Elevation Data (DTED) data and trying to load it in the following fashion.

IRasterLayer rasterLayer = new RasterLayerClass();
rasterLayer.CreateFromFilePath(myFilePath);
Globe.AddLayerType(rasterLayer , esriGlobeLayerType.esriGlobeLayerTypeElevation , true);

Nothing shows up on the globe.  Thoughts?  Ideas?


texts: Import csv data to local nominatim/photon installation

texts: I have some additional address data (.csv file) that I would like to import to my photon/nominatim installation. Is there some way to do this? I cannot upload the dataset to OSM because it is not open data.


texts: Which is the correct distribution to use?

texts: I wanted to know that if the pre-build js of OpenLayers-3 available at http://openlayers.org/en/v3.2.1/build/ol.js is self-sufficient or does it need to use the internet for any of its functioning?


texts:
I am using windows10 64 bits, my os.name is 'nt', I used pip to  install but its complaining about the spatialindex_c.dll

 File "C:\Users\Bachir\AppData\Local\Temp\pip-build-td64lrth\rtree\rtree\core.py", line 101, in <module>
   raise OSError("could not find or load spatialindex_c.dll")
 OSError: could not find or load spatialindex_c.dll

I downloaded rtree-0.8.2 and then run python setup.py install and i have the same message

PS C:\Users\Bachir\documents\Python Scripts\Rtree-0.8.2> python setup.py install
Traceback (most recent call last):
File "setup.py", line 4, in <module>
import rtree
File "C:\Users\Bachir\documents\Python Scripts\Rtree-0.8.2\rtree\__init__.py", line 1, in <module>
from .index import Rtree
File "C:\Users\Bachir\documents\Python Scripts\Rtree-0.8.2\rtree\index.py", line 6, in <module>
from . import core
File "C:\Users\Bachir\documents\Python Scripts\Rtree-0.8.2\rtree\core.py", line 101, in <module>
raise OSError("could not find or load spatialindex_c.dll")
OSError: could not find or load spatialindex_c.dll
PS C:\Users\Bachir\documents\Python Scripts\Rtree-0.8.2>

I donwloaded the spatialindex dll files ' libspatialindex-1.8.1-win-msvc-2010-x64-x32.zip ' . this file contain both 32  and 64 bits ,  unzip and put in the installation folder, when installing using python setuo.py install it still complaining about the spatialindex dll file


texts: How to specify HGT file name in gdal_translate command?

texts:
gdal_translate command is generating a HGT file with user defined name but its giving an error as `N5464347W74.hgt' not recognized as a supported file format. Below is the command and the entire output.

**gdal_translate -of SRTMHGT 1e9815c1-b74e-4a75-8a98-844d730707e7.tif user_def.hgt**

Input file size is 1201, 1201
Warning 1: The corner coordinates of the source are not properly aligned on plain latitude/longitude boundaries.
Warning 1: Expected output filename is N48W002.HGT.
0...10...20...30...40...50...60...70...80...90...100 - done.
ERROR 4: `N5464347W74.hgt' not recognized as a supported file format.

But if we provide the file name suggested by the command then it runs error free. Not sure how can we extract the file name. following is the gdalinfo for the tif file.

**gdalinfo 1e9815c1-b74e-4a75-8a98-844d730707e7.tif**

Driver: GTiff/GeoTIFF
Files: 1e9815c1-b74e-4a75-8a98-844d730707e7.tif
Size is 1201, 1201
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-1.686258165938839,48.111195526365400)
Pixel Size = (0.000002416422364,-0.000001681829966)
Metadata:
  AREA_OR_POINT=Point
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  -1.6862582,  48.1111955) (  1d41'10.53"W, 48d 6'40.30"N)
Lower Left  (  -1.6862582,  48.1091756) (  1d41'10.53"W, 48d 6'33.03"N)
Upper Right (  -1.6833560,  48.1111955) (  1d41' 0.08"W, 48d 6'40.30"N)
Lower Right (  -1.6833560,  48.1091756) (  1d41' 0.08"W, 48d 6'33.03"N)
Center      (  -1.6848071,  48.1101856) (  1d41' 5.31"W, 48d 6'36.67"N)
Band 1 Block=1201x1 Type=Float32, ColorInterp=Gray
  NoData Value=0
  Unit Type: metre

Any suggestions???


texts: Resampling 13 arc second DEM to 1/3 arc second

texts: I need to download DEMs for almost half of the U.S. but not all the locations have 1/3 arc second DEMs but the have 13 arc second. Can anyone help me how to resample 13 arc second DEM to 1/3 arc second?


texts: Reloading OpenLayers 5 WMS automatically with new data

texts:
I am making a weather application. It is almost finished now, but I just have one problem left. I want a WMS-layer to automatically reload every few seconds. I thought I had the solution, but it seems that it doens't quite work. Below (some of) my code.

It seems to redraw the WMS, but not with changed data. Only when the datetime or weathertype CHANGES, then the WMS-gets new data and that is shown on the map. But I want it to keep reloading and show the new data for that datetime. (Because the datetime only changes once in a while, but the data keeps rolling in.)

I hope I made my question clear.

I searched and found some solutions, but they didn't work for me. Besides most solutions I found were for OL2 or OL3 and I am working with OL5.

Please help me.


texts: Point in Polygon operation for 100G of data

texts:
I have lat/lon point locations for over a billion points, all based in the US. For each point I would like to generate a dataset that links them to their county and to their MSA according to US Census Shapefiles.

My question is how to do this efficiently? I have access to a research computing server, postgis, shapely and statistical software.


texts: Open source solution to convert terrain data to collada?

texts:
I am exploring how to convert terrain information (DTED, SRTM, etc) to a 3D mesh file.  Collada is preferrable as it is very compatible.

I have read this: Converting dem and geotiff to collada model with texture?

and followed the process of exporting terrain in QGIS using Qgis2threejs, selecting the File Export template, but no .dae file object is created - only several .js files.

I get that broad questions are frowned upon on SE, but I'm not sure how else to phrase it.  Does anyone have a recommended method to convert a terrain raster to .dae or other widely compatible 3D file format?


texts: PgRouting coverage area with node restrictions

texts:
I have a water network in which there exist valves at nodes between two edges (pipes). When the valve is open water can flow, if not, it can't.

How can I calculate the catchment area from a start point in the network(e.g. a break/pipe burst) that takes into account not just pipe/edge costs (e.g. pipe diameter), but also the cost at the stop valves/nodes (e.g. high cost when closed, no cost when open)?


texts: Include polygons that don't contain line segments in Summarize

texts:
I would like to get a summary of road lengths per township, but also include townships that have no roads.

Using the first answer from this question Calculate line distances inside of polygons using ArcGIS for Desktop?

...how would I include the polygons (in my case, townships) that don't have any polylines (roads)? The Summarize spreadsheet that outputs only includes that do.


texts: Adding column to attribute table indicating polygon each data point is in using ArcGIS Pro

texts:
I am working on a project linking ideas in Geology and Ecology.
I currently have a map showing the geology of an area (polygon layer), with a points layer on top indicating the locations of data collection.
I have added an empty column to the attributes table for the points layer (geology_where_captured) and would like to know if there's a way - using arcpy -  to fill in this column with the name or id of the geology that each data point is sitting within?
I have looked into the SummarizeWithin function but when I try to use this (I am using ArcMap 10.7.1 on my home computer) I get an error saying the object has no attribute?
I'm pretty new to Python.
I wrote this small bit of code following advice from https://pro.arcgis.com/en/pro-app/tool-reference/geoanalytics-desktop/summarize-within.htm
and
https://pro.arcgis.com/en/pro-app/tool-reference/analysis/summarize-within.htm
import arcpy
arcpy.env.workspace = 'C:\Users\User\Documents\University\4th Year Project\Data\mydata_arcgisOnline/XYLEMUR_DATA_AMBATOVY2.dbf'
summarizedLayer = "XYLEMUR_DATA_AMBATOVY2"
summaryPolys = "Geology"
arcpy.SummarizeWithin_analysis(summarizedLayer, summaryPolys, 'geology_where_captured')

Runtime error Traceback (most recent call last):   File "",
line 1, in  AttributeError: 'module' object has no attribute
'SummarizeWithin_analysis'

No joy! (I'm new to this)


texts:
I want to clip a polygon using Qgis "clip" tool. The problem is that the resulting polygon presents areas that where not clipped off.
The first screenshot shows the layer I want to clip using the polygon on the second screenshot. In the resulting layer (3rd screenshot) you can see the areas of the clipping layer that were not clipped.

How can this be fixed?


texts: How to add Lat/long columns to a point layer which is not in WGS84?

texts:
I need to extract the Lat/long coordinates of the centroids of districts.
My QGIS project SCR and the different layers have all been set on "WGS 84 / UTM zone 18N".

I've created the centroids using the usual tool.
I could then extracted the coordinates using the "Vector > Geometry Tools > Export/Add Geometry Columns", or saving the layer in CSV typing "GEOMETRY=AS_XY" in the "Layers" part of the "OGR creation options".

Nevertheless, both methods output Cartesian coordinates, whereas I need Lat/long coordinates.


texts:
I have a series of raster map images that I would like to use in QGIS and overlay some simple GPX routes over the top to produce a custom map image.
This is all quite new to me. I am struggling to get the raster image and GPX layers to display on the same map. My suspicion is that the two layers are using different projections or coordinate systems (excuse me if this the incorrect terminology).
When both layers are displayed in QGIS, they are not even visible in the same display. The coordinates displayed for each are wildly different and are not of the same magnitude.
I'd like to know how to covert the images I have so that can use them with GPX files. For information, these are Australian NSW Lands Department maps as provided with the TopoView 2006 DVD product.
gdalinfo in the raster image reports:
./gdalinfo ~/Desktop/89313s.ecw
Driver: ECW/ERDAS Compressed Wavelets (SDK 3.x)
Files: /Users/jon/Desktop/89313s.ecw
Size is 6168, 3808
Coordinate System is:
PROJCS["MGA56",
    GEOGCS["GEOCENTRIC DATUM of AUSTRALIA",
      DATUM["GDA94",
          SPHEROID["GRS80",6378137,298.257222101]],
      PRIMEM["Greenwich",0],
      UNIT["degree",0.0174532925199433]],
  PROJECTION["Transverse_Mercator"],
  PARAMETER["latitude_of_origin",0],
  PARAMETER["central_meridian",153],
  PARAMETER["scale_factor",0.9996],
  PARAMETER["false_easting",500000],
  PARAMETER["false_northing",10000000],
  UNIT["Meter",1]]
Origin = (220464.635367600014433,6304336.451381999999285)
Pixel Size = (4.000000000000000,-4.000000000000000)
Corner Coordinates:
Upper Left  (  220464.635, 6304336.451) (149d59'44.79"E, 33d21'50.07"S)
Lower Left  (  220464.635, 6289104.451) (149d59'27.75"E, 33d30' 3.98"S)
Upper Right (  245136.635, 6304336.451) (150d15'38.36"E, 33d22'12.14"S)
Lower Right (  245136.635, 6289104.451) (150d15'22.81"E, 33d30'26.17"S)
Center      (  232800.635, 6296720.451) (150d 7'33.40"E, 33d26' 8.35"S)
Band 1 Block=6168x1 Type=Byte, ColorInterp=Red
Overviews: 3084x1904, 1542x952, 771x476, 385x238
Band 2 Block=6168x1 Type=Byte, ColorInterp=Green
Overviews: 3084x1904, 1542x952, 771x476, 385x238
Band 3 Block=6168x1 Type=Byte, ColorInterp=Blue
Overviews: 3084x1904, 1542x952, 771x476, 385x238

Extract from GPX is as follows:
<?xml version="1.0" encoding="UTF-8"?>
<gpx version="1.1" creator="Garmin Connect"
  xsi:schemaLocation="http://www.topografix.com/GPX/1/1        http://www.topografix.com/GPX/1/1/gpx.xsd http://www.garmin.com/xmlschemas/GpxExtensions/v3 http://www.garmin.com/xmlschemas/GpxExtensionsv3.xsd http://www.garmin.com/xmlschemas/TrackPointExtension/v1 http://www.garmin.com/xmlschemas/TrackPointExtensionv1.xsd"
  xmlns="http://www.topografix.com/GPX/1/1"
  xmlns:gpxtpx="http://www.garmin.com/xmlschemas/TrackPointExtension/v1"
  xmlns:gpxx="http://www.garmin.com/xmlschemas/GpxExtensions/v3"  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <metadata>
    <link href="connect.garmin.com">
      <text>Garmin Connect</text>
    </link>
    <time>2010-08-30T21:03:24.000Z</time>
  </metadata>
  <trk>
<name>Canberra to 'Castle Day 3</name>
    <trkseg>
      <trkpt lon="150.05063133314252" lat="-33.67591420188546">
        <ele>1030.0</ele>
        <time>2010-08-30T21:03:24.000Z</time>
        <extensions>
          <gpxtpx:TrackPointExtension>
        <gpxtpx:atemp>2.0</gpxtpx:atemp>
          </gpxtpx:TrackPointExtension>
        </extensions>
      </trkpt>
      <trkpt lon="150.05062630400062" lat="-33.67588419467211">


texts: How can i display the color of my labels into the content table?

texts:
The thing is that i need to display labels based on colors but without the symbol. In this case  i'm working with points in Arcgis 10 and the maps needs to display only labels with color-based categories and in the content table the layers display only the symbol.

My question is if its possible to, somehow, the color of my labels to be displayed in the content table?

I hope you can help me on this.


texts:
It may seem like a simple question but I haven't found a solution.
I have a structure of directories with many shapefiles inside. Is there a method to open all these shapefiles at the same time with ArcMap?

The function "add data" forces me to do it for every directory separately.


texts: Google earth engine: function to apply code to sequential dates

texts:
0
I have a code which creates a value for several different variables for a specific time period (usually over 1 month) and exports these as a csv to my google drive. You can find a link to the code here.
At the moment, if I want to create the code for different dates I just do this manually, by altering the startDate and endDate line in my code.
// select dates
 var startDate = ee.Date('2016-01-01');
 var endDate = ee.Date('2016-02-01');

 var startyear = 2016; //
 var endyear = 2016;

 Export.table.toDrive({
 collection: Dat,
 description: 'Dat',
 selectors: cols,
 fileFormat: 'CSV'
});

The above would then create the data values for January 2016, and export it as a task which I then run in the google earth engine console.
I'm looking to automate this so I can get monthly values over a number of years, or values for specific months within each year, with out having to do this manually. So rather than changing the startDate, endDate, startyear and endyear manually, it would run the code and start again.
I'm aware in some coding systems you can do this with for loops, but I'm not familiar with them at all, and in google earth engine these are not recommended.
Is it simple enough to set up a function that calculates the score for four months (e.g. for Jan, April, July and October) for a set number of years (say 2003-2019), so it would run the score for Jan 2003, then for April 2003, July 2003, etc?


texts: Know the shadows of a house before building it

texts:
Imagine a house/building, seen from outside, at sunrise. As the the sun moves in the horizon, a) the shadow of the house is displaced/deformed, and b) some faces of the house have more/less sunlight in the afternoon than in the morning.

I need to predict those changes along the day and along the year. Besides the external 3D shape of the house, this only depends on space & time. Space: the latitude, longitude and orientation of the house/building. Time: the hour, day, month of interest.

I have the exterior measurements of the house. At this point I treat the house as a solid object (ie I don't care about the light inside the house; I also don't care about cases where light which passes thorough the building eg two opposite windows).

I only need to predict two interrelated things (along the day/year) for planning and visualization purposes:

1) the shadow of the house

2) the faces of the house which will be receiving light

With (1) I can predict how the adjacent buildings/parks/spaces would be affected by sunlight decrease and with (2) I can plan to put more/less and/or smaller/larger windows in different places of each wall of the house.

I know that many years ago this was possible with a plugin for Google Sketchup, but for different reasons I would like not to use this option but a free application (if FOSS, even better). Given the recent surge of 3D modeling/design programs I thought maybe there are other programs which are capable of doing this.


texts: Lining up cells from two raster datasets

texts:
I am using ArcGIS 10.4.1 to convert a polygon feature class to a raster dataset and I would like the cells to line up with an existing raster layer.  I am running the feature to raster tool and using the existing raster layer as a reference for the cell size and as the snap raster.  However, my output raster cells are still offset as seen here:

Any other suggestions for getting the cells to line up?  I have double checked that my existing raster and feature class are both in the same projection.

Update - I also tried the suggestions in the answer listed here: https://geonet.esri.com/thread/98483

I turned off background processing and used data from straight from the catalog and am receiving the same result.


texts: Using API.query() to get the query give me an error

texts:
I have a problem when implementing the API.query () to get the following query it give me an error.
I saw  this solution from stack exhange and I saw that its recommended to put this (._;>;);
line at the end of the query but still didn't work
res=api.get("""
([bbox:40.72797,-73.90540,40.73293,-73.89534];
  node["highway"="secondary"];
  way["highway"="secondary"];
  relation["highway"="secondary"];

  node["highway"="tertiary"];
  way["highway"="tertiary"];
  relation["highway"="tertiary"];

  node["highway"="unclassified"];
  way["highway"="unclassified"];
  relation["highway"="unclassified"];

  node["highway"="residential"];
  way["highway"="residential"];
  relation["highway"="residential"];
);
(._;>;);
""")

when I run the code it gives me this error
ERROR:root:An unexpected error occurred while tokenizing input
The following traceback may be corrupted or invalid
The error message is: ('EOF in multi-line string', (1, 0))


texts: How to display only certain features in a layer in Global Mapper?

texts:
Is there a way to restrict the records displayed on the screen from a shapefile or other layer type in Global Mapper?

In ArcGIS this is called a definition query. For example if I have a text attribute called type I want a query that says: "TYPE" = 'SIDE TRACK'

If I double click the layer in the Overlay Control Center, choose Line Styles and Apply Styling Based on Attribute/Name Values, I can then apply a Null style. This method works, but there are lots of different attribute combinations this would be quite tedious. I'm looking for a way to apply a SQL type query.


texts: How to add vertex on the center of a polyline object?

texts:
I've got an arcpy.Polyline() object and I have to insert a new vertex in it, which has to be placed on a center of a line.
Right now I can calculate the center point coordinates but how to add it to an existing Polyline object in the right place of points order?
The existing vertices order should not be changed, so the output line has to have previous vertices plus the new one somehere in the middle.

with arcpy.da.UpdateCursor(line_in,'SHAPE@') as uc:
    for row in uc:
        Midpoint = row[0].positionAlongLine(0.50,True).firstPoint
        mid_p = arcpy.Point(Midpoint.X,Midpoint.Y)

ArcGIS 10.2 has no segmentAlongLine method so is it possible to update line's vertices here?

Due to the suggestion that this question has been already asked I would like to clarify that I don't need a new layer with points. The goal is to add a vertex to existing line.


texts: Blocky Raster from TIN. How to increase resolution?

texts:
So i am trying to create a seamless DEM using a DTM created in pix4d with a 5cm resolution (Raster).  I have then collected RTK bathymetric data which i have turned into a TIN below
To create a seamless DEM i have to merge two raster DEMs together.....as i understand it. When i convert the river TIN into raster, the output is sub-standard. The cell size is 2.8m which is why i think it apears so blocky but when i try to resample it, there is no change. As you can see from above, this does not represent the river very well, there is overlap onto dry land and exposed water showing. How can i get this to look more like the DTM quality in the first image which is a lot smoother?


texts: New Network Dataset Error: "There are no feature classes in the network container that can be added as feature sources"

texts:
I have three feature classes in a database: a line feature class with roads, a point feature class with road intersections, and another point feature class that I would like to use as starting points for a service area calculation using Network Analyst in ArcMap 10.3.

When I try to create a new Network Dataset, I get an error that reads:

  There are no feature classes in the network container that can be added as feature sources.

What exactly is a feature source, and how can I create a network data set using these shape files?

What qualities must a feature class possess in order to be added as a feature source?


texts: Does anyone know how to get the QGIS Web Client measure tool to show units in feet and miles?

texts: Does anyone know how to change the QGIS Web Client measure tool (distance and area) so that it shows feet or miles rather than meters?


texts: How is data formatted in the Environment Agency's LiDAR Composite DSM 1m files

texts:
I haven't been able to find how the level data is formatted in the files and hopefully somebody knows.
The files have an asc identifier and data is typically formatted as follows (from file tq2515_DSM_1M.asc):
   ncols 1000
   nrows 1000
   xllcorner    525000
   yllcorner    115000
   cellsize     1
   NODATA_value -9999
   33.182 33.232 33.322 33.392   etc (1000 values ?)
       |
       |
       V  (1000 rows)

My best assumption is they contain ground level data at 1m intervals for a 1000m grid. For instance in the file above the bottom left OS grid reference is 525000, 115000 and then the level data for each point is listed from left to right and top to bottom ie data starts at GR 525000, 116000 and finishes at 526000, 115000.  There are a thousand lines of data in the file and (I assume) 1000 values in each line.
Can anybody confirm as I have checked the levels at each corner against the contours on the 1:25,000 OS maps and (in one instance) the level was not within the contour interval?
I am not a computer geek and only starting to learn GIS, but this seemed an ideal way to refine the ground levels from the 5m contour interval shown on OS mapping (UK).


texts: Calculate conditional in line feature using an expression

texts:
I'm trying to calculate a conditional in a line feature but I can't see were my sintax is wrong. Could you guys help me?

def Reclass(TIPO_PAV):
    if TIPO_PAV=='TERRA':
        return !shape.length@kilometers! /40*60
    if TIPO_PAV=='ASFALTO':
        return  !shape.length@kilometers!
    else:
        return 0


texts: QGIS "Statistics by Categories" tool greyed out at 3.8.3

texts:
I have both 2.18.28 and 3.8.3 loaded on my computer.  I have a Geopackage containing a line layer of hiking trails.  The attribute table contains two fields:

surface_type (text)
length_miles (decimal)

I am trying to sum length_miles by surface_type, using the Statistics by Categories tool.

Using 2.18.28 the tool works perfectly.

However, at 3.8.3 the Field(s) with categories option is greyed out (please see highlighted screenshot).  This is where I would normally select the surface_type field.  Now that it is greyed out I cannot calculate miles by surface_type.

The Help section offers no clue.  What am I doing wrong?


texts: QGIS HTML in print composer, how to use "if" statement when building an expression?

texts:
I'm trying to include some images in a map layout to make an atlas.
Currently I'm using the following HTML
<img src=[%"FILEPATH_TO_IMAGE"%] height="400"/>
However, there are some situations where the atlas feature does not have a file path associated with it.
In this case it's shown as a question-mark instead.
Is there a way to add an IF-statement in the form of: IF filepath does not exist, nothing happens (do not execute the HTML img-block).


texts: Did leaflet-bing-layer change? (Angular 6)

texts:
So for the development of my app I've been using Esri's basemap.  I wanted to switch to bing.  This is my import as i saw in another question asked here.

import 'leaflet-bing-layer';

Then this is in my ngOnInit? The map loads fine but none of the other types besides the default Aerial will show up?  Any ideas

L.tileLayer.bing('key', {
      type: 'Road'
    }
).addTo(this.map);


texts: LiDAR data denoising algorithm complexity

texts:
I want to remove outliers (noise) from LiDAR data.
Currently I'm using the same algorithm as LAStools:

  The LAStools software version 171215 implements a straightforward outlier
    removal technique (lasnoise) based on the spatial frequency (SF) of the
    points.As described in the documentation, for each point in the dataset,
    a 3D structure - voxel - formed by 3x3x3 elements is generated having the
    current point as centre. If the number of points within this voxel
    (27 cells) is inferior to a certain threshold, the current point is
    considered an outlier (isolated point).

The problem is even if I subdivided the original dataset into smaller files, I still end up with around 4M data points. And because the algorithm complexity is O(n²) (for each point, check if any of the points is inside the voxel), the running time explodes around 500000 data points.

How do they manage the performance aspect and is there a better way to remove outliers ?


texts: How to find closest point on a line from a point in Oracle Spatial?

texts:
I want to find the nearest point on a line from a point in Oracle Spatial. I know how to do in SQL Server and PostgreSQL but could not find anything about Oracle.
I have used MDSYS.sdo_nn to find the closest line from a point, but I need the closest point on line as well to calculate distance from the start of that line.


texts:
I am attempting to call do a GetFeatureInfo request on it a counties layer which is stored withing ArcSDE and hosted via ArcGIS Server 10. The request looks like

WMS GetFeatureInfo Request

The county web service was created by way of ArcMap with the save service definition file. However, no data comes back with attempting to do run this command. Are there any settings on ArcGIS Server specific to whether a service can have GetFeatureInfo run on it? I was expecting values from the attribute table to come back.

Thanks,
Derek


texts: Clipping a raster from quickmapservices or XYZ tiles

texts:
I'm trying to clip a google satellite XYZ tile over a layer mask of my DEM layer to make a 3D map.
the problem is that after making the connection and adding it as a layer to my project, im unable to save it as a geotiff. It keeps returning an error saying "data source "filename.tif" is an invalid data source.

Moreover, I've noticed another issue, if i uncheck the vrt option while trying to export> save as geotiff, the file does not get saved at all!


texts: PostGIS returns an empty GeometryCollection out of an intersection between a MultiPolygon and a Point

texts:
Context

I've narrowed down an issue I faced involving many buffers and so on, to this in PostGIS (version 2.4):

SELECT St_geometrytype( St_geomfromtext(
    'MULTIPOLYGON(((0 0,4 0,4 4,0 4,0 0),
                   (1 1,2 1,2 2,1 2,1 1)),
                  ((-1 -1,-1 -2,-2 -2,-2 -1,-1 -1)))'
))

which answers: ST_MultiPolygon
It's OK, I guess.

But this;

SELECT St_geometrytype(St_intersection(St_geomfromtext(
    'MULTIPOLYGON(((0 0,4 0,4 4,0 4,0 0),
                   (1 1,2 1,2 2,1 2,1 1)),
                  ((-1 -1,-1 -2,-2 -2,-2 -1,-1 -1)))'
), 'POINTZ (30 40 50)') )

answers: ST_GeometryCollection
A GeometryCollection, ok, but out of what?

I even tried to force casting it into a Multi*Geometry but without any success, it still outputs a ST_GeometryCollection:

SELECT St_geometrytype(St_Multi(St_intersection(St_geomfromtext(
    'MULTIPOLYGON(((0 0,4 0,4 4,0 4,0 0),
                   (1 1,2 1,2 2,1 2,1 1)),
                  ((-1 -1,-1 -2,-2 -2,-2 -1,-1 -1)))'
), 'POINTZ (30 40 50)')) )

Question

I cannot figure out why this is giving a GeometryCollection where, obviously, it should return NULL or a "NULL first input geometry" (i.e. a NULL MultiPolygon)?


texts:
I am working on two different types of data. One is related to mapping of Genome structure and another is creating a layout of the inventories inside Internet Data Center.

What I am confused about is what kind of projection do we use (if any) is such conditions. The reason using GIS APIs here is very helpful is that these projects also require the same layering of the information. So how do I create or structure the data?

Should it also be shapefiles? What should their projections be? Can I use a shapefile to be hosted in Geoserver without projections?

Please let me know if you need any more info.


texts:
I want to import a ply file in GRASS using v.in.ply . However, I don't know what I should write in the "input =" field . Should I paste the way to my file on my computer (C:\Users\Marion\Documents\Photogrammetrie\test)? I tried doing so and it send me the error : "cannot open input file".
If not what should I put?


texts: How do I get the average point value in a raster with QGIS?

texts:
I have a point layer with standard deviation values and i want to get the average value of points in a 1 square-kilometer raster. I used the rasterize tool from QGIS to make a raster but this tool uses only one value of all points that lay in the raster.

Is it possible to do this?


texts: is there any tool to open and process MCD12Q1(modis landcover product)? apparently ENVI doen't work for this. any help highly appriciated.


texts: openFeatureForm in PyQGIS Standalone Application

texts:
I am trying to implement a desktop application with PyQGIS 3. When I try to use openFeatureForm to display the data of the selected layer object it shows an error, however when I use this code in the Python Console in QGIS it works correctly.
How can I use it for the application I'm doing outside of QGIS?
def info(self):
    self.vlayer = QgsVectorLayer('provincias.shp', "Provincias", "ogr")
    QgsProject.instance().addMapLayer(self.vlayer)
    layer = self.vlayer
    features = layer.getFeatures()
    feat = layer.getFeature(2)
    iface.openFeatureForm(layer, feat, True)

This is the console error :
Traceback (most recent call last):
  File "D:/proyectos/catastroesp/catastro.py", line 111, in info
    iface.openFeatureForm(layer, feat, True)
AttributeError: 'NoneType' object has no attribute 'openFeatureForm'


texts: Forcing Python Toolbox tool to break loop and do cleanup when user clicks Cancel?

texts:
I am writing a Python tool that perform heavy tasks in a loop. When the user clicks "Close" I want to break the loop, perform some cleanup, and then finish.

Consider this example tool that counts down from ten:

import arcpy
import time

class Tool(object):

   def __init__(self):
      self.label = "Cancel Test"
      self.description = "Tests how the cancel function works."
      self.canRunInBackground = False

   def execute(self, parameters, messages):
      for i in range(10, -1, -1):
         arcpy.AddMessage(i)
         time.sleep(1)
      arcpy.AddMessage("We have lift-off!")
      return

When I click "Close" in the beginning of the execution it does not brake the loop or give me an opportunity to perform any clean up. Instead it loops to the very end and adds the message "We have lift-off!" before it outputs this:

Completed script CancelTest...
Cancelled function
(CancelTest) aborted by User.
Failed at [TIME] (Elapsed Time: 10,01 seconds)

I have found documentation for ArcGIS Pro (that uses Python 3.4) that describes how you can do this there using isCancelled:

import arcpy
import time

#Make sure it does not auto cancel.
arcpy.env.autoCancelling = False

class Tool(object):

   def __init__(self):
      self.label = "Cancel Test"
      self.description = "Tests how the cancel function works."
      self.canRunInBackground = False

   def execute(self, parameters, messages):
      for i in range(10, -1, -1):
         arcpy.AddMessage(i)
         time.sleep(1)
         #Check if the user clicked "Close"
         if arcpy.env.isCancelled:
            arcpy.AddMessage("Launch aborted!")
            return
      arcpy.AddMessage("We have lift-off!")
      return

However, when I run this from ArcCatalog 10.3 (i.e. not using ArcGIS Pro) it gives me the following error:

Traceback (most recent call last):
  File "H:\Mina Dokument\DGD\Python\CancelTest.py", line 19, in execute
    if arcpy.env.isCancelled:
AttributeError: 'GPEnvironment' object has no attribute 'isCancelled'

Is there anyway to mimic the behavior that is available in ArcGIS Pro in an ordinary Python toolbox using Python 2.7?


texts:
When I execute a grassgis alogrithm in Qgis I always get following error:

  'ascii' codec can't encode character u'\xfc' in position 23: ordinal
  not in range(128) See log for more details

and

  This algorithm cannot be run :-(  It seems that GRASS GIS 7 is not
  correctly installed and configured in your system. Please install it
  before running GRASS GIS 7 algorithms.

In processing opitions I set the path for grassgis to
C:\OSGeo4W64\apps\grass\grass-7.2.0, same for grassgis docs.

There is no option for setting the path to msys.

I reinstalled whole OSGeo4W with coustom options and installed all available options. No result.


texts:
I have a large area with lots of small 'sites' within it.  I want to label each of these 'sites' but when I zoom out to 1:95,000 the labels overlap so I cannot read some of them.  Is there a way to prevent the labels from overlapping with each other but still show the area?

I have added an example below.  I really want to be able to read all the labels, no matter if the label is directly above the 'site'.


texts:
I want to prevent the user from starting an edit session under certain conditions.

How can I do that?

void m_editEvents_OnStartEditing()

        {
            // Code
        }


texts:
I get geoserver 2.8-SNAPSHOT code from github. I delete some madules of it.

Geoserver developed with Spring framework, Now I want to add my Controllers to this project. Some of My controllers have view pages.

How do I do?

Where do I add my Controller? No body can Help me please?


texts: Can I use esri-leaflet together with arcgis api?

texts:
I am not sure I have understood the co-existence of arcgis javascript api with the  esri-leaflet library.

Can I use leaflet together with the javascript api.

For example can I have my map built with the javascript api (new esri.Map(...) )
and then add an L.esri.Dynamic layer ?

If not, I suppose I will have to choose which of the two different paths I want to follow.


texts: How to choose workspace when uploading SLD to Geoserver with the REST API

texts:
I have 200+ SLD's I would like to upload to Geoserver using the REST API. I successfully uploaded an SLD using the gsconfig Python module.

However, I would like to be able to set the workspace when I upload. Is this possible? Below is my script so far.

from geoserver.catalog import Catalog

cat = Catalog("http://chchchch:8080/geoserver/rest", "username", "pw")

with open("C:\Users\Documents\SLD\my_sld.sld") as f:
    cat.create_style("my_sld", f.read())


texts: Changing time units of direction window in ArcGIS network analyst?

texts:
I want to change direction window unit from minutes to seconds . I set travel time unit to

second and the result of the direction window show second but the unit remains min .For

example in the below image, the numbers are in second but unit remains min.How can i change unit from min to sec ?


texts: Linking demographics data to zip code, census data

texts: Where do I find population data to link to DMA / MMA or zip code? I need to be more granular then state. Any free data source? ArcGIS?


texts: Creating shapefiles with attached/hyperlinked PDF forms in ArcMap?

texts:
What is the best way to attach pdf forms to point shapefiles to fill out while out in the field?

I know how to create a hyperlink to the specific folder but is there a way to have the shapefile created with the link or attachment already attached?


texts:
I don't have the acceleration tab in QGIS:

QGIS: 3.10.1-A Coruña

Nvidia driver: 440.44

Graphics card: RTX 2070 super

Cuda: 10.2.89-3

Operative system: manjaro - Linux

Edit, i also don't see any option related to GPU acceleration in the rendering tab:

Edit 2, I am also don't see any OpenCL in Advanced > core, using the search bar doesn't show any OpenCL options either.


texts: How to make the Overlays section in an OpenLayers2 LayerSwitcher exclusive?

texts:
At the simplest level I need to figure out how to create a customized OpenLayers LayerSwitcher Overlays section which is not using Checkboxes, but uses RadioButtons instead.

For my application its not desirable for the user to be able to select multiple layers.  Eventually I'm going to need to be able to create group sections in the Overlay section of the LayerSwitcher, so if anyone has any hints on that I'd appreciate it.

I've seen the customized LayerSwitcher example on the OpenLayers site, but I can only get it to work in Chrome.  I need to be able to support Safari on an iPad.


texts: How to symbolize points using a numbering system

texts:
Using the ArcGIS Server JavaScript API version 4, is it possible to symbolize points using a numbering system (1-n) as in this mockup?

A simple marker symbol or picture marker symbol can be used to generate the circles, but how can the numbered labels be generated within the circles?


texts: Find nearest neighbour for set of points

texts:
Using PostgreSQL and PostGIS. I have two sets of point data (geoms), one is a node graph indicating roads and the other represents street crossings. For my set of crossing data I want to find the single closest node to all the crossings. And use these geoms for further analysis. However, I am having trouble with nearest neighbour techniques as they mainly seem to only function for one point instead of for a set of points e.g.
SELECT
  streets.gid,
  streets.name
FROM
  nyc_streets streets
ORDER BY
  streets.geom <->
  (SELECT geom FROM nyc_subway_stations WHERE name = 'Broad St')
LIMIT 10;


texts: What do the quantitative values for the colour ramp gradient for heatmap style in the layer properties for QGIS signify?

texts: I would like to understand whether any quantitative values can be extracted when using the heatmap style in the layer properties? I assume areas with colours related to 1 are areas which have the highest risk and those with 0 are at least risk, however what can be said about the values in between? For example, in the image here, what do colours above 0.5 in the plot signify?


texts: How to calculate marker area on OpenLayers 3?

texts:
I have the following scenario: a polygon (or few) and a marker which is displaying some data of its related geometry/ies:

What I want is to manage the map's zoom level in order to decrease the size of the marker when it is too large for the background polygon which is hovering. I thought to calculate the area of the marker every time the zoom changes. I had this algorithm working fine for Google Maps but I can't get the solution.

This is the algorithm for Google Maps:

// Get the position of the marker
var pos = this.marker.getPosition();

// Calculate the height and the width to compare with
var height = 32 / 2;
var width = 60 / 2;

// Get the scale and the bounds taken up by the map indicator
var scale = Math.pow(2, this.map.getZoom());
var proj = this.map.getProjection();
var wc = proj.fromLatLngToPoint(pos);
var bounds = new google.maps.LatLngBounds();
var sw = new google.maps.Point(((wc.x * scale) - width) / scale, ((wc.y * scale) - height) / scale);
bounds.extend(proj.fromPointToLatLng(sw));
var ne = new google.maps.Point(((wc.x * scale) + width) / scale, ((wc.y * scale) + height) / scale);
bounds.extend(proj.fromPointToLatLng(ne));

// Get the locations delimited by the map indicator bounds
var sw = bounds.getSouthWest();
var ne = bounds.getNorthEast();
var southWest = new google.maps.LatLng(sw.lat(), sw.lng());
var northEast = new google.maps.LatLng(ne.lat(), ne.lng());
var southEast = new google.maps.LatLng(sw.lat(), ne.lng());
var northWest = new google.maps.LatLng(ne.lat(), sw.lng());

// Calcualte the area
var area = google.maps.geometry.spherical.computeArea([
    northEast, northWest, southWest, southEast
]);

And this is my adaptation to OpenLayers 3:

var pos = [item.lng, item.lat];
var scale = Math.pow(2, zoomLevel);
//var scale = 1;

// Calculate the height and the width to compare with
var height = 32 / 2;
var width = 60 / 2;

var wc = map.getPixelFromCoordinate(ol.proj.transform(pos, 'EPSG:4326', 'EPSG:3857'));

var sw = [((wc[0] * scale) - width) / scale, ((wc[1] * scale) - height) / scale];
var nw = [((wc[0] * scale) - width) / scale, ((wc[1] * scale) + height) / scale];
var se = [((wc[0] * scale) + width) / scale, ((wc[1] * scale) - height) / scale];
var ne = [((wc[0] * scale) + width) / scale, ((wc[1] * scale) + height) / scale];

var posSW = ol.proj.transform(map.getCoordinateFromPixel(sw), 'EPSG:3857', 'EPSG:4326');
var posNW = ol.proj.transform(map.getCoordinateFromPixel(nw), 'EPSG:3857', 'EPSG:4326');
var posSE = ol.proj.transform(map.getCoordinateFromPixel(se), 'EPSG:3857', 'EPSG:4326');
var posNE = ol.proj.transform(map.getCoordinateFromPixel(ne), 'EPSG:3857', 'EPSG:4326');

var extent = ol.extent.boundingExtent([posSW, posNW, posSE, posNE]);

var wgs84Sphere = new ol.Sphere(6378137);
var area = Math.abs(wgs84Sphere.geodesicArea([posSW, posNW, posNE, posSE]));

Is there any other algorithm to calculate the surface of a marker on OpenLayers 3 or a better solution for my final goal?


texts: Conflicting WGS84 coordinates between Trimble and Garmin GPS receivers

texts:
I am a contractor at a government facility in Nevada. We have our own in-house survey department who, for a number of reasons, still uses NAD27 State Plane Coordinates (Nevada Central Zone 2702) in survey feet (referred herein as SPC27).  They also provide lat-long coordinates in NAD83 degrees, converted using Corpscon, I believe. For my project, they have established several local control points.  When I occupy these same control points with my 15-year old but WAAS-enabled hand-held Garmin 76, I get WGS84 coordinates (and SPC27 coordinates converted using Corpscon) very close to those of the surveyors (within several feet).

When my co-worker occupies the same control points with his Trimble GeoHX (differential GPS with post processing), he gets WGS84 coordinates that are consistently 0.0009 degrees farther east, or when converted to SPC27 with his older CoordTrans software, are 250 feet farther east (and 18 feet farther north).

But when he converts his seemingly wrong WGS84 coordinates to SPC27 using his newer Trimble Pathfinder Office software, his resulting SPC27 coordinates plot closely with mine and those of our surveyors.

So, occupying the same control point, my Garmin and his Trimble give WGS84 readings that are off by 0.0009 degrees longitude (or 250 feet when converted to SPC27 using Corpscon). Yet this difference seems to be corrected when his WGS84 coordinates are converted to SPC27 using his Pathfinder Office software.

The error (250 feet) is about what I've read to expect in my geographic region between SPC NAD27 and SPC NAD83.  But why would our raw WGS84 GPS readings be this far apart?


texts: Obtaining print quality PDF or PNG output from QGIS?

texts:
I work with people who use ArcGIS Desktop and they have found the best resolution output images are made by creating PDFs. They tell me that *.png and the rest out of ArcGIS Desktop are not satisfactory and that they do not bother with them. Typically we are making the output appear in a report and so the PDFs of the maps and words are assembled in a PDF editor package to make the final document.

As a result I have been toying around with QGIS on my Windows 7 machine and notice that the best output comes by printing directly out of the map composer.

This output on my Docuprint C2120 printer is as good as any of the ArcGIS Desktop PDFs and is truly impressive. QGIS makes very attractive maps in this configuration.

However, I do not seem to be able to create PDFs or .pngs and then print them which are anywhere near as good. I find the .svg output to Inkscape unrealiable where, for instance, the grid/graticule may come out oversized.

By way of example, the good quality output is where a sub-6 pt font can be read easily off the printed page.

Any thoughts?

I am using QGIS trunk.


texts:
How to show a gif image for a point in GeoServer SLD like blinking point?

I found here animated wms layer. Please help me.

I tried SVG Icon. Its working without animation. My svg code is

<svg width="500" height="500" viewBox="0 0 320 320">
        <path
            fill="#FF0000" stroke="#FF0000"
            d="M40,50a10,10 0 1,0 20,0a10,10 0 1,0 -20,0">
            <animate
                attributeType="XML"
                attributeName="fill"
                values="#900;#f00;#800;#800"
                dur="0.8s"
                repeatCount="indefinite"/>
            </path>
 </svg>


texts: ArcGIS Server connection error for feature layers

texts:
I get the following error when trying to access my hosted feature services from ArcGIS Server: "Connection refused. Check that the hostname and port are correct and that the postmaster is accepting TCP/IP connections." Other service types (Image layers, Scene layers) are working fine. This error appeared suddenly and everything was working fine before.

What might be causing it?


texts: Grouping map points into fixed cluster sizes?

texts:
I have a dataset of 655 lat/long pairs which I'd like to split into about 100 groups.  A group should have 5-10 pairs which are geographically close to each other.  Dense groups should have more points, sparse groups should have less. For example urban groupings should be larger, rural ones smaller.

Is there an established algorithm for doing this kind of grouping, or am I going to have to design one from scratch?

I'm using google maps v3 api to display this data, but as it's a fixed dataset I'm prepared to do some offline number crunching.


texts: Appending map to another map using GRASS?

texts:
I'm trying to simply append a map to other map, without any sort of geometric operations as merge or intersects, just the original geometries with theirs attributes.
This is what I already tried:

# input maps are poly2 and poly3
# First, i change the categories of poly3 to avoid categories conflicts:
# see http://grass.osgeo.org/grass64/manuals/v.patch.html
v.category in=poly3 option=sum cat=10000 out=poly3_category --overwrite
g.remove vect=poly3
g.rename vect=poly3_category,poly3
# I update the cat column of the table as well (why this is not done
# automatically by v.category??)
v.db.update poly3 col=cat qcol="cat+10000"

# now the to maps are patched together
db.droptable -f table=polygons
v.patch -e input=poly2,poly3 output=polygons --overwrite

# export to result as geojson:
v.out.ogr in=polygons olayer=polygons format=GeoJSON dsn=/tmp/polygons.json type=area

The geometries of the result map polygons.json look good, but the attributes are messy: I can not distinguish which records come from poly2 and which ones from poly3. Many records have only empty attributes. What I'm doing wrong?


texts:
I'm writing a short script that iterates over a list of feature classes and sets values for a specific field ("Contour_Type"), based on the values in two other fields ("Type" and "Depression"). It seems my script does not pass the 'with arcpy.UpdateCursor(i) as cursor' line. I've searched for an answer, but can't find anything that pertains to my issue.

EDIT: The following code now works!

inputloc = arcpy.GetParameterAsText(0)
arcpy.env.workspace = inputloc
inputfcs = arcpy.ListFeatureClasses()

    for i in inputfcs:
        with arcpy.da.UpdateCursor(i,["Type","Depression", "Contour_Type"]) as cursor:
            for typ, dep, con_typ in cursor:
                if typ == "Intermediate" and dep == "Yes":
                    arcpy.AddWarning("Adding Values..")
                    row = (typ, dep, "Intermediate_Depression")
                    cursor.updateRow(row)
                elif typ == "Index" and dep == "Yes":
                    row = (typ, dep, "Index_Depression")
                    cursor.updateRow(row)
                else:
                    row = (typ, dep, typ)
                    cursor.updateRow(row)


texts: Error while running PySAL SpaceTimeEvents function

texts:
I'm new to PySAL and trying to run Knox test to find space-time correlation in my dataset. My original dataset is in the form of csv including attributes such latitude, longitude, date and other related features. I converted the csv file into shapefile (point-type shapefile) using QGIS. I stored the .shp, .dbf and .shx in a folder in pysal/examples directory.

I'm trying to run these lines of code:

import numpy as np
import pysal
import pysal.spatial_dynamics.interaction as interaction
np.random.seed(100)
path = "...pysal/examples/Shapefile/trial.shp"
events = interaction.SpaceTimeEvents(path, 'Date', infer_timestamp = True)

(source: http://pysal.readthedocs.io/en/latest/users/tutorials/dynamics.html#space-time-interaction-tests)

I'm facing the following error even though the trial.dbf file is present in the same folder:

  KeyError: 'trial.dbf not found in PySAL built-in examples.'

Can anyone guide me as to what I'm doing wrong here?

What is the solution for this error?


texts: How to get all ID's for a line in a network of lines in order

texts:
I have a path (pink) in a network (green) and I want the id's of all the lines that make up that path. They are two separate layers.

What I need returned is a list of id's in their order starting from one end of the path going to the other.
A short path might be: Starting Id 118395-->133285-->283495-->283490-->293505 End Id
I have been able to succesfully select the id's using a Buffer around the path then using the Select by Location tool to find the network lines inside that buffer. The issue here is I cant work out what order the selected lines are in.
The path is calculated using the network analysis shortest path (point to point) tool if anyone knows how I can get this to return the id's in order this will solve the problem.
The network (green) doesnt have any direction in its attributes.
The path (pink) only has a start and end point and a 'cost'.
Does anyone know how to do this?
Thanks


texts: Scripting a population sum based on union

texts:
I'm having trouble automating a process to sum one field into a new field, summing whenever the value of an FID field from my Union is the same.

I'm looking to sum the population, stored in field "TheissenPopNonSum", based on polygons which share the same value in "FID_shape", into a field called "TheissenPopSum".

Below is an example of all the polygons that share "FID_shape"=543.

Each of the selected has a population associated that I'm trying to sum into the Theissen polygons generated from the point data.

I know this can be programmed but I don't know how to nest it, so it goes through each FID_shape, sums the associated Population, then moves on.

Here is what I was thinking, but fails:


texts: Downloading data sets during installation

texts: I can successfully install QGIS 2.8.1 on a Windows 7 machine but I get an error when downloading the Carlina, South Dakota and Alska data sets. The error message is http/1.0 504 Gateway Timeout. Looks to me like as I am behind a firewall, I need to supply login credentials. But how would I do that?


texts: Select all columns except geometry using virtual layers

texts:
Is it possible to select all columns except geometry using virtual layers in QGIS?

For example I have this (simplified) query:

Select *, st_buffer(geometry, 100) From mypolygons

It does not work. It does not buffer at all. But using the (simplified) query as follows does buffer all polygons as intended:

Select id, name, st_buffer(geometry, 100) From mypolygons

Now the issue is, that if I have polygons with a lot of columns I need them to type in all manually. Is there a way to avoid this?


texts: Hillshading makes crater-like depressions look elevated?

texts:
I have a 1 m resolution DEM of the seafloor of a field of depressions. Hillshading the DEM makes the depressions look like mounds instead.

What I can do to fix this?

My raster is projected in UTM, horizontal datum is GCS 1984, my data frame is projected, and I'm using ArcGIS Desktop 10.5.


texts: Can X-Ray for ArcCatalog be used without Excel? Our organization went with Libre Office as a budget savings move. Calc in that suite can read/write .xls files, but X-Ray appears to require Excel itself (presumably it's calling Excel.exe). Thanks.


texts: Convert sgrd (SAGA GIS format) to raster or ascii or with R

texts:
I need to convert a sgrd file to tif or ascii using R. I converted a big sgrd file (1.4 GB, 376521831 pixels) to ascii using the package RSAGAbut it takes too long (approx. 45 minutes):

rsaga.sgrd.to.esri(in.sgrds = "path\\dem_filename",
                 out.grids = "output_filename",
                 out.path = "path\\foldername",
                 format = "ascii",
                 georef = "corner",
                 prec = 5,
                 env = work_env)

How could I convert sgrd format with rgdal or other R packages?


texts:
from osgeo import gdal
ds = gdal.Open('/Data/F182010.v4/F182010.v4d_web.stable_lights.avg_vis.tif', gdal.GA_ReadOnly)

The code showing Runtime Error Permission denied How to fix this?

/Data/F182010.v4/F182010.v4d_web.stable_lights.avg_vis.tif: **Permission denied**


texts:
I'm creating a simple print layout with Python.
I've created a legend but I want to remove a few items from it.

Here is my code to create a blank Print Layout called "LayoutWithLegend" and add a legend to it:

from qgis.core import *

projectInstance = QgsProject.instance()     #creates a project instance
manager = projectInstance.layoutManager()

layout = QgsPrintLayout(projectInstance)    #needs a QgsProject instance as an argument
layout.initializeDefaults()
layout.setName('LayoutWithLegend')

manager.addLayout(layout)                   #add layout to manager. layout is a project instance

itemLegend = QgsLayoutItemLegend.create(layout)

itemLegend.attemptResize(QgsLayoutSize(60, 40,
QgsUnitTypes.LayoutMillimeters))

layout.addLayoutItem(itemLegend)

When created, my legend looks like this :

I want to remove the 'Patient Data' and 'Google Maps' basemap from the legend. I can't figure it out but some of my best efforts to try are:

-Looking through QgsLayoutItemLegend methods but none of those appear to be relevant

-Playing with QgsLayerTreeGroup. This class has a removeLayer() method.

Any help?


texts: I want to download image from QGIS. What should i write in Resolution setting(2 dialog boxes horizontal and vertical)? What do they mean?


texts: How to show a map with queries from a map file connected to postgres in browser using javascript?

texts:
So, the problem is how to show a map in browser using JavaScript.

I have a map file with 2 layers in it: "peta" is a shp data from PostGIS and "korban" is a numeric data from PostgreSQL, and then I have a HTML file, in that file I use JavaScript to show a toolbar but I don't know how to show the map from map files.

I have tried to call map file with "http://my.site/cgi-bin/mapserv.exe?MAP=C:/ms4w/Apache/htdocs/webgis.map&LAYERS=all&MODE=MAP" and the map shows but it cannot be zoomed, so I need to know how to show map with a query?

The query I mean is something like this

var query = 'select ' + 'admin_90_des(edit_pande,jawa)_2010."gid", ' + 'data_korban."2009", ' + 'admin_90_des(edit_pande,jawa)_2010."kel_desa", ' +
'admin_90_des(edit_pande,jawa)_2010."kecamatan", ' + 'admin_90_des(edit_pande,jawa)_2010."geom"' +
'FROM "admin_90_des(edit_pande,jawa)_2010" ' + '"data_korban"' +
'where "admin_90_des(edit_pande,jawa)_2010".gid="data_korban".gid ';

new_dbd = new OpenLayers.Layer.WMS("Peta", "http://my.site/cgi-bin/mapserv.exe?MAP=C:/ms4w/Apache/htdocs/webgis/webgis02.map&LAYERS=ALL&MODE=MAP&korbanDBD=geom FROM (" +
query +
") AS foo USING UNIQUE gid and USING SRID=-1", {
    layers: ['korban', 'peta'],
    format: 'image/png',
    transparent: true
}, {
    isBaseLayer: false
});

map.addLayers([new_dbd]);

But I don't understand how it works so please anyone can explain to me or at least give me some clue how can I show my map in browser?


texts: Changing Coordinate system from custom CRS?

texts:
Is it possible to reproject a shapefile from a custom CRS?

My example is as below but the EPSG is not a valid EPSG supported file.

    ogr2ogr -s_srs EPSG:100000 -t_srs EPSG:27700 output.shp input.shp


texts: Programmatically open Viewer/Magnifier window in ArcMap 10.2+

texts:
Is there a way to programmatically open Viewer/Magnifier etc. windows of ArcMap from Add-In?
I am using 10.2, but anything 10.2+ will be fine.


texts: Converting coordinates in metres to Longitude & Latitude (CRS conversion) in OGC WKT format?


texts: Distinguish between the types of cultivated crops in remote sensing

texts: I recently produced map of LULC but there is the difficulty of distinguishing between the types of  crops. Is there a method to distinguish between the types of cultivated crops such as wheat, barley, maize etc.?


texts: How can I activate interactivity on cartoDB.js with createLayer?

texts:
I've made a map with cartoDB.js with the createLayer method. I've defined the sublayers,inside a variables called layerSource, giving them an interactivity field like this:

sublayers:[ {sql: "SELECT * FROM db_name",
interactivity: "cartodb_id",
cartocss: "styles"},...more sublayers...

After this I've added the layers to the map

cartodb.createLayer(map,layerSource)
                    .addTo(map)
                    .done(function(layer) {
                       for (var i = 0; i < layer.getSubLayerCount(); i++) {
                           sublayers[i] = layer.getSubLayer(i);
                       }
                    })
                    .error(function(err) {
                        console.log("error: " + err);
                    });

But the sublayers are not interactive. I've seen in the documentation that I have to enable the interactivity with setInteraction:true, but where I should activate this option?? Inside the createLayer (as an argument?), inside the .done??.. I assume this should be something very obvious but I'm new to javaScript and not being able to get this done.

EDIT:

After @Andy_Eschbacher answer.
I'm trying to do what you said but with the first block of code

.done(function(layer) {
                      for (var i = 0; i < layer.getSubLayerCount(); i++) {
                           sublayers[i] = layer.getSubLayer(i);
                       };

                        sublayer[0].on('featureClick',function(event, latlng, pos, data, layerIndex) {
                            console.log(data);
                        });
                    })
                    .error(function(err) {
                        console.log("error: " + err);
                    });

the console sends an error because sublayer is not defined, so I define a sublayer = layer.getSubLayer(0);as you do in this tutorial but no action is triggered.


texts: Proper way to deal with repeated observations in IDW interpolation

texts:
I am trying to apply IDW interpolation in a situation where I have multiple clusters of points, and all points in each cluster have identical coordinates.
My understanding is that different IDW implementations have different ways to deal with repeated measurements (see for instance this and this posts)
Currently I am using the R package gstat v2.0-6, and I can't really figure out how it addresses this issue.
So my question is: what is a proper/accepted way to deal with points with identical coordinates when using IDW? I guess that my options are:

average all values from repeated observations before running the IDW algorithm
add some minute spatial noise so that the distance between any two points is always > 0
?

Switching to other interpolation methods is probably not at option, given the very limited extent of the data set.


texts: Difference between gdalwarp and projectRaster

texts:
I am trying to project a Raster. In R there is the projectRaster() function to to this (below a fully reproducibly example) :

# example Raster
require(raster)
r <- raster(xmn=-110, xmx=-90, ymn=40, ymx=60, ncols=40, nrows=40)
r <- setValues(r, 1:ncell(r))
projection(r)
# project to
newproj <- "+init=epsg:4714"

# using raster package to reproject
pr1 <- projectRaster(r, crs = CRS(newproj), method = 'bilinear')

Which works fine. However it is quite slow.

In order to increase speed I though to use gdalwarp instead (with a SSD the cost of reading and writing from/to disk/R are not very high).

However, I cannot reproduce the results of projectRaster() using gdalwarp:

# using gdalwarp to reproject
tf <- tempfile(fileext = '.tif')
tf2 <- tempfile(fileext = '.tif')
writeRaster(r, tf)
system(command = paste(paste0("gdalwarp -t_srs \'", newproj, "\' -r bilinear -overwrite"),
                       tf,
                       tf2))
pr2 <- raster(tf2)

It seems to work, however the results are different:

# Info
system(command = paste("gdalinfo",
                       tf))
system(command = paste("gdalinfo",
                       tf2))

# plots
plot(r)
plot(pr1)
plot(pr2)

#extents
extent(r)
extent(pr1)
extent(pr2)

# PROJ4
proj4string(r)
proj4string(pr1)
proj4string(pr2)

# extract value
take <- SpatialPoints(matrix(c(-100, 50), byrow = T, ncol = 2), proj4string = CRS(newproj))
plot(take, add = TRUE)
extract(pr1, take)
extract(pr2, take)

What am I missing / doing wrong?

Are there other (faster) alternatives to projectRaster()?


texts: Python script kills the data of QGIS Project

texts:
i´ve a problem with a script, working fine until last month and there is no idea coming up what to do.
I use PyCharm 2017.3.7 for this script because every Version above 2018.3.3 do not work with QGIS

It is a complete easy basic script, loading and saving a python project as a batch function to handle bad layers.
But now the script kills all layers of my project so i get an empty one back.
Any ideas for solving the problem. Maybe another code....

from qgis.core import *
import os
QgsApplication.setPrefixPath(r"C:\OSGeo4W64\apps\qgis-ltr", True)
qgs = QgsApplication([], False)
qgs.initQgis()

ordner = r"C:\Users...."

for root,dirs, files in os.walk(ordner):
    for file in files:
        if "QGZ" in file.upper():
            project = QgsProject.instance()
            project.read(os.path.join(root,file))
            QgsProject.instance().write()

qgs.exitQgis()


texts: Classification result from SAGA GIS shifted by 5 meters if used in ArcMap?

texts:
I use a band composite of 10 Sentinel-2 bands for my classifications. If I use the same composite in SAGA to make a classification this classification result is shifted by 5 meters to the left and to the bottom.

Since I used the same band composite I do not understand why is the result from SAGA shifted when should be exactly above the other result from ArcGIS.


texts: Extract specific geometry type from a QuickOSM query in QGIS graphic modeller

texts:
I am trying automate the extraction of various features from OSM using QuickOSM in GQIS graphic modeller. I have the query etc running fine but I want to end up with just the lines and export to a shapefile.
How best to select the lines and save as a shapefile?

My query is set to "type="way"" but I still get other geometry types
<osm-script output="xml" timeout="25">
<union>
    <query type="way">
        <has-kv k="highway" v="cycleway"/>
        <bbox-query {{bbox}}/>
    </query>
</union>
<union>
    <item/>
    <recurse type="down"/>
</union>
<print mode="body"/>


texts:
I'm adding a custom form to one of my attributes fields, I want to create a drop-down menu where there is predecided values to choose from.

I'm going into the layer property>Forms>Selecting the target field, then I chose Widget type: Value map.

This works and I get a drop-down menu as I wanted. But this doesn't stay permanent, whenever I close and reload the field it goes back to Widget-type "text editing".

Is there a way I could keep my form permanent?

Using QGIS 3.4.5


texts: Seeking algorithm in QGIS Process Modeller to create vector grid?

texts:
I am using Process Modeller to automate merging several vector layers and then creating a grid to cover the extent of the merged layer. There are functions which already exist which are not listed in the modeller such as :

Vector -> Research Tools -> Vector Grid

MMQGIS -> Create -> Create Grid Layer

The QGIS geoalgorithm Create grid in the modeller is not useful to me as I have to manually determine the position and set the CRS etc beforehand.

Is there any other plugin or function which allows me to accomplish this?

I am now trying to get the Processing plugin (formerly Sextante) to recognise the Vector Grid script. The file is named "doVectorGrid" which I renamed to simply "VectorGrid". I followed this advice: http://planet.qgis.org/planet/tag/sextante/ where it mentions at the bottom how to add your own modified algorithm.

I have added the new terms into the QGISAlgorithmProvider.py (the equilavent to FToolsAlgorithmProvider.py) but I get the following error: ImportError: cannot import name VectorGrid

I'm guessing I have to change the code of the Vector Grid script but I am new to Python. I have copied the ui_frmVectorGrid.py file into the same directory.


texts: Aligning points labels to another line in QGIS?

texts:
These days I've found very difficult to position my label points parallel to lines such as rivers and streets.

I have these two shapefiles, one is the point shape with names and the other is the polygon shape with the Streets.

This is the visualization:

Then when I'm trying to align the "Name" of the street parallel to the "Street" Shape with the label settings and there is no option to do it.

The result  I am expecting is something like this:

Can anyone help me on this?


texts:
I am trying to use pyproj and transform but I have an error name 'pyproj' is not defined.
I check my environment and I can see the pyproj version 2.6.0. My python version is 3.7.3. So, what's going wrong there? Why pyproj is not defined?

project = partial(pyproj.transform, pyproj.Proj(example1.crs), pyproj.Proj(example2.crs))


texts: How do I get rid of edge effects while using focal in R to smooth a raster?

texts:
I am trying to use R with the focal function from the raster package to smooth raster images.

I am working with Radar images that have already been preprocessed. As part of that, some calculations introduce unrealistic values that are then deleted from the raster.
So basically, I am working with images with values ranging from -2 (black) to 7 (white) like this and want to smooth over the small datagaps:

To do so, I tried the focal function from the raster package with     the following code:

raster1 <- raster("example_raster.tif")

mat <- matrix(1/25,ncol=5, nrow=5)
raster_res <- focal(raster1, mat, FUN="mean", na.rm=T)

writeRaster(raster_res, filename="smooth_raster.tif", format="GTiff")

The problem with this is, that the NAs seem to introduce some kind of edge effect. na.rm is set to TRUE to exclude the NAs from the actual calculations, but if you look at the edge at the east or how spotty the areas around former NA values appear, it looks like R still somehow includes them in the moving window calculations, somehow influencing the sourrounding values.
Here is what the output looks like:

To check, if this is a normal focal-problem, I also used r.neighbors in the QGIS processing toolbox set to a window size of 5 to smooth my data and it does exactly what I expected. Unfortunately I can't post a third link due to being new, but the raster is very smooth and does not have the spotty appearence of the output R produces.

I already tried using the pad and padValue arguments in focal but did not manage to make it work the way I want. Unfortunately, I don't know enough about GRASS to understand, what exactly r.neighbors does differently than focal.
Since this will be part of a script that automatically processes several rasters, doing everything manually with QGIS is not really an option.


texts: Can PostGIS functions be written in Java?

texts:
I am wondering if there is a possibility of programming new functions for PostGIS using Java, so that you can refer them same way as classic PostGIS functions are. Something like this:

CREATE OR REPLACE FUNCTION public.st_line_interpolate_point(geometry, double precision)
  RETURNS geometry AS
'$libdir/postgis-2.0', 'LWGEOM_line_interpolate_point'
  LANGUAGE c IMMUTABLE STRICT
  COST 1;

I know that most of PostGIS functions are written in C however I am not very familiar with that language.


texts: Transfering attributes from polyline to polygon features using ArcGIS Desktop?

texts: How can I transfer the fields and attributes from a polyline feature to a polygon feature in ArcGIS 9.3?


texts: Select empty Postgis (Postgresql) geometry column

texts:
I try to select rows in Postgis that have an empty geometry column.

This works but gives me the opposite:

SELECT * FROM foo where geom <> '';

And this fails:

SELECT * FROM foo where geom = '';

>>>ERROR: parse error - invalid geometry
>>>SQL state: XX000
>>>Hint: You must specify a valid OGC WKT geometry type such as POINT, LINESTRING or POLYGON

And this returns nothing:

SELECT * FROM foo where ST_IsEmpty(geom);

Does anyone know the correct syntax?


texts:
Is there a way to enclose my conic shaped map in graticules on ArcMap? I can only figure out do this which doesn't look very nice:

I want something like this, I think this was done in GMT though:


texts: Looking for formula for the Robinson projection

texts: I'm looking for a source that describes how the Robinson projection is calculated. I'm going to use this for converting long/lat into Robinson coordinates in both R and in javascript.


texts: How do I know the vertical/horizontal units of my DEM? I'm creating a slope using gdaldem

texts:
I'm using gdaldem to create a slope map.

I can apply a scale. The default is 1.

The notes say:

Ratio of vertical units to horizontal.
If horizontal unit of source DEM is degrees (e.g WGS84 projection), use scale=111120
if the vertical units are meters (or scale=370400 if they are in feet)

My DEM is in NAD27 UTM 18N and in the projection properties it says linear units Meters and angular units degrees.

How do I know what the vertical units are?


texts: Seeking shapefile for postal codes/zip codes of India?

texts:
I have been searching for free/paid (preferably free) sources for postal code shapefiles for India.

Does anyone know where I can find one?


texts: How to specify units when using SelectLayerbyLocation_management WITHIN_A_DISTANCE?

texts: Can anyone help by suggesting how to specify the units when using the SelectLayerbyLocation_management tool. I am selecting by "WITHIN_A_DISTANCE" and the ArcGIS resources do not suggest how to do this. I am getting some funny results at the moment.


texts: How to fix the broken MGRS Tools plugin?

texts:
Here the case: I have a plugin installed in .qgis2/python/plugins that remains broken and I can't make it work. I've found it here https://github.com/volaya/mgrs-tools. If any of you know the way to fix it I'll be so glad.

Hint: from forum I've heard it can come from dash characters in the plugin's folder name.

I use Windows 7 64-bits and QGIS 2.6

UPDATE:

It appears that running QGIS 2.10 32-bits solves my problem.
Meanwhile as the issue persists with 64-bit versions of QGIS I've just partly found an answer.


texts: Display Popup from postgis db over WMS layer file using leaflets?

texts:
i want to display popup information over WMS service file that i published from geoserver and using it Leaflets. Below is the picture for illustration.

Now i just want when user click any polygon the information related to it will be display as popup from postgis database, what are the possible ways to fetch information db ?

I am new to WebGIS.


texts: Using Python to field calculate Lat/Long with decimal degrees rather than meters

texts:
The below is a snippet of a python script that I am writing that basically adds two fields - Latitude and Longitude, then calculates them. I could use the Add XY Tool which gets me to the same point. The way that I did it with the expression3 and expression4 variables should work but doesn't. The projection is WGS 84 the script works fully from start to finish. The columns are calculated in meters rather than decimal degrees which is what is needed for the project.

Any thoughts?

...
...
...

## Add Lat/Long Fields
inFeature = "Point_Feature_Class"
fieldName1 = "Latitude"
fieldName2 = "Longitude"
fieldType = "DOUBLE"
arcpy.AddField_management(inFeature, fieldName1, fieldType)
arcpy.AddField_management(inFeature, fieldName2, fieldType)

## Execute CalculateField of Lat/Long Fields
expression1 = "{0}".format("!SHAPE.extent.XMax!")
expression2 = "{0}".format("!SHAPE.extent.YMax!")
#expression3 = "{0}".format("!SHAPE.extent@decimaldegrees!")
#expression4 = "{0}".format("!SHAPE.extent@decimaldegrees!")

arcpy.CalculateField_management(inFeature, fieldName1, expression1, "PYTHON_9.3")
arcpy.CalculateField_management(inFeature, fieldName2, expression2, "PYTHON_9.3")
#arcpy.CalculateField_management(inFeature, fieldName1, expression3, "PYTHON_9.3")
#arcpy.CalculateField_management(inFeature, fieldName2, expression4, "PYTHON_9.3")

...
...
...


texts: Having trouble to exportmap in OpenLayer3

texts:
I'm trying to export the map, but I'm always getting the same error:

Failed to execute 'toDataURL' on 'HTMLCanvasElement': Tainted canvases may not be exported.

The funny is that It works the first time... But when I do zoom It appears the error. I have the feeling that The map isn't refreshing or something like this... Really I'm not sure.

I have the following code:

var canvas = event.context.canvas;
var exportPNGElement = document.createElement('a');
exportPNGElement.download = 'Mapa.png';
exportPNGElement.href = canvas.toDataURL('image/png');
document.body.appendChild(exportPNGElement);
exportPNGElement.click();
document.body.removeChild(exportPNGElement);

What is wrong? Do you have any idea?


texts: erase one line layer by another line layer

texts:
I keep on having this issue and I can't find the "perfect answer":

First of all, I work in ArcGIS 9.3.1.
I have 2 layers: roads1950, roads2000. I need to make a layer with roads which were built/changed between 1950-2000. The most perfect thing was to have a "erase-like" tool that would cut out the identical to roads1950 features in roads2000. The problem is that the erase tool only works with polygon layers. I could do a merge and sow out the doubles but it's a (very) long process and it isn't accurate.

Do you know of a tool/method that could help me ?


texts:
Is there a straightforward way of randomly sampling a raster so that the output of the process is a raster?

I'm using an example that I found on the r-sig-geo list and I have also tried the sampleRandom function in the raster package. Both of these approaches produce an output that I am not certain how to transform into a raster. I was not able to find an approach after searching for several combinations of "SpatialPointsDataFrame raster".

library(raster)

# read in raster
rasterSource <- 'landsat.TIF'
r <- raster(rasterSource)

# convert to spatial points data frame
r.spgrd<-as(r,"SpatialPointsDataFrame")

# elminate NA values
r.spgrd = r.spgrd[!is.na(r.spgrd[[1]]),]

# sample points
selectedPoints = sample(1:length(r.spgrd[[1]]), 1000)
r.sampled = r.spgrd[selectedPoints,]

# try to make spgrd into a raster
r.test <- raster(r.sampled)

When I run r.test I get the output:

class       : RasterLayer
dimensions  : 10, 10, 100  (nrow, ncol, ncell)
resolution  : 28617, 14766  (x, y)
extent      : 1838505, 2124675, 2328685, 2476345  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
values      : none

So that the following line which tries to write a raster produces the message:

# write out as ascii file
writeRaster(r.test, filename="test1.ASC", datatype="ascii", overwrite=TRUE)

Error: hasValues(x) is not TRUE

My main objective is to produce some type of raster after the sampling process. I'm also fine with just changing the values within my raster (I'm just not certain how to do that).


texts: Why does QGIS image sometimes appear on an angle?

texts:
Why does the satellite imagery that I bring in and the KMZ file appear to be on an angle, as if it's angled back rather than being a top down view?

Here's a screenshot that should better explain what I'm seeing.


texts:
I am new at QGIS and am having a problem aligning layers. Possibly I just don't understand how projection works.

I have two layers:

GPS waypoints collected in a GPS UTM zone 55N
Polygon shapefile from ArcGIS CRS Guam 1963

The project CRS is EPSG:3857

The background is OpenStreetMap from the OpenLayers plugin and it does align correctly with the GPS UTM layer.

I have tried to save as the layers into what I thought was the correct CRS but I can't get them to line up. Below is the result.

What steps do I need to take to align these two layers?


texts: ArcPy replaceDataSource to ArcSDE raster catalog error?

texts:
I'm trying to re-path a multitude of rasters currently the workspace is a folder, the intended workspace is a raster catalog in an SDE database - using the replaceDataSource method I'm getting the error as below. Does anyone know if this method can't output as a Raster Catalog or what the actual problem may be?

lyr = arcpy.mapping.Layer(r"C:\Folder\MyLayer.lyr")
lyr.replaceDataSource(r"pathToSDEConnectionFile.sde","SDE_WORKSPACE","rasterCatalogName", True)

  Runtime error : Layer: Unexpected error


texts:
Possible Duplicate:
  Problems with OpenLayers and WFS - proxy.cgi issues

I have seen some example related to WMS Get Feature Info like :

http://openlayers.org/dev/examples/getfeatureinfo-control.html

In this example they are using

 OpenLayers.ProxyHost = "proxy.cgi?url=";

My question is :

How we can setup our own proxthost for our geoserver layers?

Help would be appreciated.


texts: Writing shapefile to SQL Server database with ogr2ogr Error

texts:
I have been using ogr2ogr to write point shapefiles to a SQL Server Database of mine.

I used it 2 weeks ago and it worked fine, but, today when I tried to run it, I get the following error:

ERROR1: INSERT command for new feature failed. [Microsoft][SQL Server Native Client 11.0]Memory allocation failure`

ERROR1: Unable to write feature 0 from layer points.`

I used -skipfailures to see what output is produced, and only some of the results are written to a sql table, in this case I have 8 points, and only 5 have been written to a sql table... what could be causing this to happen? Here is the ogr2ogr command I am using,
ogr2ogr -f MSSQLSpatial "MSSQL:server=servername\SQL2012;database=db;UID=username;PWD=password" C:\Pathtoshape\points.shp

I am wondering if I have fatally altered my code by mistake, or is there some other reason and now it is causing this error...

I am working on a solution, and hopefully I or someone else can find the solution, this question hasnt worked for me: ESRI Shapefile not loading into SQL Server


texts:
Do you have any experience with adding custom polygons in to the google maps?

I have made polygons from raster layer and my goal is to show this layer in google maps.

For now i have experience with showing points from mysql in google maps.

I have used raster to vector so they are in shapefile.


texts:
I installed QGIS server and QGIS web client, and I can see

localhost/web/site/qgiswebclient.html?map=/var/www/web/projects/helloworld.qgs&format=image/png&visibleLayers=Hello,Country

And I installed PostgreSQL, PostGIS, and geoserver. i connected import shape files and i see

localhost:8080/geoserver/GIStest/wms?service=WMS&version=1.1.0&request=GetMap&layers=GIStest:Country&styles=&bbox=-176.248494673366,-67.5929956881973,182.893410989597,83.6210855559975&width=783&height=330&srs=EPSG:4326&format=application/openlayers

How to read qgis web client from geoserver wms?


texts: QGIS does not react when using large expression in the Expression Dialog

texts:
QGIS does not respond when pasting a large expression in QGIS with If-statements following each other:
If("A", 'B',
  If("C", 'D',
    If("E", 'F',
      etc)))

I have checked my statement, it is correct and it works in QGIS 3.4 or lower. I have done some testing, and it seems that with each extra bracket the system takes longer to process or evaluate the expression. After about 10 repeating If-statements the waiting time becomes multiple minutes.


texts: Installing plugins for all users and on multiple servers

texts:
I'm running QGIS on a set of high powered terminal server for our students. I want to install a set of plugins for all users on all the servers, making the plugin manager a bit sluggish to use.
Is there a way I can install plugins for all users?

Is there a command line option for the plugin manager?
Is there a folder where I can unpack zip-files so that all users get the plugins on start?


texts:
I use FME to make GML file based on SQLServer data. When I create GML, I cann't see anything in FME Data Inspector or Universal Viewer. Like something is missing in file.
I have schema, and XMLValidator validates file against schema, but I still cann't see graphical data.

Link to Result.gml

Link to schema.

Does anyone knows what is the problem and what is missing in schema?

EDIT:

As requested!

This is just a part, but main part. All features come into XMLTemplater, like on the picture. Some of them (with non-point geometry) goes first in GeometryExtractor and then into XMLTemplater. These StringReplacers works fine. XMLValidator validates document successfully.


texts: Qspatialite not loading in QGIS 2.6.1 for Mac

texts: I am running QGIS 2.6.1 on MacOS (Mavricks). I have reinstalled the Qspatialite plugin, and am still encountering the same issue. The plugin will not open. Has anyone seen/over come this issue?


texts:
I work mainly with ArcGIS in backoffice, but use QGIS for some field data maintenance, the challenge is that:

The info worked in office has about 4.5M records, but for maintenance purposes, we create a small project (with a few thousand records) just with the data needed for the area we're focusing during a given week. So then, when there's new info to add to the main database, we have to be careful with the ID's, some will be duplicated.

Is there a way to ensure the new records added out the field respect the full ID list (of the 4.5 records in backoffice)?


texts: ArcObjects ICommandBar is not refreshing

texts:
I am migrating a 9.3 based Desktop extension, to 10.1
The Extension, is building toolbars on the fly.
It creates the toolbar, and add the commands to it.
This is the snippets of the code that does this

ICommandBar commandBar = Document.CommandBars.Create("name", ESRI.ArcGIS.SystemUI.esriCmdBarType.esriCmdBarTypeToolbar);

commandBar.Dock(....);

int indexAsObject = commandBar.Count;
commandBar.add((ToolContainer) item, ref indexAsObject);
.....
.... (repeat the above two statements to add all commands)
....

ICommandItem toolbarAsCommandItem = (ICommandItem) commandBar;
toolbarAsCommandItem.Refresh();
toolbarAsCommandItem.Tag = tagvalue;

The above code is working on 9.3

But on 10.1, it is adding the commands correctly, but the commandBar is not refreshed, and it is empty
The toolbars are shown empty, and without any name.
The product I am using is ArcMap, and I am calling the above code from

IDocument.NewDocument event

Any idea?