texts:
I've been working on a project, in QGIS 3.16 Hannover (we started the project a while ago, in the Bonn version, but upgraded during the duration of project). Upon project completion we've been asked to hand over the shape files also, not just the layout pdfs.
As it turned out, the last version of the shp is missing geometry of most attributes. The attribute table contains up-to-date state, but most of features now have no geometry, i.e. they exist as entries in the attribute table, but have null geometry.
Does anybody have a solution? Is it possible that this situation arose, because me and my colleague worked on the same project, often simultaneously? I.e. could it have happened because we inadvertently wrote over same shapefiles at the same time?
Also, what is the latest on archiving QGIS projects? Is there an automatic function available?
I guess I'll live with the loss of data in this case, probably it'll take me a day to recover (that is draw missing geometries by hand).


texts: Multiple conditional statements for raster calculation in R

texts:
In the context of a RUSLE erosion modelling I would like to calculate the sum of precipitation at the entry date of phenological phases.

Suppose I have 13 grids, one of which is the phenophase entry date PhenoDate with values from 1 to 356, the other 12 represent average daily shares of the annual precipitation sum for the individual months [dPrecJan,...,dPrecDec].

I would like to calculate the precipitation sum using a set of conditional statements:

If the phenophase begins in january, multiply the phenophase entry date with the january average daily precipitation, else if the phenophase begins in february, multiply [...]

I already tried to implement my conditional statements using the sum of products, where the condition returns either 0 or 1:

PrecSum <- ((PhenoDate<32) * PhenoDate * dPrecJan) + (
(PhenoDate>=32 & PhenoDate<60) * ((31*dPrecJan)+((PhenoDate-31) * dPrecFeb))) + (
(PhenoDate>=60 & PhenoDate<91) * ((31*dPrecJan)+(28*dPrecFeb)+(
(PhenoDate-59) * dPrecMar))) + (( [etc.]))

Actually not too complicated construction, but since I use very high-resolution raster layers, the calculation goes beyond my memory capacities.

I'm looking for a memory-efficient construction, which allows to set multiple conditional statements on my raster calculation.


texts: I am trying to compute the TPI from a raster file (geotiff with elevation data) but when I run the tool it works forever. I have left the PC working over the weekend without any results. The file is 14mb... What can I possibly do wrong?


texts: Reading XYZ with tiles in tif format using QGIS?

texts:
I am trying to set up a XYZ server for scientific data. The server is able to send both PNG and TIF images. The TIF format is used to get numeric values, whereas the PNG is used for visualization purposes.

I would like to be able to interact, within QGIS, with the TIF layer (change color palette, perform some computations in areas, etc.). However, I am not even able to see the TIF layer (works perfectly for the PNG layer).

Is there any way to achieve that?


texts: How to create styled markers with custom marker tooltip

texts:
How to create custom markers and when mouseover on marker I need a tooltip. Can any one suggest to me how to do these type of markers?

I need markers like in the below image:


texts: How do I connect a point in one layer to a polygon in another layer, and have them appear when the point is selected or identified?

texts: In QGIS I am creating a map of nonprofit organizations in Costa Rica, and I would like link each nonprofits headquarters to its zone of impact (not necessarily in the same area), so that when the headquarters is selected or identified, the zone of impact appears as well. I currently have a point layer with all the nonprofit headquarters mapped out, and a polygon layer with all the different areas of impact mapped out. How can a link a specific nonprofit in the headquarter in the point layer to its associated area of impact  polygon?


texts: Export Symbolised Point Feature to Text File in ArcGIS

texts: Is there a way to save symbols of a feature class from Arc Map as XML or text file format? I know one can set up representations on a feature class in a Geoodatabase and Export XML Workspace, but I thought there was an easier way to do this..


texts: Convert Sentinel-5P netCDF O3_TCL (Tropospheric Ozone) subdataset to GeoTiff by GDAL?

texts:
I've used these steps below to extract the subdataset of Sentinel-5P netCDF file:
First, extract the latitude and longitude variables as geolocation arrays
gdal_translate -of VRT HDF5:"S5P_OFFL_L2__AER_AI_20190101T013729_20190101T031859_06309_01_010202_20190107T010005.nc"://PRODUCT/latitude lat.vrt
gdal_translate -of VRT HDF5:"S5P_OFFL_L2__AER_AI_20190101T013729_20190101T031859_06309_01_010202_20190107T010005.nc"://PRODUCT/longitude lon.vrt

And create a customized VRT file for a subdataset
<VRTDataset rasterXSize="450" rasterYSize="3245">

 lon.vrt
 1
 lat.vrt
 1
 0
 0
 1
 1

<SimpleSource>
  <SourceFilename relativeToVRT="1">HDF5:S5P_OFFL_L2__AER_AI_20190101T013729_20190101T031859_06309_01_010202_20190107T010005.nc://PRODUCT/aerosol_index_340_380</SourceFilename>
  <SourceBand>1</SourceBand>
  <SourceProperties RasterXSize="450" RasterYSize="3245" DataType="Float32"/>
  <SrcRect xOff="0" yOff="0" xSize="450" ySize="3245" />
  <DstRect xOff="0" yOff="0" xSize="450" ySize="3245" />
</SimpleSource>

Then, project the subdataset to GeoTIFF
 gdalwarp -geoloc -t_srs EPSG:4326 -srcnodata 9.96921e+36f aerosol_index_340_380.vrt aerosol_index_340_380.tif -tr 0.069 0.069 -tap

The problem is this O3_TCL (Tropospheric Ozone) (e.g: file name: S5P_OFFL_L2__O3_TCL_20190101T120803_20190107T125357_06315_01_010105_20190115T235953.nc) data doesn't have latitude and longitude variables.
Does anyone know how to process this data to GeoTiff?
Below is what it looks like when I open with Panoly software

It has these values for lat and long variables:
float latitude_ccd(latitude_ccd=80);
  :_FillValue = 9.96921E36f; // float
  :long_name = "pixel center latitude for CCD data";
  :units = "degrees_north";
  :standard_name = "latitude";
  :valid_min = -20.0f; // float
  :valid_max = 20.0f; // float
  :_ChunkSizes = 80U; // uint

float latitude_csa(latitude_csa=8);
  :_FillValue = 9.96921E36f; // float
  :long_name = "latitude center for CSA data";
  :units = "degrees_north";
  :standard_name = "latitude";
  :valid_min = -20.0f; // float
  :valid_max = 20.0f; // float
  :_ChunkSizes = 8U; // uint

float longitude_ccd(longitude_ccd=360);
  :_FillValue = 9.96921E36f; // float
  :long_name = "pixel center longitude for CCD data";
  :units = "degrees_east";
  :standard_name = "longitude";
  :valid_min = -180.0f; // float
  :valid_max = 180.0f; // float
  :_ChunkSizes = 360U; // uint

float longitude_csa(longitude_csa=18);
  :_FillValue = 9.96921E36f; // float
  :long_name = "longitude center for CSA data";
  :units = "degrees_east";
  :standard_name = "longitude";
  :valid_min = -180.0f; // float
  :valid_max = 180.0f; // float
  :_ChunkSizes = 18U; // uint


texts: raster::intersect not producing expected results

texts:
I am having trouble with R's raster::intersect method. I give intersect a SpatialPolygonsDataFrame and SpatialPoints as documented in the raster reference manual (p. 111). The result acts like it ignores second parameter. My code follows.

Make spatial points, note pcts_10_sp is a SpatialPolygonsDataFrame

temp_coords <- coordinates(pcts_10_sp)
temp_crs <- crs(pcts_10_sp)
sp2010 <-SpatialPoints(temp_coords , proj4string=temp_crs )
summary(sp2010)

Results of SpatialPointsDataFrame summary

Object of class SpatialPoints
Coordinates:
                min       max
coords.x1 -97.03242 -96.52624
coords.x2  32.55765  32.98299
Is projected: FALSE
proj4string : [+proj=longlat +datum=NAD27+ellps=clrk66+nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat]
Number of points: 694

Read polygon shapefile and project

pcts_12_sp <- rgdal::readOGR(dsn="Precinct12.shp", layer="Precinct12", verbose = TRUE)
pcts_12_sp <- sp::spTransform(pcts_12_sp, CRS("+proj=longlat +datum=NAD27"))
summary(pcts_12_sp)

Results of SpatialPolygonsDataFrame summary

Object of class SpatialPolygonsDataFrame
Coordinates:
        min      max
x -97.03810 -96.5166
y  32.54508  32.9892
Is projected: FALSE
proj4string :
[+proj=longlat +datum=NAD27 +ellps=clrk66
+nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat]
Data attributes:
    2012_PCT         LD
 1000   :  1   100    : 72
 1001   :  1   114    : 68
 1002   :  1   103    : 66
 1003   :  1   108    : 66
 1004   :  1   109    : 61
 1005   :  1   105    : 56
 (Other):768   (Other):385

Run raster:intersect function

RasterShp <- raster::intersect(pcts_12_sp, sp2010)
summary(RasterShp)

Results of intersect summary

Object of class SpatialPolygonsDataFrame
Coordinates:
        min      max
x -97.03810 -96.5166
y  32.54508  32.9892
Is projected: FALSE
proj4string :
[+proj=longlat +datum=NAD27 +ellps=clrk66
+nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat]
Data attributes:
    2012_PCT         LD
 1000   :  1   114    : 57
 1001   :  1   100    : 55
 1003   :  1   109    : 52
 1004   :  1   108    : 49
 1005   :  1   103    : 46
 1006   :  1   111    : 45
 (Other):610   (Other):312

The barebones plot of polygons without any points.

plot(RasterShp)


texts: Greek characters are not displayed correctly in QGIS on Windows 8

texts:
I am working on a laptop with two operation systems. One Windows 8 and one with Ubuntu 12. In both systems I have installed QGIS and in Windows I also have the FME 2014. I am trying to load a shapefile which has attributes in Greek characters but while everything is displayed correctly in my Ubuntu installation, I get "strange" characters in my Windows system (both in FME and QGIS). The encoding of the shapefile is in UTF-8 and I have tested it in another PC with Windows 7 and it was displayed correctly.

I tried to save the shapefile and loading it again but this doesn't fix the problem. I am convinced that this must be related with my Windows installation but I don't know what I can do to fix it?


texts: How to add a Group Layer control to a Table of Contents sample

texts:
I'm still trying to find the perfect Table of Contents sample for the ArcGIS Server JavaScript API 2.3. My previous question gave some hints but I now have a more specific question regarding checkboxes for layer control.

I need to provide a Table of Contents containing a variety of Tiled, Feature and Dynamic layers, which may be sourced from multiple ArcGIS Server instances. The TOC should allow the user to toggle individual layers on/off, and toggle multiple layers on/off simultaneously using a Group layer checkbox (ie, one high-level checkbox to control multiple lower-level layers).

This sample contains a nice Group layer control, but presumes that all the component layers may be switched on/off as VisibleLayers of a single Dynamic layer. This logic seems to be integral to the sample, so it doesn't look like a simple matter to modify it for feature and tiled layers. (Please let me know if I'm on the wrong track, and this sample would prove suitable)

The Legend Widget sample is more suitable as it simply switches any layer on/off (see the Natural Disasters panel on the right-hand side).

My question

How would I add a Group Layer control to the Legend Widget sample? eg, add a higher-level checkbox which will switch both Fire and Earthquakes layers on/off.


texts: How to iterate and combine bands in model builders

texts:
Similar question here:
How to combine more iterations of a for loop to bands of a single raster?

I have downloaded over 100 Landsat MSS and TM images with separate bands, and I want to combine bands as a single file with multi-bands by the composite bands function.
However, I got stuck in using the ArcGIS model builder. As an example for TM images, my idea was to use six submodels (iterate raster) to gather 6 TM bands (except thermal), and combine bands, but each iterate raster function only outputs one raster, instead of multiple rasters.
The files were organized by years. In other words, each year may have MSS and TM bands together, but MSS and TM bands did have specific file names.

Any further suggestion?


texts: How to determine the accuracy of getFeatureInfo on a WMS layer?

texts:
I have added a WMS layer to the ArcGIS Server JSAPI and am attempting to query it using getFeatureInfo when the user clicks on the map.

This is mostly working, except that the accuracy of the results is very suspicious. If I click within a few hundred metres of a polygon boundary, I often get back the neighbouring polygon's results:

In this case I have clicked within Staffordshire Moorlands 010B (as shown by the infoWindow's anchor) but the result shows the neighbouring polygon 011C.

In order to run the getFeatureInfo, I am calculating the map extent in its native web mercator, converting this into the units of the WMS layer, and using this as the bounding box.

Is this a valid approach?

The question ArcGIS Server WMS GetFeatureInfo tolerance implies that there may be some trickery required to get the tolerance right?


texts: Interpolating sea level data using ArcGIS Desktop?

texts:
I have monthly sea level data obtained from PSMSL in the form of excel sheet for the region Chennai.

I want to interpolate the data in ArcGIS Desktop. The shapefile of the region has Polygons and I want to interpolate 26 Points on it to calculate the CVI.


texts: Determine QGIS version in a startup script that PYQGIS_STARTUP points to

texts:
The environment variable (Win10) PYQGIS_STARTUP points to a python script that is executed before QGIS starts (cp. Central deployment of startup.py in QGIS).

Now I need to determine the QGIS version in the python script PYQGIS_STARTUP points to. Of course, this does not work (Answer to Getting user's QGIS version using PyQGIS?):

try:
    from qgis.core import Qgis
except ImportError:
    from qgis.core import QGis as Qgis

if Qgis.QGIS_VERSION_INT >=30000:
    #do something
else:
    #do something else or nothing

resulting in

ModuleNotFoundError: No module named 'qgis'

Is there another way to determine QGIS version in such script or am I lost?


texts: Map not rendering at zoom level 18+ using openlayers3 and leaflet

texts:
I've configured an OSM map server with mod_tile/renderd/postgres to display a basic map and everything works fine until I try to zoom into level 18 or more. At any zoom level less than 18, I can see renderd -f accepting connections and generating new map tiles. However, when I zoom into level 18 or more, I just get 404's in my weblogs and renderd service is completely quiet.

Also this only seems to be a problem through the browser. When I use render_list to pregenerate tiles, it works at every zoom level. And even when I point the browser to the area of the map that's been pregenerated at zoom 18+ the tiles don't display and again I get a bunch of 404's. I've used both openlayers3 and leaflet libraries and the problem is the same. Here's the relevant configs

renderd.conf

[style2]
URI=/osm_tiles2/
TILEDIR=/var/www/mod_tile
TILESIZE=256
XML=/var/www/osm-bright/osmbright.xml
HOST=localhost
TILESIZE=256
MINZOOM=4
MAXZOOM=20

leaflet

<!doctype html>

    <link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet/v0.7.7/leaflet.css" type="text/css">
    <link rel="stylesheet" href="leaflet/L.Control.Zoomslider.css" type="text/css">

    <title>Maps</title>

The openlayers3 html is basically the same. Anybody seen this before or can assist?


texts: Choosing Open GeoSuite package to install on cloud hosting server?

texts:
Our company has asked me to implement a system to share cartography on the Internet through Open Source GIS technology with our operators.

We have acquired the CloudServer hosting, and we have installed GeoServer from the Open GeoSuite, and I have created an html where I have my map viewer / application created using a geoexplorer link and another with a link from the Open Layers application .(attached image)

How do I configure the web server to share this html with other clients through the web?


texts: How to resize a raster SUMMING the values of pixels in QGIS?

texts:
I would like to be able to resize a raster such that the values of the resized pixels are based on the SUMS of the original raster pixels, rather than the averages. The output should be generated thus: each new pixel value equals the sum of any old pixel values contained fully within its "zone" of coverage, plus the corresponding fraction of the values of any pixels which are partially within its "zone" of coverage.

This would be very useful for resizing rasters containing certain types of data for which averaging the values is not necessarily appropriate, and preserving exactly the same total sum of pixel values may be important--population counts, for example. Is there a straightforward way to do this in QGIS?


texts:
I using QGIS on Ubuntu and Windows, I have some questions for QGIS on Ubuntu, first I follow official details to install QGIS on Ubuntu.

On the QGIS on windows I have processing tab to show me toolbox, graphical model and more.

In that toolbox I can use algorithms from SAGA GIS,GDAL GRASS GIS and more.

My problem is on the QGIS Ubuntu I dont have that toolbox and algorithms SAGA, GRASS, GDAL, any idea why ?

Any idea how to install QGIS complete to use that algorithms on Ubuntu?

Finally in the plugins I don't have more plugins like with the Windows version, the plugins on QGIS Ubuntu is specific.


texts:
How to execute for each table in the list one fixed query replacing just %table%  in PostgreSQL (like variables and loops)
sql = 'select * from $table$ where...some_filter'

$table = {tb1, tb2, tb3...}

Somenone could clarify?


texts: Tiles are blank when Using GeoWebCache TMS to get zoom/x/y.png files

texts:
I have a similar requirement to Using Geoserver to serve Tiles in zoom/x/y.png format?, so my existing application works with OSM server and I've been trying to switch it to use our local GeoServer for the tiles.

I've implemented the Y coordinate conversion detailed in the question above, but that has not solved my issue, in that I am being returned empty PNGs from the server. Specifically i'm using this URL:

http://localhost:8080/geoserver/gwc/service/tms/1.0.0/osm:osm/16/32764/43740.png

which I am expecting to be a tile in Isle of Dogs, London UK. This tile from OSM.

I'm using an OSM download deployed in potgis, and previews from the web interface do seem to work so it just seems like using the TMS URL's are going wrong somehow.

The one thing I'm unsure of if is the server will generate the tile if it's not already in the cache, or if it only serves tiles if they already exist?

Is there a way to make it generate them if this is the case?

I've just tried this:

http://localhost:8080/geoserver/gwc/service/tms/1.0.0/osm:osm/1/1/1.png

And it produces this:

So it looks like the layer is correct, I should have mentioned I only have the UK so my guess is that my co-ordinate conversion is not correct, so I'm looking in to that but I thought I had converted Correctly.


texts: Importing South African NGI OSM imagery xml file now unsupported by QGIS

texts:
Here in South Africa, at our very small company, we often make use of this simple xml file (downloadable from http://aerial.openstreetmap.org.za/aerial-osm-org-za-tiles.xml) to link to this imagery: aerial.openstreetmap.org.za  since any other OSM imagery simply doesn't enjoy our archaic internet speeds.

Problem is in January this year suddenly my QGIS (back then version 2.12 i think) started to refuse to load the xml file (I used to find the xml file via the load raster button and wham it would load imagery for a polygon which was already loaded, I never changed the process of loading the imagery). Now I get this error and then it refuses to load the imagery:

  Unsupported Data Source:
  C:\Users\ColinF\Downloads\aerial-osm-org-za-tiles.xml is not a
  supported raster data source.

The interesting thing is that my colleagues are also using the same xml file (just copied onto their PCs) and their QGIS is linking just fine to get the imagery (via the same process). I downloaded a new xml file and still get the same error.

So what am I doing wrong? Or is there a small setting somewhere I have changed on my QGIS?

We use very similar PSs spec wise and highly similar operating systems, one thing is they are using a QGIS 2.8.1 version and I am using 2.14.0 (however recently colleagues have used 2.14.0 with the imagery so I'm baffled). I have also tried reinstalling it back to 2.8.1 without any success.


texts: Setting Halo property of label using ArcPy and CIM

texts:
I don't know how to set halo property in Python. I try it by CIM (code below) but I don't know, what line "lc.textSymbol.symbol.haloSymbol" should contain.
Can you advise me how to achieve halo effect using ArcPy?
aprx = arcpy.mp.ArcGISProject("current")
m = aprx.listMaps("Map15")[0]

show_labels = m.listLayers("okolni_listy_ztm5*")[0]
show_labels.showLabels = True

for lyr in m.listLayers("okolni_listy_ztm5*"):
if lyr.supports("SHOWLABELS"):
if lyr.showLabels:
for lblClass in lyr.listLabelClasses():
lblClass.expression = "$feature.OBJECTID"

label_show = m.listLayers('okolni_listy_ztm5*')[0]
l_cim = label_show.getDefinition('V2')
lc = l_cim.labelClasses[0]

lc.textSymbol.symbol.haloSize = 2
lc.textSymbol.symbol.haloSymbol =

label_show.setDefinition(l_cim)


texts: Converting coordinates from GCS_WGS_84 Linear_Meter to lat/long?

texts:
I have a LiDAR .laz file with points from Toulouse (France). I convert this file to a .las file using LAStools's las2las function. Then, using lasinfo on this .las file gives me some information, including:

GeoKeyDirectoryTag version 1.1.0 number of keys 3
  key 1024 tiff_tag_location 0 count 1 value_offset 2 - GTModelTypeGeoKey: ModelTypeGeographic
  key 2048 tiff_tag_location 0 count 1 value_offset 4326 - GeographicTypeGeoKey: GCS_WGS_84
  key 4099 tiff_tag_location 0 count 1 value_offset 9001 - VerticalUnitsGeoKey: Linear_Meter

So, I assume that my data is in the WGS84 format.

Regarding points, I have this information using lasinfo on my .las file:

scale factor x y z:         0.001 0.001 0.001
offset x y z:               1574926.999698425875977 3157519.458761402405798 142.690597958476445
min x y z:                  1574775.710 3157384.005 125.497
max x y z:                  1574981.625 3157614.866 226.186

Looking, for exemple, at the minimum point which should be located in Toulouse: (1574775.710, 3157384.005) and following the .las documentation, my real point is :

Xcoordinate = (Xrecord * Xscale) + Xoffset
Xcoordinate = (1574775.710 * 0.001) + 1574926.999698425875977 = 1576501.774698425875977
Ycoordinate = (3157384.005 * 0.001) + 3157519.458761402405798 = 3160676.842761402405798

However, as I saw on many sources, WGS_84 should be between [-90, 90] and [-180, 180], right ?

I saw this question on the forum that may be related to my problem, but even playing with the 3second-->100m that he talks about, I can't find any way to go from my data point (1574775.710, 3157384.005) to the (43.603266, 1.442261) latitude/longitude point that Google Maps gives me for Toulouse in France.


texts: How to batch merge rasters using FOSSGIS tools?

texts:
I'm working with QGIS 1.8 (gdal 1.9.2) and GRASS 6.4.3...but on Windows 7.

I "just" want to merge my raster.

I've got a folder with 50 GTiff raster from Lidar data.

The aim is to merge raster before or after an import into GRASS in order to do some contours...

Before GRASS, I tried gdal_merge with QGIS's gui, without any option. I had an error message : "permission denied".

So I tried it in GRASS with r.patch but :
First the gui don't allows you to select multiple raster at one time...

Then I don't know how to batch with GRASS functions in Windows...I tried the GRASS command line but can't locate in the current mapset (didn't understand how). When I had the windows installer for QGIS, this commandline putted you in the current mapset but it's not the case with the basic installation of QGIS.

I'm quite locked...Any help ?


texts:
I need to style an ArcGIS Server 10.0 WMS Service with an SLD File. I have picked a very simple point style SLD which should style all points with one circle.

http://docs.geoserver.org/latest/en/user/styling/sld-introduction.html

I have then simply defined the SLD-Path in ArcGIS Server by editing the appropriate properties. This is however not being rendered. The points are still being symbolized with the symbols defined in the mxd with which it was created. Has anyone any experience with binding SLDs to ArcGIS Server wms services who can possible see where I am going wrong?


texts: Creating polygon features on scanned historical maps to measure areas in ArcMap?

texts:
I have a collection of historical scanned topographic maps from the 19th century that have already been georeferenced (NAD 1927).
I'm uploading the maps as JPGs into ArcMap and would like to use the editing toolbar to create polygon features which I can then measure the area of. When I click "Start Editing" I receive a message that there are no editable layers.

How can I create polygon and line features so I can get measurements? Is my problem uploading as a JPG? Also, will the area measurements even be accurate in this sense? I would like area displayed in meters or feet squared, but will just drawing a polygon around a lake on the map provide a correct measurement (assuming I can get to the point where I can draw a polygon)? I'm relatively new to GIS and feel as if this should be as easy enough task but I am struggling a lot with it.


texts: Is there any way to rotate a polygon or a feature by entering a numeric value in QGIS


texts: How to get map units to find current scale in OpenLayers?

texts:
I need my OpenLayers application to be able to determine the current scale.  This is not for display (there is already the scale line for that).  This is because I have a scale-dependent style function, to display styles differently at different scales.  I need to be able to do it according to scales, not resolutions, because the styles and their respective scales, are configured by the user, and the user enters scale values similar to how they would in a GIS for layer display constraints.  (Ie, although scales are not great for on-screen, they are what users are used to using to constrain styles/layers to specific resolutions).

This is all working well in my app, changing styles as I zoom in and out.  Except that my calculation of the current scale is incorrect.  For testing purposes, I've simply calculating the scale (incorrectly) as:

resolution * DPI

I think that the correct way to calculate the scale would be:

resolution * meters-per-map-unit * inches-per-metre * DPI

The only part of this I cannot figure out is the metres-per-map-unit.  there is the ol.proj.METERS_PER_UNIT lookup, but I can't use it because I don't know what the current unit is.  My users may have swapped in any map of any projection in the application, so I need to be able to determine the current unit dynamically.

How can I get the current map unit from OpenLayers?

Or can I get it from proj4.js for the current projection? (and if so, is that the same thing?)


texts: ArcGIS Desktop 10.4 crashes when opening a 10.3 MXD

texts:
Freshly installed ArcGIS desktop v10.4 over v10.3
ArcMap crashes (sometimes) when opening a 10.3 MXD.
I updated 2 machines from 10.3 to 10.4.  One machine is crashing when trying to open the MXD's we have been using without problems up to now.  The thing is, it only crashes when opening some of the MXD's, not all.

On the 2nd machine...no problems at all.
These machines are identical.

Ideas?


texts: What is the best FME Transformer for Text search and replace in attribute fields

texts:
I have a number of feature classes which have fields containing text which has somehow become improperly encoded.  The fields contain german text which should be either ä, ö, ü, ß but instead contains for instance KindergÃ¤rten instead of Kindergärten ...etc

What transformer best allows me to check the field for "Ã¤" and replaces it with ä?


texts: Select by attributes ModelBuilder - text string partial values for parameter

texts:
I have a feature class of suburbs that I want to do a select by attribute on the Locality (=suburb names) field, passing a user input string parameter to selection expression.

I have the query working when the string parameter is the full Locality name - the working query is UPPER("LOCALITY") = UPPER('%String%'). The Locality string values are title case.

I want the user to be able to pass a string that is only part of the full locality string, as the parameter for the expression. In the data I have localities Mount Jeff, Jeffs Point, Jefferson - I want the user to be able to input "Jeff" as the parameter text string and select these three occurences.

I think my confusion is with how to add the wildcard to the UPPER('%String%') to indicate that the parameter string may only be part of the text field value. My previous use of  wildcards has been using the % character but cannot get my head around this expression where the % character is also used to identify the parameter.

Currently at home working off a shapefile copy of the data, (finished model will be pointing to an SDE/enterprise feature class).


texts: Clear task list in the Earth Engine command line tool

texts:
I would like to reduce the output which is printed to the terminal when I run:
earthengine task list

Is there some way to clear the task list in the earthengine command line tool?
UPDATE
I don't know whether the earthengine task list can be cleared. However, the output which is printed to the terminal can easily be reduced.
This command will only print the first 5 tasks:
earthengine task list | head -n 5;
while this command will print the last 5 tasks:
earthengine task list | tail -n 5


texts: How to extract coordinate pairs of crossing points of two polygons in PostGIS?

texts:
Extracting coordinate pairs of crossing points of two polygons in PostGIS

Is it possible to extract x,y-coordinate pairs of crossing points of two overlapping polygons in PostGIS?

Suppose there are only 2 overlapping points of the polygons the ST_Intersection function extracts more than two coordinate pairs (also the shared vertices of polygons are included).

SELECT ST_Intersection(geometry geomA , geometry geomB);

I am interested only in the 2 boundary crossing points of polygons.


texts:
I posted this on the OSM help page but it was suggested I try here instead:

I'm new to OSM data and I'm using the ArcGIS editor for OpenStreetMap within ArcMap. The data downloads correctly but the table doesn't seem to include the 'name' field (which is the field I really need).

When I follow the instructions listed here, I receive "ERROR 0008000: The value is not a member of ALL." It works when 'ALL' is there but if I try and add 'name' it fails.

Any suggestions?


texts: clarification of "Section a line in Arcmap from a known point"

texts:
Prior to my previous question: https://gis.stackexchange.com/questions/217031/section-a-line-in-arcmap-from-a-known-point

This point represents where a gas pipeline crosses a river (or road, rail for other points). The length of this crossing is given as 8 metres. The other point crossings vary from 1 - 100metres.

I want to split the lines so that lengths of crossings are shown as new line sections. So from this point, a new line section would be created on the gas pipeline, 4 metres long either side of this point.

I would like to do it automated as possible. I have 60 points to work through, all with varied crossing lengths.

Original question in full:

I have a line layer of pipes and a point layer of where these pipes cross roads, rivers etc. I have already snapped the points to the lines. In the point layer I have a length of these "crossings" e.g a bridge.

I want to be able to "re-section" the pipe line layer so that from the point a bridge of 10 metres would have a new section 5 metres either side of the point of crossing.

I am using Arcmap 10.3 and DO NOT have access to Spatial Analyst.


texts:
With PostGIS, I need to perform the union of two polygon layers that have intersection areas and distinct areas while preserving all limits. I need to extract separated polygons for each intersected and distinct areas from both layers.

Using QGIS, the tool is called "Vector > Geoprocessing > Union"

here is a visual of the problem:

Two Input polygon layers (green vs red dashes) =

Output polygon layer (one distinct color per polygon)=

I need to do it with PostGIS. st_union seems to do a different job. Any idea?


texts: Creating hundreds of maps showing different layers but same area using QGIS?

texts: I want to create bird species location maps for a single area (and there are multiple areas). There are a lot of species, every one needs an individual map, so I'm facing 100+ maps per area. Is there any way beyond manually creating the same map extent per every species? I imagine i can deal with visibility of specific species with style based on $id = @atlasfeatureid, but I'd hate to copy-paste the same rectangle over 100 times for each area.


texts: I have a massive amount of tiff files (+prj. all together 100k). For test purpose I have loaded a small bit of these into geoserver and worked fine. But when I try to add to the store the data folder what contains all of the images, it is just loading nothing else happening. Is there any way to load these amount of rasters into geoserver?


texts: How to find names of Unincorporated parts of a county in Navteq data?

texts:
I am currently using Navteq Data and am unable to find how to label unincorporated places within the LA County. Most of these spots are small and spread throughout the entire county however they do have well known names, most times these are same as neighboring cities.

Has any one had similar problem? i.e. how to label the unincorporated parts of any county?

Thanks,
Shakira


texts: Clipping but using outer border of grouped polygons only

texts:
I am buffering, dissolving, and then clipping a 100 mile radius around each POI over a US map. I got the map from the US Census, and it represents every Zip Code Tabulation Area. Unfortunately, there are no ZCTAs in some parts of the country (anyone know why that is?). So, when I add my POI and then buffer, dissolve and clip, I get a splotchy result, and there are inexplicable lines running through the radii:

Here's a trial I did using the county shapefile. The entire map is covered in counties, so I didn't have the above problem here. How can I get this kind of result for the ZCTA shapefile?


texts: DEM file in .tiff format but Esri Grid Format (.adf) is needed

texts:
The DEM that I received is not in .adf format; it is in .tiff format.
How do I get that supporting file using .tiff.
I am using ArcMap.


texts: Converting ArcMap feature layers to GeoPackage using

texts:
Is it possible to directly convert feature layers in ArcGIS to GeoPackage?

In ArcGIS Desktop 10.5's ArcToolbox, under "Conversion Tools/To Geopackage" the only option that appears to me is "Add Raster to GeoPackage". There is no "Add Features" or something alike.

In Support for OGC GeoPackage specification in ArcGIS, a blog post written in 2014, Esri claimed that "At 10.2.1 or with 10.2.2 ArcGIS desktop, you can create an empty GeoPackage and populate the GeoPackage by copying feature data into it". However, I cannot see how. In ArcCatalog, there is no option to create a GeoPackage file.


texts: Importing 3D geometries (Collada DAE) into PostGIS?

texts:
I have been trying to import some 3D geometries (from Blender, Meshlab or Sketchup ) into PostgreSQL Database (with PostGIS extension). I found this similar unanswered question - Importing 3D geometries from DXF, VRML or X3D into PostGIS? - which led me to this PostGIS & 3D blog post.

I found this method to export X3D files to postGIS but I get this error:

Error : Command not found

I have a sample file in COLLADA (COLLAborative Design Activity)'s format (.dae) extension: https://transferxl.com/08vMDqqmtsWg6N

How can I import 3D geometries, specifically Collada DAE format, into PostGIS?


texts: How to Extract Data from U.S. National Landcover Dataset (NLCD) dataset?

texts:
I am working with U.S. National Landcover Dataset (NLCD) to classify habitat type at more than 150 sites across the northeast U.S. The dataset is very large (15GB) so I cannot upload it here, but it comes in .img format at 30m resolution. I have GPS coordinates in decimal degrees  (WGS1984) for the center point of all the sites. I would like to be able to extract the proportion of landcover classes in a 1 square kilometer around the point. My questions are:

How do I upload .img files into ArcGIS?
How do I extract the information from around the GPS coordinates as proportions of the different habitat classes?
How do I project my sites onto the NLCD map?
Should I change the projection of my coordinates?


texts: How to switch between Advanced and simplified list of algorithim in QGIS 2.14.1?

texts: There is no option to switch between advanced Advanced and simplified  list of algorithm in this version 2.14.1 at the bottom of processing tool


texts: Measuring Continuous Distance of GPS track in ArcGIS Desktop

texts:
I am looking for a way to quickly measure the distance covered of GPS tracks of foraging Seabirds. A tool that could give accumulative distance of point to point from the long and lat coordinates would be ideal.

I have tried a few different tools in the distance toolbox but they don't seem to measure continuous distance from point to point more from a single source. And the Hawth's tools is not supported on arcmap 10.1, (http://www.spatialecology.com/htools/pntdistmatrix.php) which sheds some light on potential tools.

Thanks in Advance


texts:
This is an issue I've been trying to solve for a while and am told it's not possible.

I have made a PDF form with an input box for coordinates. As the form is for use on mobile devices, I would rather have a button to "grab" the coordinates directly from the device.

Within the PDF there is an area to insert Javascript but the script that works for HTML or some other apps does not seem to work.  Importing a map from ArcMap with geospatial data isn't the answer.


texts: Geoserver SLD Line styling with external graphic instead of WellKnownName

texts:
I have problem with styling in geoserver.
Using tutorial from https://docs.geoserver.org/stable/en/user/styling/sld/cookbook/lines.html#id7
everything is ok but when I want to use ExternalGraphic instead of WKN symbol nothing happend.

My sld file look like:

<?xml version="1.0" encoding="UTF-8"?>
<sld:UserStyle xmlns:sld="http://www.opengis.net/sld" xmlns="http://www.opengis.net/sld" xmlns:gml="http://www.opengis.net/gml" xmlns:ogc="http://www.opengis.net/ogc">
    <sld:Name>Default Styler</sld:Name>
    <sld:FeatureTypeStyle>
        <sld:Name>name</sld:Name>
        <sld:Rule>
            <sld:LineSymbolizer>
                <sld:Stroke>
                    <sld:GraphicStroke>
                        <sld:Graphic>
                            <sld:ExternalGraphic>
                                <sld:OnlineResource xmlns:xlink="http://www.w3.org/1999/xlink" xlink:type="simple" xlink:href="http://127.0.0.1/svg2.svg" />
                                <sld:Format>image/svg</sld:Format>
                            </sld:ExternalGraphic>
                            <sld:Size>16</sld:Size>
                        </sld:Graphic>
                    </sld:GraphicStroke>
                </sld:Stroke>
            </sld:LineSymbolizer>
        </sld:Rule>
    </sld:FeatureTypeStyle>
</sld:UserStyle>

Someone of you know why it happens?

BTW I try make it work because I need line on my wms looking like
http://wms.e-mapa.net/cgi-bin/mapserv?map=/home/www/egmina/implementation/czerwonak/pln/wms/mpzp.map&VERSION=1.1.1&SERVICE=WMS&REQUEST=GetMap&LAYERS=lzb&SRS=EPSG:2180&BBOX=364177.88820526295,519023.6162198745,364271.63820526295,519054.7197354995&WIDTH=1920&HEIGHT=637&TRANSPARENT=TRUE&FORMAT=image/png&styles=,,,,,&unique=3971

Thanks for any help.


texts: Initiate field value before adding new feature

texts:
I have two layers:

Chambre(obj_id)
Ouvrage(obj_id, fk_chambre)

The attribute Ouvrage.fk_chambre is referencing the attribute Chambre.obj_id.

I'm trying to create an action on the Chambre layer that permit to add a feature in the Ouvrage layer and directly add the obj_id of the chambre feature in the new created Ouvrage feature.

Here is my actual code:

ouvrage_lyr = QgsMapLayerRegistry.instance().mapLayer('od_ouvrage20171121140457844')
chambre_lyr = QgsMapLayerRegistry.instance().mapLayer('od_chambre20171121140457773')

feature = chambre_lyr.getFeatures( QgsFeatureRequest().setFilterExpression( 'obj_id = \'{}\''.format( '[% "obj_id" %]' ) ) ).next()

#Get the value of the clicked feature obj_id
fk_chambre = feature["obj_id"]

#create the new feature in the ouvrage layer
ouvrage_lyr.startEditing()
qgis.utils.iface.setActiveLayer(ouvrage_lyr)
qgis.utils.iface.actionAddFeature().trigger()

My question: How can I automatically fill the attribute Ouvrage.fk_chambre with the fk_chambre variable?


texts: Difference in GeoTiffs QGIS Export vs GDAL Warp

texts:
If I export a GeoTiff from PostGIS (which is essentially a gdal operation), it loads into QGIS using the VRT driver but works fine.

When I right click and export without making any changes it imports back into QGIS using GTiff driver and this file is 0.6kb larger.

Gdalinfo gives exactly the same output for both files.

The reason for the question is when I send the file from PostGIS to someone they cannot load it into their old software, but no worries with the one exported from QGIS.

Is there some extra header information that gdalinfo does not show? I want to be able to generate a GeoTiff file similar to the exported type from QGIS.

EDIT: SQL:

SELECT lo_from_bytea(1234,(
WITH
rast as ...
SELECT ST_AsGDALRaster((SELECT * FROM rast), 'GTiff', ARRAY['PROFILE=GeoTIFF','TFW=YES']) as tif))

SELECT lo_export(1234,'/tmp/tif.tif')

Driver: GTiff/GeoTIFF
Files: RGB_2017-2-24 (5)__.tif
Size is 512, 160
Coordinate System is:
PROJCS["WGS 84 / UTM zone 55N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",147],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32655"]]
Origin = (815896.631598158506677,-3239186.335398061200976)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  815896.632,-3239186.335) (150d14'59.93"E, 29d14'31.50"S)
Lower Left  (  815896.632,-3240786.335) (150d15' 1.57"E, 29d15'23.40"S)
Upper Right (  821016.632,-3239186.335) (150d18' 9.30"E, 29d14'26.85"S)
Lower Right (  821016.632,-3240786.335) (150d18'10.97"E, 29d15'18.75"S)
Center      (  818456.632,-3239986.335) (150d16'35.44"E, 29d14'55.13"S)
Band 1 Block=512x5 Type=Byte, ColorInterp=Red
  NoData Value=0
Band 2 Block=512x5 Type=Byte, ColorInterp=Green
  NoData Value=0
Band 3 Block=512x5 Type=Byte, ColorInterp=Blue
  NoData Value=0


texts: Problematic file size after georeferencing in QGIS

texts:
I have to georeference a few hundred JPEG maps in QGIS. Each of these maps has an average size of 20 MB. After georeferencing, these maps must be loaded into ArcGIS on another computer, of course correctly georeferenced.
I just have access to QGIS, while my colleague, who has to work with these maps later, has only access to ArcGIS.

Up to now, I tried several attempts in QGIS. Georeferencing itself isn't problematic, but the outcoming filesizes are too big for each single map (mostly greater than 100 MB, sometimes up to 300 MB). While georeferencing I tried the compression options "deflate", "lzw" and "packbits", but all resulting files are still too big for handling them on my colleque's machine.The output files are all in GeoTIFF format.

I tried to georeference the maps without producing a GeoTIFF but only a wld file. After doing so, the map could not be loaded into ArcGIS without showing the warning "file has no spatial reference" even if the original JPEG map and the world-file are stored in the same location on my colleague's PC.

A quick search on this site did not lead to satisfactory results.

Is there a way in which I can georeference JPEG files (20 Mb) in QGIS in a way, that the resulting file shows similar sizes as the original file and can be correctly displayed in ArcGIS on another computer afterwards?


texts: Adding column with area values in attributes table with standalone PyQGIS

texts:
I would like to add an area column to the Attributes table using Python. I prefer not to use QGIS Python console. I would like to do that using cmd. How should I change my code?
My code and error below:
  import sys, os
  sys.path.append(r'C:/PROGRA~1/QGIS3~1.16/apps/qgis/python/plugins')
  os.environ['PROJ_LIB'] = 'C:\\Program Files\\QGIS 3.16\\share\\proj'

  from qgis.core import *
  from qgis.utils import *
  from qgis.analysis import QgsNativeAlgorithms
  from qgis.core import QVariant

  qgispath = r'C:/PROGRA~1/QGIS3~1.16/apps/qgis'

  QgsApplication.setPrefixPath(qgispath, True)
  qgs = QgsApplication([], False)
  qgs.initQgis()

  import processing
  from processing.core.Processing import Processing

  Processing.initialize()
  QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
  feedback = QgsProcessingFeedback()

  processing.run("native:joinattributesbylocation", {'INPUT':r'C:/Users/ja/Inne/Desktop/malebudynki.shp|layername=malebudynki',
                                                     'JOIN':r'C:/Users/ja/Inne/Desktop/adresy/adresyzukladem.shp',
                                                     'PREDICATE':[1,5],
                                                     'JOIN_FIELDS':['PNA','SIMC_id','SIMC_nazwa','ULIC_nazwa','Numer'],
                                                     'METHOD':0,
                                                     'DISCARD_NONMATCHING':True,
                                                     'PREFIX':'',
                                                     'OUTPUT':r'C:/Users/ja/Inne/Desktop/treningowawarstwa.shp'})

  #add area column to attributes table
  layer = r'C:/Users/ja/Inne/Desktop/treningowawarstwa.shp'
  provider = layer.dataProvider()
  area_field = QgsField("Area", QVariant.Int)
  provider.addAttributes([area_field])
  layer.updateFields()

  idx = provider.fieldNameIndex('Area')
  for feature in layer.getFeatures():
      attrs = {idx : int(feature.geometry().area())}
      layer.dataProvider().changeAttributeValues({feature.id() : attrs})

>>> #add area column to attributes table
... layer = r'C:/Users/ja/Inne/Desktop/treningowawarstwa.shp'

>>> provider = layer.dataProvider()

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: 'str' object has no attribute 'dataProvider'

>>> area_field = QgsField("Area", QVariant.Int)

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
NameError: name 'QVariant' is not defined

>>> provider.addAttributes([area_field])

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
NameError: name 'provider' is not defined

>>> layer.updateFields()

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: 'str' object has no attribute 'updateFields'

>>> idx = provider.fieldNameIndex('Area')

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
NameError: name 'provider' is not defined

>>> for feature in layer.getFeatures():
    ... attrs = {idx : int(feature.geometry().area())}
    ... layer.dataProvider().changeAttributeValues({feature.id() : attrs})


texts: GeoServer TMS with leaflet coordinates are different

texts:
I use the example tiger-ny layer from geoserver with leaflet like this:

var map = L.map('map').setView([62.92578, 32.08523], 13);
L.tileLayer('http://localhost:8080/geoserver/gwc/service/tms/1.0.0/tiger-ny/{z}/{x}/{y}.png', {
  maxZoom: 18,
  tms: true,
  crs: L.CRS.EPSG4326,
  attribution: false
}).addTo(map);

I can see the Map and all works but not with the correct latlng.
The coordinates latlng are completely different.
This map shows basically the central park in new york and the correct latlng should be (40.78546, -73.96583) but in this geoserver tms the central park is at (62.92578, 32.08523).

Leaflet with geoserver tms should show the location with the correct latlng, like the central park at (40.78546, -73.96583) and not at (62.92578, 32.08523).

I don't understand why the coordinates are wrong when i use the tms of the GeoServer?


texts: Controlling field names with arcpy.Combine?

texts:
If I run the arcpy.Combine command as such:

result = arcpy.Combine(["CDL/raster07", "NHD/raster07", "OtherRaster"])

The resulting raster has an attribute table with fields

[f.name for f in arcpy.ListFields(result)
['OTHERRASTER', 'RASTER0701', 'RASTER0702']

The field order doesn't seem to correspond to the order in which the fields were input into the tool and they aren't descriptive in a way that I can automatically parse in a script. I'm also unable to rename the input rasters.  Is there a way to control the field names, or to control the order in which they are written?


texts:
I am getting these errors when someone attempts to open my app in IE. I get a popup asking to Allow Blocked Content, and when this is clicked everything works fine. Is there anyway I can go about enabling CORS to bypass clicking Allow Blocked Content?

SEC7118: XMLHttpRequest for http://services.arcgisonline.com/ArcGIS/rest/info?f=json required Cross Origin Resource Sharing (CORS).
EDDM_API.html

SEC7120: Origin file: not found in Access-Control-Allow-Origin header.
EDDM_API.html

SCRIPT7002: XMLHttpRequest: Network Error 0x80070005, Access is denied.
EDDM_API.html

SEC7118: XMLHttpRequest for http://static.arcgis.com/attribution/World_Street_Map?f=json required Cross Origin Resource Sharing (CORS).
EDDM_API.html

Searching the web I found this bit of code, just unaware of where to put it in my script, or if it will even solve my problem.

// Script running on http://eddmapi-dev.usps.gov/eddmapi/eddm_api.html
var xhr = new XMLHttpRequest();
    xhr.onerror = _handleError;
    xhr.onload = _handleLoad;
    xhr.open("GET", "http://eddmapi-dev.usps.gov/eddmapi/eddm_api.html", true);
    xhr.send();


texts: Reclassify NDVI and count number of pixels above certain threshold in Google Earth Engine

texts:
Following on from a previous question (NDVI reclassify with mean +1SD as threshold in Google Earth Engine) I have managed to show the pixels that have value higher than NDVI mean+1SD at the estuary area on the map, which also overlap with the vegetation area. But I do notice the coloured regions are not universally black, which means I failed to do NDVI reclassify. I meant to use NDVI mean +1SD as a minimum threshold for reclassification, above the threshold to value 1 reclassify as 1 and below the threshold to value -1 reclassify as 0.
The purpose of doing ndvi reclassify is :1. test if the pixels above the threshold overlap with the estuary vegetation area (which I have achieved); 2. count total number of pixels above and below the threshold.
The code I have used:

var mean = ee.Number(stats.get('nd_mean'))
var std = ee.Number(stats.get('nd_stdDev'))
var maskImage = ndvi.updateMask(ndvi.gt(mean.add(std)));
Map.addLayer(maskImage, {}, 'mask image');
//count number of pixels above the mean+1SD
var count1 = maskImage.reduceRegion({
  reducer: ee.Reducer.count(),
  geometry: table
  });
print(count1);
//count number of pixels below mean+1SD
var unmaskImage = ndvi.updateMask(ndvi.lte(mean.add(std)));
var count2 = unmaskImage.reduceRegion({
  reducer: ee.Reducer.count(),
  geometry:table
});
print(count2);
//count total number of pixels
var count3 = ndvi.reduceRegion({
  reducer: ee.Reducer.count(),
  geometry:table
  });
print(count3);

I have got number of pixels above the mean+1SD:4274; number of pixels below mean+1SD: 52773;total number of pixels: 57047. I have got the number of pixels in ArcGIS for the satellite image in the same day with the same shapefile for the estuary but the total number of pixel is 57060. Why the total number of pixels calculated by GEE and ArcGIS is different?


texts:
I'm trying to work an example from the "PostGIS Cookbook" and I'm getting an error I don't understand.  I'm running the example using Postgres.app on a MacBook.  Postgres.app is running PostgreSQL 9.1.12 with PostGIS 1.5.3.

Here's the PL/Python function:

-- load.sql
CREATE OR REPLACE FUNCTION Geocode(address text)
    RETURNS geometry(Point, 4326)
AS $$
    from geopy import geocoders
    g = geocoders.GoogleV3()
    place, (lat, lng) = g.geocode(address)
    plpy.info('Geocoded %s for the address: %s' % (place, address))
    plpy.info('Longitude is %s, Latitude is %s.' % (lng, lat))
    plpy.info("SELECT ST_GeomFromText('POINT(%s %s)', 4326)" % (lng, lat))
    result = plpy.execute("SELECT ST_GeomFromText('POINT(%s %s)', 4326) AS point_geocded" % (lng, lat))
    geometry = result[0]["point_geocoded"]
    return geometry
$$ LANGUAGE plpythonu;

When I try to load the function into the database with this command:

psql geodjango < load.sql

I get this error:

ERROR: type modifier is not allowed for type "geometry"

I suspect this error has something to do with the fact that the book example assumes you're using PostGIS 2.0 and I'm using 1.5.3.  However, it still doesn't make sense to me.  I'm read about the ST_GeomFromText function in the PostGIS 1.5 Manual and this code appears to be correct.  What am I doing wrong?

Thanks.


texts: Update Attachments in a File Geodatabase dynamically

texts:
I am working in ArcGIS 10.5 in a file geodatabase.

I have a polygon feature class that currently has pdf attachments for various features.

Periodically, the folder that contains the pdfs gets updated, some features now have multiple attachments. I am trying to avoid manually adding new attachments using the Attachment Manager because sometimes there are hundreds of updates.

Currently, to batch update all the attachments within the GDB, I need to delete the existing attachments, regenerate a match table, and re-run the Add Attachments tool.

I was considering building a model to take care of this workflow but was wondering if there is more efficient way to handle this without having to run a tool or series of tools every time there are changes to the pdf folder?


texts:
I would like to know if anyone has worked with stereoscopic analysis (e.g extracting features such as drainage, forest an so on) using the following softwares: QGis, gvSIG, ILWIS or MicroImages TN Tmips.

If so, how can I do that?


texts: Using a URL instead of a form to add a shapefile to a layer in JavaScript

texts:
I have the following code that is used to add a shapefile to the map and it works perfectly.

var genRequest = request({
                    url: portalUrl + '/sharing/rest/content/features/generate',
                    content: myContent,
                    form: dom.byId('uploadForm2'),
                    handleAs: 'json',
                    load: lang.hitch(this, function (response) {
                        if (response.error) {
                            return;
                        }
                        addShapefileToMap(response.featureCollection);
                    }),
                });

However currently it is using a form where the user has to manually select a file from on his pc.
Is there any way so that instead of using the form I can pass a link which is passed through the code? I can't seem to code the form dynamically to automatically accept the link to avoid the user having the choose the file either.


texts: 3D Symbology for a Model Farm in ArcScene

texts: I am working on visualization of a model farm I have developed consisting of various types of grazing land, cropping areas and other. As I have constructed it myself, I don't have an image with it for display. I am looking for symbology options I can use in ArcScene to give it a realistic touch. I have tried to use picture but it doesn't work well. Here, the screenshots below shows the things I've tried with 3D texture fill sumbol and picture fill symbol. I don't have many options in the texture one and it looks all the same if I use the grass1 or grass2 as symbols for grazing or cropping land. As this is a created farm, not an actual one on earth, so I obviously don't have any aerial images or satellite images for the farm area, which would have been great, therefore, I am not sure how I can make it look real. I want to make a fly over animation from this once I make it look more real. I have draped a DEM on the farm landscape in ArcScene. Any ideas?


texts: ESRI QueryLayer + Javascript API - map displays inconsistent results

texts:
We recently migrated one of our map services from FGDB to an SQL Server spatial db backend. We are noticing a problem now in that the map appears to display locations inconsistently. I created a JSFiddle to demonstrate the issue and to simplify the code as much as possible.

http://jsfiddle.net/acurtiss/pbuLot1g/24/

If you zoom in or out, or simply re-run the application, you should see that some zip codes will randomly go blank. (The zip codes on this map should be nearly 100% filled in.) I have checked the JSON object returned from the map service and verified that the results are consistent and correct, so I'm beginning to suspect a problem with the Javascript API itself. I have contacted ESRI support and have a ticket open with them, but so far they haven't identified the source of the problem so I thought I would open this up to a broader set of GIS experts here on Stack Exchange.

As requested, here is a screenshot showing a significant number of blank zip codes. (Due to some limitations with my stack exchange account, apparently I can't post multiple screenshots, but hopefully this example shows the problem I'm trying to address.)


texts:
I'm working in AutoCAD with a layer from a shapefile.

The problem occurs when Autocad opens the file : the coordinates are the same than in QGIS, but when I try to draw a line 1 meter is 0,00001 in autocad.

I do not know how to get it in a workable scale. I saw similar questions, and tried to apply the solutions shown (modify CRS,...), but them didn't work. I am working with a coordinate system (WGS84)


texts: Is adding a new table for every layer a good design idea?

texts:
Is it a good design idea to maintain separate tables for each layer in a PostGIS database?. Every layer is defined by users. A user creates his/her layer say ATM location. That's meant by defined by users (I don't have control over the table definition of the layers much) This makes it impossible to predict fields(attribute fields and type are not known initially) used in the table/layers. School vs. roads have different geometry and set of attributes.

The problem is traditional databases won't advise having dynamic tables very much. Does the same apply to spatial databases?


texts: Drawing polygon around a raster (only where pixels have data) QGIS

texts:
I have a raster (grey scale image with white background [can change to e.g. green]).
I want to make a polygon of the extent of the pixels with data (outline)- not the whole raster file.
How can I do this in QGIS?
Also can you do this for many rasters in a group (e.g. a mosaic of rasters making a larger image)


texts: Determining Wet Well Service Area using ArcGIS for Desktop?

texts:
I am using ArcMap 10.1

I have a wet well with a known ground elevation. I am trying to determine the service area of that area based on the assumptions that the wet well is 20 feet deep and that my minimum slope is .5%. I do not know which way my pipe will go so I want to create a buffer service area with the possible distance in all directions.


texts: How to set bounding box in geoserver WCS request

texts:
When making the below geoserver WCS request it seems that the BoundingBox param is ignored: &BoundingBox=-86,20,-76,40,urn:ogc:def:crs:EPSG::4269

Is my request written wrong or is BoundingBox not supported with WCS (I do see mention of it in the getcapabilities)? For WMS &bbox=-86,20,-76,40 works.

Here's the full WCS request but it's returning everything in the full extent of the raster (not just the part in the bounding box):

https://geoserver-dev.usanpn.org/geoserver/wcs?service=WCS&version=2.0.1&request=GetCoverage&coverageId=si-x:average_leaf_best&BoundingBox=-86,20,-76,40,urn:ogc:def:crs:EPSG::4269&SUBSET=time(%222003-01-01T00:00:00.000Z%22)&format=geotiff


texts: saga_gui.exe identified as malware by antivirus when downloaded from OSGeo4W

texts:
I downloaded QGIS 2.6.1-2 (Brighton) yesterday from OSGeo4W-setup-x86.exe using the advanced installation option on a machine running on windows 7 professional. The download worked Ok except for saga 2.1.2-1 which I was really looking for. At first the saga package was identified as malware and its download was blocked by my free Antivirus system (Avast). When I decided to disable the Antivirus for a few minutes to allow for the download it then worked fine. However, I felt uneasy about the whole thing so to check, I have run an antivirus scan and it identifies saga_gui.exe as Win32:Malware-gen. If instead I choose saga 2.0.8-5, I do not have this problem but I believe I need the 2.1 version to do add polygon attributes to points.
Has anybody else experienced these issues? Is there a way around this?


texts: Google Earth Engine: Javascript to Python conversion: syntax error with `image1.and(image2)`

texts:
I need to convert my Earth Engine javascript code into the Python API.
Most of it works, but when using the following in Python (2.7):

mask = image1.and(image2)

I get:

File "<ipython-input-65-e84ccdaf9093>", line 41
    mask = image1.and(image2)
                    ^
SyntaxError: invalid syntax

Though it works well with javascript. Any idea how to make this work?
More generally, how to know which methods available in GEE docs are/are not implemented in the Python API?


texts: Intersect maps in Rstudio similar to what is done in QuantumGIS

texts:
I want to make the intersection of two shapefiles in Rstudio as I usually do in QuantumGIS (due to the intersection function that belongs to the vector tools of this software)

I tried to use the left join tools, but in the table of attributes of these shapefiles there is no field in common to make the intersection.

I can do this for QGIS, but I would like to do it for Rstudio, do you have any suggestions?


texts: How to calculate the distance of points to line in R?

texts:
I need the distance from each point of the object found to the transect line.
'id' is a research transect, and 'obj.id' is the object found. in this case, 1 object was found in transect A and 2 objects in transect B.
So how do I calculate the distance between objects towards its transect lines?
I have only come this far:

points <- data.frame(
  'id' = c('A','B','B'),
  'obj.id' = c('x1', 'x2', 'x3'),
  'x' = c(60, 65, 66),
  'y' = c(130, 134, 135)
)

lines <- data.frame(
  'id' = c('A','B'),
  'x1' = c(58, 64),
  'y1' = c(128, 132),
  'x2' = c(61, 67),
  'y2' = c(131, 135)
)

rows <- split(lines, seq(nrow(lines)))
lines <- lapply(rows, function(row) {
  lmat <- matrix(unlist(row[2:5]), ncol = 2, byrow = TRUE)
  st_linestring(lmat)
})

lines <- st_sfc(lines)
lines <- st_sf('id'=c('A','B'),'geometry'=lines)

join <- inner_join(points, lines, by='id')

I tried dist2Line, st_distance with sf packages, but haven't found the right method.


texts: Map is broken on GeoEXT when dir=rtl in chrome (pics inside)

texts:
I've created a GeoExt app with simple map and east panel.
When I change the dir to rtl it works fine in FF but in Chrome the map isn't scrollable.

LTR:

RTL:

Here is the link to my site (with the bug)

http://77.235.53.170/LehavimExt/

I created a JSFiddle example and in the fiddle window it works fine.
http://jsfiddle.net/uvY3r/

So now I'm clueless what is the problem.

This is the code I use :

 var map;
var format = 'image/png';
var epsg4326 = new OpenLayers.Projection('EPSG:4326');
var epsg900913 = new OpenLayers.Projection('EPSG:900913');

Ext.Loader.setConfig({
    enabled: true,
    disableCaching: false,
    paths: {
        GeoExt: "js/GeoExt",
        Ext: "js/ext-4.2.1.883/src"
    }
});

Ext.require([
    'Ext.container.Viewport',
    'Ext.state.Manager',
    'Ext.window.MessageBox',
    'GeoExt.panel.Map'
]);

Ext.application({
    name: 'Lehavim Ext',
    launch: function () {

        var extent = new OpenLayers.Bounds(34.79, 31.36, 34.834, 31.38).transform(epsg4326, epsg900913);
        map = new OpenLayers.Map('map', {
                units: 'm',
                restrictedExtent: extent,
                controls: [
                new OpenLayers.Control.PanZoomBar(),
                new OpenLayers.Control.Navigation(),
                new OpenLayers.Control.ScaleLine(),
                new OpenLayers.Control.MousePosition()  ,
                new OpenLayers.Control.LayerSwitcher()
            ],
        projection: epsg900913,
        displayProjection: epsg4326 //Is used for displaying coordinates in appropriate CRS by MousePosition control
        });

        var resolutions = OpenLayers.Layer.Bing.prototype.serverResolutions.slice(16, 20);
        var osm = new OpenLayers.Layer.OSM("OSM", null, { zoomOffset: 16, resolutions: resolutions });

        map.addLayers([osm]);
        var pos = new OpenLayers.LonLat(34.8161900, 31.370).transform(epsg4326, map.getProjection());

        map.setCenter(pos, 16);

        Ext.create('Ext.container.Viewport', {
        layout: 'border',
        rtl:true,
        bodyBorder: true,
        defaults: {
            collapsible: false,
            split: true,
            //margins: '0 0 0 0',
            bodyPadding: 15
        },
        items: [
        {
            //title: 'Header',
            region: 'north',
            height:50,
            html: 'Menu'
        },
        {
            title: 'Navigate',
            region:'east',
            collapsible: true,
            width: 250,
            html: 'Secondary content like navigation links could go here'
        },
        {
            region: "center",
            id: "mappanel",
            xtype: "gx_mappanel",
            map: map,
            stateful: true,
            stateId: 'mappanel',
        }
        ]
        });
    }
});


texts:
I have a raster (geotiff) and I would like simply to know some parcentile values, specifically the 80th, 85th, 90th and 95th, because I need to highlight all values above these thresholds in QGIS.

First I thought that in Propreties--> Style-->Cumulative count cut I had simply to write the percent value in Min or Max boxes (80 and 95 respectively, for example). But I discovered that if I write 80 in the first box or in the second one (next to the cumulative count cut), the value will be different.

So I tried to convert the raster in an Ascii grid, I opened the Ascii with text editor, I copied the cell values and I pasted in an Excel file. I deleted all the no data values (-9999 in my case) and I used the Percentile function in Excel to extract from the resulting matrix my values.
I thought that it would the right way to find the values, but I wanted to be sure, so I used the Max function in Excel just to see if the value was the same to the one I can find in Propreties-->Metadata-->Maximum value.... and the to values were different, so I don't trust in what I did in Excel.

I think that what I need is a very simple opration, but I can't figure how to do it. Any help? I hope it is all clear


texts:
Is it possible to open .ECW files with QGIS 3.2?

I have read Opening ECW files in QGIS?, but that applies to QGIS 2.14.3.

Is it possible in QGIS 3.2 without downloading anything as we cannot install anything at work or access the program files?


texts: Clip a large raster (TIFF) by a shapefile with a lot of polygons

texts:
I want to clip a large raster file (1.8 GB) by a shapefile with more than 100 smaller polygons. The polygons have the size of under 1% of the TIFF image.

I'll tried to use the clip function of QGIS, but when I start it, my whole system freezes in one second.

I tried it with only one smaller polygon and it worked. My workstation should manage the workload. I have 98 GB of RAM, a Xeon six core CPU and even a Titan GPU (if QGIS can utilize it).

In the best case I want to loop over every polygon and save every polygon as one file by it id name.


texts: ogr2ogr won't read directory of GeoJSON into GPKG

texts:
I have a directory of GeoJSON files, which I want to load as layers into a single new GeoPackage file. According to the GDAL GeoPackage documentation, I can name my destination GeoPackage file, enter a directory of shapefiles as the source, and it will load them each as their own layer/table. Indeed, I've tested this with shapefiles, and it works.

When I try the same thing with GeoJSON files, however, I get this:

Unable to open datasource ./seattle_washington_osm with the following drivers.
   …
-> GeoJSON
   …

Is this simply not possible with GeoJSON files, or am I doing something wrong?


texts: Use geometry binary predicate list to subset simple feature object?

texts:
I am trying to create a new simple feature from two overlapping features. I want to subset all polygons that are fully or partially located within my buffer (red border).

However, the st_overlaps generates  Sparse geometry binary predicate list. I have honestly no idea how to use this list to subset my original fc to create new geometry?. This suggestion using st_within does not work: nbrs.buff <- left[st_overlaps(buff, left) %>% lengths  > 1, ], as it subset all of the files, not ones overed by buffer. What I am missing?

Here is my dummy example:

# Load data
shp = system.file("shape/nc.shp", package="sf")

my.sf <- st_read(shp, quiet = TRUE)

# Convert crs to projected system to make buffer
my.sf.web<- st_transform(my.sf, 3857)

# Subset the data to create two independent shps
i = 10

# Split datasets in two files
one  = my.sf.web[i, ]
left = my.sf.web[-i,]

# Create buffer
buff = st_buffer(one, 40000 ) # distance

# CHeck which polygons overlaps with my buffer
out.overlap = st_overlaps(buff, left)

# Subset the overlapping polygons to create new geometry: DOES NOT WORK??
nbrs.buff <- left[st_overlaps(buff, left) %>% lengths  > 1, ]


texts: SQL Server query, simillar to Symmetrical Difference in ArcMap

texts:
i have two tables with geometry data in SQL Server. How can I get all geometry differences with spatial query? I tried:

select k.geom.STSymDifference(r.geom).ToString() from table1 k, table2 r

But it returns to many rows. I would like to get same result as I get it if I use Symmetrical Difference function in ArcMap.


texts: Using sf to combine polygons that share border with a polygon, but dropping the central polygon?

texts:
I have rows of respondents within census polygons and am trying to generate first-order aggregated polygons, excluding the census unit a respondent resides in. Put another way, I want to combine neighboring polygons but exclude the polygon that touches all its neighbors.
st_union has been helpful in combining a polygon with its neighbors, but my trouble now is with dropping the central polygon.
I hope to remain within the sf environment.


texts: Techniques for Filtering the Attribute Table in QGIS's Print Composer

texts:
For the legend of a trail map I'm working on I'd like to show the trails and their lengths. I have this data in the table. The problem is that not all of the features (lines) have a name or length. Some are ancillary trails that shouldn't be include in the attribute table in the Print Composer.

Are there any suggested techniques for filtering the list or creating some kind of "virtual" vector layer that will filter the data and I can point my Attribute Table to that table?

Sorting by name descending and limiting the number I'm showing is not really an acceptable solution.


texts: We put the location of the MXD at the bottom of every map so people can find it when looking at the final exported map (PDF).  Is there any way to make this a hyperlink so people can click on the link in the PDF and it will just open the MXD?


texts: How to convert CSV to point shapefile using Python?

texts:
import csv,sys
import osgeo.ogr, osgeo.osr #we will need some packages
from osgeo import ogr
from osgeo import *
import osgeo.osr as osr

def test():
    boolValue = True
    canvas = qgis.utils.iface.mapCanvas()
    spatialReference = osr.SpatialReference() #will create a spatial reference locally to tell the system what the reference will be
    spatialReference.ImportFromEPSG(4326) #here we define this reference to be the EPSG code
    driver = osgeo.ogr.GetDriverByName('ESRI Shapefile') # will select the driver for our shp-file creation.
    export_shp = r'E:\Folder'
    shapeData = driver.CreateDataSource(export_shp) #so there we will store our data
    outLayer = shapeData.CreateLayer(str('WaypointLayer'), spatialReference, osgeo.ogr.wkbPoint) #this will create a corresponding layer for our data with given spatial information.
    outDataSource = driver.Open(export_shp, 1)
    outLayer = outDataSource.GetLayer()
    layer_defn = outLayer.GetLayerDefn() # gets parameters of the current shapefile
    index = 0

# Reading the CSV file
    cnt = 0
    input_file = r'E:\Folder\CSV\points.csv'
    with open(input_file, 'rb') as csvfile:
        readerDict = csv.DictReader(csvfile)

# Create the fields in the layer
        for header in readerDict.fieldnames:
            field_name = ogr.FieldDefn(header, ogr.OFTString)
            field_name.SetWidth(50)
            outLayer.CreateField(field_name)

        for row in readerDict:
            print (row['latitude'], row['longitude'])
            point = osgeo.ogr.Geometry(osgeo.ogr.wkbPoint)
            point.AddPoint(float(row['longitude']), float(row['latitude'])) #we do have LATs and LONs as Strings, so we convert them
            feature = osgeo.ogr.Feature(layer_defn)
            feature.SetGeometry(point) #set the coordinates
            feature.SetFID(index)

            for field in readerDict.fieldnames:
                i = feature.GetFieldIndex(field)
                feature.SetField(i, row[field])
            outLayer.CreateFeature(feature)
            index += 1
    shapeData.Destroy()
    outDataSource.Destroy()
    conLayer = QgsVectorLayer(export_shp,'pointLayer',"ogr")
    QgsMapLayerRegistry.instance().addMapLayer(conLayer)
    canvas.refresh()

The problem is that initially it is not updating the data in the attribute table but it creates the features when I executed once. When I executed again it is updating the field data.  It sometimes create duplicates also.

What I am doing wrong here?


texts:
I have a large number of compressed GeoTIFF (600,000, Orthophotos) which I want to convert into an MBTiles (into a file).

I have summarized the TIFFs via GDAL in a VRT file but QGIS is not able to create an MBFile from it. GDAL_Translate somehow doesn't get any further either.

Has anyone ever processed similar amounts of data and maybe has a tip for me? I work with Windows if that is somehow important.


texts:
I have a geotiff file in UTM and units in meters. I used gdal_translate with this command :

gdal_translate -of GTiff -a_nodata 0 -a_srs EPSG:4326 old.tif new.tif

gdalinfo on new.tif shows that the file units are in degrees while the data (corners,origin ...etc) are given in meters which is wrong.
Is it possible to reverse/correct this action (on the new file) so I can get the file with either UTM in meters or go on and make the file in "real" degrees because I deleted the original file?

Thanks.


texts: What is ArcGIS's lightest-weight data type for web mapping apps to minimize network io?

texts:
I am wanting to move geometries/attributes to and from arcgis server and web client.  If, I've already thought through my network io vs app architecture and have concluded with architecture Z (really I just have enough time to optimize where easiest!..low hanging fruit!)

My specific use case is in regard to handling feature records varyingly based on attributes which would be boolean.  I will either put 1 or 0 in certain fields and do something useful for each of these?  If I am doing field creation via a dbf template on ArcGIS Server with python, what should my boolean fields be if I want to minimize traffic?  Text with 1 of Length 1 or Short of Precision 1 or other?

Please provide some explanation/link to help or other.  Thanks.


texts: Preserve coordinates when converting TIFF image to numpy ndarray

texts:
I have satellite image that I want to convert into numy array (and then to Pandas).
I already know how to do that thanks to Stefan, but th eproblem is that it does not preserve the coordinate data.
This is how I do it now:
import xarray
src=rasterio.open('14072020.tif')
array = src.read()

pd.DataFrame(array.reshape([13,-1]).T)

I have also tried to follow this answer (keeping the coordinate system of raster files in the resulting raster file after operation with numpy)
but I get all the time error

NameError: name 'gdal_array' is not defined

which doest allow me to open the image.
My end goal is to have Pandas table that contains  the bands values toegther with the coordinate values.


texts:
I have two polygon Layers, A and B correspondingly. The polygons of B are lying on top of layer A and also overlapping the polygons of A (I made an example of the attribute table for better understandings).

I want to keep track of how much percentage of polygon A has "value_1" in different categories in it ("value_1" is originally from polygon B).

Example. Let's say my attribute table looks like this:

A_ID  B_ID  area_a  area_b  %area  value_1
1      1     100     30       30     2
1      2     100     45       45     5
1      3     100     25       25     16

2      4     50      10       20     12
2      5     50      30       60     5
2      6     50      10       20     5

The Output should look like this:

ID_A  area_a  value_1 <10  value_1 <20
1       100     75%            25%
2       50      80%            20%

I guess I have to dissolve "A_ID" but I don't know how to calculate the field 'value1_<10' and 'value1_<20' while dissolving. Any suggestions?


texts: OpenLayers - TileStache Overlay: tiles loaded but not visible

texts:
I have been trying this since yesterday. If I load only the tilestache layer, it is visible. If I load both bing and tilestache layer, the tilestache layer is not visible. When I look at the firebug console, the tiles are being loaded.

_map = new OpenLayers.Map({div: 'map'});
bing = new OpenLayers.Layer.Bing({
    key: 'AqR9gcXm14-PKVQsWrL8H9A_RnusHZCzQKX_tbJegzU1EOZGafcjzCYbcpkwukQh',
    type: 'Aerial'
});

tilestache = new OpenLayers.Layer.XYZ(
        'hul',
        'http://www.hyderabadurbanlab.com/tiles/tiles.py/hul2/${z}/${x}/${y}.png',
        {
          tileSize: new OpenLayers.Size(256, 256),
          tileOrigin: get_point(17.4099, 78.4761),
          sphericalMercator: true,
          isBaseLayer: false,
          transparent: true,
          opacity: 0.6,
          transitionEffect: 'resizing',
          buffer: 1,
          wrapDateLine: true
        }
    );
_map.addLayers([bing, tilestache]);

Attached is the screen shot of how it looks.

When I load just the tilestache layer, it looks like this:

I want the Bing maps as the base layer and the tilestache layer on top of it as an overlay.

Thank you.


texts:
I'd like to calculate the total area of several polygons along a north-south axis. A bit like this http://www.closr.it/canvas/1307/#/spot

There is a similar question with R :
Spatial graphing, distribution of polygon area over latitudes . But I was wondering if there was not a simpler solution with QGIS, and if possible in a continuous way, not a discrete one.


texts: Handling KML/CSV with geopandas - DriverError: unsupported driver: u'CSV'

texts:
I'm trying to import a dataset of UK postcode districts UK postcode districts from Google Fusion tables into geopandas with the aim of splicing it to a dataset of house prices.

The problem is that geopandas doesn't seem to like either the CSV/KML version of the data.

If I use the command geopandas.from_file(path) then I receive the error:

DriverError: unsupported driver: u'CSV'

I have the same error when I try to import a KML version of the data (You can download the data below).

I'm sure it's trivial, but can anyone work out what I'm doing wrong?

UK postcode districts - KML

UK postcode districts - CSV


texts: Selecting features by radius using PyQGIS

texts:
in QGIS there's a very useful tool "Select features by radius", I'd like to add a similar function in my plugin. The user should be able to choose the radius and the center of the circle (from a chosen feature).

Is there a way to call "Select features by radius" using PyQGIS?

I can't find anything useful.


texts: Are there GUI resources for gdal and ogr?

texts: For those who are new to the programming side of GIS, I feel that a lot of great tools and resources are inacessible. I am interested in GUI based wizards for handling geoprocessing with libraries like gdal and ogr. Any tips?


texts: What are practical means for drawing an approximately rectangular geofenced area that includes the North Pole?

texts:
A naïve geofence approach is to take a "rectangle" that has a min/max latitude and a min/max longitude. This approach fails at the meridian, but that can be fixed by inverting some logic when "min" longitude is greater than "max" longitude. However, no such simple solution works at the North (or South) Pole. If you try a min/max latitude/longitude approach, you're going to get a triangle that folds back in on itself (and that's assuming the logic understands you want to cross the North Pole).
So, one solution that occurs to me is to have a flag in the algorithm for specifying that you want to include the North Pole, but that then requires you to use a completely different coordinate system, e.g., one with "poles" at the equator and +/- 180 degrees longitude. That's a very inelegant solution, but it does allow one to specify something approximating a rectangular region at the North Pole.
Surely someone has come up with a better system. What is it?
Note: assume that the rectangles are no larger than 1000 miles in either dimension, and will likely be much smaller.


texts: I am trying to clip the outside black boundary in between the TIFF files out.  The problem is the TIFF files being on a layer above the other TIFFs.  This is what causes the border and I want the picture to look consistent all the way through.


texts: 2 Shapefile vectors will not plot properly in same space

texts:
There are slight differences in 2 shapefiles projection information. I have activated on-the-fly, but they do not plot in the same space. I tried to save as new shapefiles with the same CRS, but nothing changed. Here are the projection files followed by the properties shown in Qgis.

  GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

Extents:

In layer spatial reference system units : xMin,yMin -176.732,51.712 : xMax,yMax
-130.017,69.5665

Layer Spatial Reference System:

+proj=longlat +datum=WGS84 +no_defs

  GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS
  84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]

Extents:

In layer spatial reference system units : xMin,yMin -608773.35,761075.16 : xMax,yMax
1463871.73,2053821.27

Layer Spatial Reference System:

+proj=longlat +datum=WGS84 +no_defs


texts:
I'm new in programming c# and and arcgis server . i just use arcgis server wizard to create services and web applications but now i'd like to learn more about arcgis server and programming with asp.net . I just found programming arcgis server with asp.net book. i didn't buy this book .i coudln't find any website about learning arcgis server .

Can you recommend other resources to learn arcgis server programming for beginners?


texts: How to fix layer styling in OpenLayers 5?

texts:
I'm trying to draw a new layer every time user clicks on the map.
My styling function can be seen below;

function CatchmentLayerStyleFunction(feature) {
    var colors = ['rgba(255,0,0,0.5)', 'rgba(0,255,0,0.5)',
        'rgba(0,0,255,0.5)', 'rgba(36, 201, 149,0.5)',
        'rgba(207, 242, 9,0.5)', 'rgba(241, 137, 9,0.5)',
        'rgba(8, 240, 240,0.5)', 'rgba(237, 118, 231,0.5)',
        'rgba(237, 26, 82, 0.5)', 'rgba(237, 26, 222)'];
    var color = new Style({
        stroke: new Stroke({
            color: colors[Math.floor(Math.random() * colors.length)],
            lineDash: [4],
            width: 3
        }),
        fill: new Fill({
            color: colors[Math.floor(Math.random() * colors.length)]
        })
    });

    return color
};

As you can see I'm randomly coloring the polygons.
My function that adds the layers to map is like ;

var vectorSource = new VectorSource({
                renderMode: 'image',
                features: (new GeoJSON()).readFeatures(data[0], {
                    dataProjection: 'EPSG:4326',
                    featureProjection: 'EPSG:3857'
                })
            });
            map.addLayer(new VectorLayer({
                'title': 'Catchment',
                source: vectorSource,
                style: CatchmentLayerStyleFunction
            }));
            map.getView().fit(vectorSource.getExtent());
        });

But as you can see in the gif, everytime I pan or zoom on the map styles of these newly added polygons are changing. How can one fix styles for the layers ?


texts: Exported QGIS model does not work: Invalid input error

texts:
This question is linked to this one:
QGIS Graphic modeler: Loading in all the layers
I have a very complex model which converts a specific file (pmkl.json) which is a feature collection of points lines and polygons. To its specifics based on theme ( a combination of Refactor and attribute selection algorithms). There are three inputs:
PMKL File = Which is where I browse to the PMKL.Json File to convert.
Target Location = Where the Shapefiles will be saved
CRS = CRS of the target shapefiles.
It is in the PKML File where the filtering of geometry happens, normally filtering by geometry works only with Vector layer. But I managed to get around it to address a file by using an expression.
But when I try to export the model as a script and run it: I get an invalid input.

Here is the code below that is relevant how do I solve this? Note as a model this thing works fine, but as Python script it does not. I want to make it a standalone application.
"""
Model exported as python.
Name : ODCM KLIP Omzetter
Group : ODCM KLIP Omzetter
With QGIS : 31401
"""

from qgis.core import QgsProcessing
from qgis.core import QgsProcessingAlgorithm
from qgis.core import QgsProcessingMultiStepFeedback
from qgis.core import QgsProcessingParameterFile
from qgis.core import QgsProcessingParameterCrs
from qgis.core import QgsProcessingParameterBoolean
from qgis.core import QgsExpression
import processing

class OdcmKlipOmzetter(QgsProcessingAlgorithm):

    def initAlgorithm(self, config=None):
        self.addParameter(QgsProcessingParameterFile('PKMLBestand', 'PMKL Bestand', behavior=QgsProcessingParameterFile.File, fileFilter='Alle bestanden (*.json*)', defaultValue=None))
        self.addParameter(QgsProcessingParameterFile('Mapvanuitgepaktezipbestand', 'Doel Locatie', behavior=QgsProcessingParameterFile.Folder, fileFilter='Alle bestanden (*.*)', defaultValue=None))
        self.addParameter(QgsProcessingParameterCrs('crs', 'CRS', defaultValue='EPSG:31370'))
        self.addParameter(QgsProcessingParameterBoolean('VERBOSE_LOG', 'Verbose logging', optional=True, defaultValue=False))

    def processAlgorithm(self, parameters, context, model_feedback):
        # Use a multi-step feedback, so that individual child algorithm progress reports are adjusted for the
        # overall progress through the model
        feedback = QgsProcessingMultiStepFeedback(71, model_feedback)
        results = {}
        outputs = {}

        # Filter by geometry type
        alg_params = {
            'INPUT': QgsExpression('@PKMLBestand ').evaluate(),
            'LINES': QgsProcessing.TEMPORARY_OUTPUT,
            'POINTS': QgsProcessing.TEMPORARY_OUTPUT,
            'POLYGONS': QgsProcessing.TEMPORARY_OUTPUT
        }
        outputs['FilterByGeometryType'] = processing.run('native:filterbygeometry', alg_params, context=context, feedback=feedback, is_child_algorithm=True)


texts: Exporting map from multiple .gpx files and georeference .tif

texts: I have a .pdf trail map of the Red River Gorge (RRG), that I got from fs.usda.gov (36x36 Map at bottom of the page) Have georeference map, saved as .tif, works great with Avenza! I would like to add things from previous trips from my .gpx files. I add my .gpx files via new vector layer, and it looks great with my tracks and waypoints on top of the base RRG map. But when I try to export or new print layout it is unreadable because you can no longer zoom in, the quality is gone, or just a small part of the map is visible. I have changed the layout size and dpi output. I do not want to change the quality of the base map just add my tracks, camping spots, and other features to the RRG map and keep the same quality. I'm sure that I can do this, just not sure how to go about it.


texts:
We are using the field mapping application QField with QGIS 2.18 for on-site data collection.

I've looked through the documentation and am still unsure, so am bring this question here: Is it possible to use python scripts using the python init function specification in the layer field properties?

I currently use a script to automatically assign object modification date and time. I am looking use this, because as far as I can tell, the now() function only works when the object is created. Updating existing points does not update the now() variable.


texts: Saving Edits in ArcMap automatically smooths line feature?

texts:
Using ArcMap 10.3.1: I'm having an issue with automatic smoothing of line features after editing vertices of line features and clicking 'Save Edits'.

I start my editing session, click on a line feature, click on a vertex, delete that vertex using the Edit Vertices toolbar, click 'Save Edits', and my line feature becomes very generalized/smooth instead of just removing that vertex. I'm having the same problem when adding a vertex.

Any thoughts on how to fix this?


texts: Passing a list of strings in a QGIS attribute table with column type: String List

texts:
I have a GeoJSON file in QGIS and I created a column in its attribute table with column type: String list.
I want to each cell to have two strings. The two strings are urls.
I tried:

'https://i.imgur.com/EKJmqui.jpg', 'http://freakplaces.com/img/Cannon_Beach17.jpg'

['https://i.imgur.com/EKJmqui.jpg', 'http://freakplaces.com/img/Cannon_Beach17.jpg']

list('https://i.imgur.com/EKJmqui.jpg', 'http://freakplaces.com/img/Cannon_Beach17.jpg')

These didn't work. Is there something I'm missing?
Relevant image: https://imgur.com/S07lmwb


texts: Problem with making raster from Aspect and wind direction data

texts:
I know its not correct to ask for work flow, but I am learning by doing and now I don't know how to move on.

I'm trying to make new raster which will represent the facing sides of hill around a factory. Then I'll use wind data (.csv) file to get which of these sides are more often under blowing wind. Wind data are for 10years, and the average is north-west flow.

I was looking for some index but I couldn't find one. So is there any chance to make it with the QGIS raster calculator?


texts: Using postgis to flatten all polygons in table

texts:
I have a table with many polygons. The polygons all have a name and a value.

I want to flatten these polygons so that each polygon border remains as border to the resulting polygons. Instead of having a stack of polygons I will remain with a mosaic of polygons based on the same lines.

Resulting in this:

The names are insignificant. I want to sum the values of each original polygon within the different tiles of the mosaic. The number of overlaps would also be relevant to calculate.

I am looking for a solution using PostGIS. I have solved this using FME, but my data size (20.000+ objects with millions of overlaps) makes this too time consuming. The union function in QGIS does this fairly straightforward as illustrated form their document webpages.

For those in need of a context the data is actually multi-ring buffers around seabird nesting sites in Norway. The data sets are per species and we are looking at several thousands nesting sites per species.

You may download a test data set here.


texts: Density map of pedestrian paths and footways in QGIS

texts:
I need to identify public places within pedestrian areas in London, so I need to calculate the density of both pedestrian links and footways. I have classified my data according to the categorization method so I can see where there is the highest concentration in the map. However I would need to calculate the density.

I have a line shapefile and I could extract the nodes and then create a heat map. However I am afraid considering just the nodes would alter the dataset. Any suggestion?


texts: ArcGis null values in interpolation tools

texts:
I'm making some maps in ArcGIS for some geochemical data. I have some null values in my attribute table but when I use the interpolation tools (kriging or inverse distance)  the null values are automaticaly considered as zero.

Is there a way I could make the interpolation tool to ignore the points where there is a null value insted of considering it as zero?


texts:
im working with Landsat 7 images and Geomatica Focus. I want to use the spectral plot to distribute the different wavelengths (spectral signature) of my image bands. In the tool(right click- > spectral plot) i only have the options for band numbers, not wavelength in nm.

How can i change this?


texts: Resizing objects from Cartesian measurements to Ellipsoidal

texts:
I have been working on a project in QGIS 3.12.2 in EPSG 3857 so it works with google earth. Without realising the difference there would be between cartesian and ellipsoidal I buffered some points.
By default Q does this using a Cartesian system. The rest of my map is brought straight into QGIS but measurements are ellipsoidal. This means my buffer that I specified to 2km actually is about 1.5km in respect to the rest of the map. Not realising this at first I proceeded to do a few days worth of work mapping paths within the 2k buffers which now I've realised are all too small in relation to the rest of the maps.
I was wondering if there is a way to rescale from cartesian to ellipsoidal measures so I can get everything back in line in one mass rescaling rather than having to redo everything.


texts: gdalbuildvrt error, when using in Python

texts:
I am new to python/GDAL and am running into perhaps a trivial issue. This may stem from the fact that I don't really understand how to use GDAL properly in python, or something careless, but even though I think I am following the help doc, I keep getting a syntax error when trying to use "gdalbuildvrt".

What I want to do is take several (amount varies for each set, call it N) geotagged 1-band binary rasters [all values are either 0 or 1] of different sizes (each raster in the set overlaps for the most part though), and "stack" them on top of each other so that they are aligned properly according to their coordinate information. I want this "stack" simply so I can sum the values and produce a 'total' tiff that has an extent to match the exclusive extent (meaning not just the overlap region) of all the original rasters. The resulting tiff would have values ranging from 0 to N, to represent the total number of "hits" the pixel in that location received over the course of the N rasters.

I was led to gdalbuildvrt [ http://www.gdal.org/gdalbuildvrt.html ] and after reading about it, it seemed that by using the keyword -separate, I would be able to achieve what I need. However, each time I try to run my program, I get a syntax error. The following shows two of the several different ways I tried calling gdalbuildvrt:

gdalbuildvrt -separate -input_file_list stack.vrt inputlist.txt
gdalbuildvrt -separate stack.vrt inclassfiles

Where inputlist.txt is a text file with a path to the tif on every line, just like the help doc specifies. And inclassfiles is a python list of the pathnames. Every single time, no matter which way I call it, I get a syntax error on the first word after the keywords (i.e. 'inputlist' in inputlist.txt, or 'stack' in stack.vrt).

Could someone please shed some light on what I might be doing wrong? Alternatively, does anyone know how else I could use python to get what I need?

Thanks so much.


texts: Plotting Multiple Points on Google Map using VB.net

texts:
I am creating a database application using VB.net and MSSQL 2012. I have retrieved a list of customer from the database and want to plot them on google map. At it stands, I can only plot the first customer. Below is my code. Can someone help me out as this has been killing me for some times now.

    Dim lat As String = String.Empty
    Dim lon As String = String.Empty

    Dim queryAddress As New StringBuilder()
    queryAddress.Append("http://maps.google.com/maps?q=")

    For i = 0 To Me.DataGridView1.RowCount - 1
        If Me.DataGridView1.Rows(i).Cells(0).Value <> String.Empty Then
            lat = Me.DataGridView1.Rows(i).Cells(0).Value
            queryAddress.Append(lat + "%2C")
        End If

        If Me.DataGridView1.Rows(i).Cells(1).Value <> String.Empty Then
            lon = Me.DataGridView1.Rows(i).Cells(1).Value
            queryAddress.Append(lon)
        End If
        WebBrowser1.Navigate(queryAddress.ToString())
    Next