texts:
I am trying to change the visibility range for several layers on an ArcGIS web map - should be simple right??

This could be a case of my misunderstanding on how to use a tool or it could be a case of a bug.

I have tried to adjust the Set Visibility Range tool by moving the white arrows.  There is no save or OK button.  When I have clicked outside of the tool box and then back into the tool I see that it has not saved the scales I want.

I would love to simply be able to write a minimum scale of 1:3000 and no maximum.

Does anyone have any suggestions how I can work around this?

The layer is a web service which is used within many of our maps so I cannot simply alter the visibility range before publishing the service. This is the Service that I have been using.

I had tried using the scale range under the change symbols option.  This appeared to work the first time but the settings did not save properly and later in the afternoon it was back to default.

I am using Chrome.

19/04/2016:  I have tried using the change symbol options using Internet Explorer 11 (used Chrome yesterday) but again it does not appear to save the scale I want.

28/04/2016:  The original service I linked to was the wrong service so I have corrected that link.
I have trialed the service on a test map I have just created and this has worked fine.  The problem appears to be in my original online map.  As you can see in the screenshot above the scale bar goes off the chart.

15/09/2016:  I am still having issues with this same map.  I have added another layer that I wanted a limited visibility range for. Same errors happening, again this works fine in a blank map.  It appears that arcgis webmaps are easily corruptible.


texts:
I have calculated convergence index on my DEM with SAGA.
Output results seems to be ok (values in range -100;100), but I don't understand how it works.

Here is some sample data from my DEM:

I know that Convergence Index is calculated based on aspect, therefore I've calculated Aspect in SAGA (to understand whole process), however it didn't help me. I still don't know how CI was calculated.

I've found this formula in some paper:

"Where θ denotes average angle between the aspect of adjacent cells and the direction to the central cell."

Can someone explain to me step by step how CI was calculated (for example on values inside red frame)


texts: Is it possible to 'clip' two line layers in QGIS?

texts: Specifically, I am looking to merge my bus route network into my wider road network. I only want to work with roads which lie on the bus network. Is this type of clipping possible? Joining the attributes of the layers is proving problematic.


texts:
I have an query that produces an array of coordinates like so:
95.61,38.60 95.22,37.98 95.60,37.66 94.97,37.65  (long/lat separated by "," and coordinates separated by a space).

How can I use these arrays to dynamically draw polygons on my map?

Ideas?

edit 1

Long answer…  I have a table with rows like so:
Area 1      lot1        95.61,38.60 95.22,37.98 95.60,37.66 94.97,37.65
Area 1      lot2        95.63,38.65 95.27,37.98 94.60,39.66 92.97,37.64
Area 2      lot 3       95.64,38.63 95.24,37.95 95.66,37.62 94.94,37.61
The first column has a checkbox.  As I select a checkbox, I want to draw a polygon using the coordinates from each row.


texts: Co-registering multiple temporal orthomosaics?

texts:
I have multiple (19) orthomosaics of the same approximate region taken regularly over the course of a couple months with a UAV.  Some, but not all, include ground control points we placed in the scene.  In all cases, however, there are features we can pick out that will be in the same position in each orthomosaic.  I need to co-register these images such that they are closely spatially matched.  They're currently within about 5 meters of each other, but I need to get that down to centimeters.

How do I go about doing that?

ENVI has the SPEAR change-detection tool that works on two images at a time.  I could step through all the images, one pair at a time, using this tool.

What other, potentially better, options are out there?

Here are two example images:


texts: OpenLayers 3 don't go beyond given coordinates

texts:
In Openlayers 3 Max Extent there is mentioned how to set a max extent. But how can I set coordinates in such a way that it is not possible to drag the map beyond the some given coordinates?

Extent is used to make sure that at least one corner of the rectangle is still in visible, but I want that you can't even drag beyond the rectangles borders.


texts: using polygon batch overlay of lecos error

texts:
I´m trying to analyse corine land cover with polygon batch overlay of lecos and do get this error:

Ein Fehler trat bei der Ausführung von Python-Code auf:

AttributeError: 'NoneType' object has no attribute 'GetFeatureCount'
Traceback (most recent call last):
  File "/Users/deniseehrhardt/.qgis2/python/plugins/LecoS/lecos_dlg.py", line 876, in accept
    self.Output(results)
  File "/Users/deniseehrhardt/.qgis2/python/plugins/LecoS/lecos_dlg.py", line 951, in Output
    attr = func.getAttributeList(self.vector,self.LandID)
  File "/Users/deniseehrhardt/.qgis2/python/plugins/LecoS/lecos_functions.py", line 248, in getAttributeList
    for i in range(0,d.GetFeatureCount()):
AttributeError: 'NoneType' object has no attribute 'GetFeatureCount'

Python-Version: 2.7.10 (default, Oct 23 2015, 19:19:21) [GCC 4.2.1 Compatible Apple LLVM 7.0.0 (clang-700.0.59.5)]
QGIS-Version: 2.18.9 Las Palmas, exported

Here´s a screenshot... after pressing ok the error occurs...

Can anyone help?


texts: Wherer is wrong with transforming a WGS84 [lon,lat] point onto google maps

texts:
I want to convert a longitude latitude point in WGS84 coordinate system to Google Maps(not google earth) coordinate, But when I put the converted point on google map, it always has a roughly 10km offset(moved to south) on y axis.

First correct me if my understanding is wrong with the two projections:

Google Maps:
    spheroid parameters: a = b = 6378137
    projection: Mercator;
    To WGS84 Datum offset: no, (means 0,0,0)

WGS84:
    spheroid parameters: a = 6378137;b=6356752.314
    projection: WGS84;(like no projection, the coordinate is just longitude and latitude)

What I did:

Convert the WGS84 longitude, latitude(we call it point1) to geocentric coordinate(point2) with WGS84's spheroid;
//because no datum offset between Google Maps and WGS84, no need to convert point2 to Google Maps' geocentric
Convert point2 to geodentic coordinate(point3) with Google map's spheroid
Project point3 with Mercator projection into meter coordinate(point4).
Linearly transform point4 onto google maps.

On internet everybody is saying that Google Maps' projection is 'EPSG:3857' which means a and b are the same 6378137, but if I just use standard WGS84's spheroid parameters, which means a and b are different, the whole process is always correct, the result is very accurate.


texts: Solution for android map app with offline use?

texts:
I am in search of a way to develop an app for android with embedded maps which can work without an internet connection(i.e. can i get the tiles of the open streetmap for offline use?).

I have developed web apps before by using the google maps api v3.
I need a general answer as a specific one.I'm trying to find the concepts and technologies i need to learn :)

(Meanwhile i found osmdroid.Does it work offline?)


texts: Raster product in PyQGIS does not produce correct results

texts:
I'm trying to mask (set with no data values) pixels in a raster not fulfilling a condition based on another raster. It runs well but it produces huge numbers that are not the values in the raster that I want to mask.
If someone can try the code, can figure what is wrong with it?
import processing

lyr1 = QgsRasterLayer('C:/Users/LAPTOPJOSE/Desktop/prueba_normalizar_raster/negras.tif')
conras = QgsRasterLayer('C:/Users/LAPTOPJOSE/Desktop/prueba_normalizar_raster/Insar_13_25_02_2021_coh.tif')
output = 'C:/Users/LAPTOPJOSE/Desktop/prueba_normalizar_raster/negras_null.tif'
entries = []

con=QgsRasterCalculatorEntry()
con.ref='con@1'
con.raster=conras
entries.append(con)

ras=QgsRasterCalculatorEntry()
ras.ref='ras@1'
ras.raster=lyr1
entries.append(ras)

#Set_null
calc=QgsRasterCalculator('"ras@1"/("con@1" > 0.65)',output,'GTiff',\
lyr1.extent(), lyr1.width(), lyr1.height(), entries)
calc.processCalculation()

It seems an error of the multiplying operator because it yields in very long figures even when multiplying two 32bits float rasters with the Raster Calculator tool. I attach a screenshot of the inputs ("Insar....coh" as con@1, "negras" as ras@1) and the resulting raster ("negras_null"):


texts:
I have a table of linestring segments in Postgis. These segments make up distinct non-intersecting lines as shown by the diagram. To add confusion the direction of the linestrings are not uniform.

I want to create a simple set of linestrings that are uniformly joined in a single direction, this will enable me to shorten them if required using ST_LineInterpolate().

So, looking at the diagram you can see three lines which is what I want to get out of the query. I tried a recursive query but this ran forever on a table of only 10 linestrings:-

WITH RECURSIVE line_join(ogc_fid,wkb_geometry) AS (

 SELECT ogc_fid, wkb_geometry FROM TABLE.tll_lines_selected WHERE ogc_fid=47
UNION ALL
     SELECT n.ogc_fid, n.wkb_geometry FROM TABLE.tll_lines_selected n, line_join lj
     WHERE ST_Intersects(n.wkb_geometry,lj.wkb_geometry)

) SELECT ogc_fid FROM line_join;

Therefore I think something is not right. Is there a better strategy?


texts: Generating multiple page reports (Page Per Report) using ArcGIS Pro & Python

texts:
I have reached a point a little bit outside of my knowledge frame. Python.
I found that I can create multiple page reports with Python in ArcGIS Pro by combine 2 reports together. I am following this link from ESRI's site. I am using Report example 3.
I continue to get the following message :
Traceback (most recent call last):|File "<string>", line 6, in <module>
AttributeError: 'str' object has no attribute 'aprx'

I have tried several thing but lack an understanding.
What would my next step be?
Create a tool box that would allow me to make this a little more stream lined?
How do I define 'self'?


texts: Retrieve single band from multiband raster served by Geoserver

texts:
I have a Geoserver imagemosaic layer serving out a timeseries of 2band rasters with band names GRAY_BAND and Band2. When I query the layer both bands are returned I believe the map only shows Band2, but if I click on a point in geoservers openlayers previewer it shows values for both bands. Is there a way specify via a param which band is returned? My goal is to get the single band so I can apply an SLD color ramp to it.

I did see that you can use propertyName=GREY_BAND, which specifies the band in the GetFeatureInfo response, but even when specifying propertyName it seems like the SLD is applied to a composite of both bands.


texts: Reprojecting a UTM map to equidistant cylindrical projection

texts:
I just recently started to work with GRASS GIS (version 6.4.4)... as I imported a UTM map of Svalbard, I would like to change its projection to equidistant cylindrical. After looking up the EPSG code for equidistant cylindrical I tried changing the projection information with g.prop epsg=1029 (or 9823), but all I got was an error message:

ERROR 6: EPSG PCS/GCS code 1029 not found in EPSG support files.

Is this a valid EPSG coordinate system?

How can I apply this projection to my map?

Do I have to modify the EPSG support files?

This does not work:

yadgar@HEXACORE:~$ grass -c EPSG:54002 /home/yadgar/grassdata/svalbard/svalbard/EquiCyl/PERMANENT
Creating new GRASS GIS location/mapset...
FEHLER: ERROR 6: EPSG PCS/GCS code 54002 not found in EPSG support files.  Is this a valid
EPSG coordinate system?


texts: Clip_analysis with several in_features and several clip_features

texts:
I am trying to Clip several in_features with several clip_features and name the output using the name of the clip features.

Till now i made it with a loop for many in_features but only one clip_feature. I used a for loop. But now i have no idea, how to that with several clip_features...

That's my code:

import arcpy
import os
arcpy.env.workspace = "path of in_features"
out_clipped_Objects = "path of output"
clip_fc = "path of clip_features"

fcs = arcpy.ListFeatureClasses()

for fc in fcs:
     arcpy.Clip_analysis(fc, clip_fc, os.path.join(out_clipped_Objects, fc)+ "name of clip_feature")

Does somebody know how to alter the code, to do the clip_analysis with more then one clip_feature?


texts: QGIS small install package for use in low-bandwidth settings?

texts:
I'm a cartographer and data analyst with experience using ArcGIS and QGIS to format data and produce maps for a global public health contract.  I am working with a group of QGIS learners from Africa and India who would like to be able to map health data. I've written instructions on some basic QGIS functions and have been helping them online.  But several of them have had trouble just downloading the current installer for QGIS, which is around 300Mb in size, prohibitive to those with slow/unreliable internet connections.

I was reading another post about other ways to do more advanced installations of QGIS:
Is there a Small Windows install of QGIS?

...and it appears a better way to do a smaller download/install might be through the use of the OSGeo4 installer, which lets you install only the pieces of the software you need.  My trouble is, I do not know which pieces are crucial.  (I need for the group to be able to handle really only vector data files.  The most common data types they will be working with at this point include CSV files and shapefiles. We have not gotten into any raster data analysis.)

Can anyone make any suggestions?  Which packages, at a minimum, will they need?


texts: Gray tiles when using OpenTopoMap on R-leaflet

texts:
I'm trying to use OpenTopoMap on a R Leaflet map. The basemaps from OSM, Esri and Stamen loads fine (thus no proxy problem) but when I try to load OpenTopoMap, all I get are gray tiles.

Anyone else has this problem? Bellow is the R code I use:

require(leaflet)

m <- leaflet() %>%
  addTiles(group = 'OSM') %>%
  addProviderTiles(providers$OpenTopoMap, group='Topo') %>%
  addProviderTiles(providers$Esri.WorldTopoMap, group='Esri') %>%
  addProviderTiles(providers$Stamen.Toner, group = "Toner") %>%
  addLayersControl(
    baseGroups = c('OSM', 'Topo', 'Toner', 'Esri'))

m


texts:
How can I see in Carto the density of data points geolocalized if some of them are located in the same coordinates? In an example:

If I have 5 coordinates as follows:

50.64815,12.96872
48.65920,8.94052
48.65920,8.94052
48.65920,8.94052
48.65920,8.94052

In Carto I will be able to see 2 datapoints (because four coordinates are the same). However, the dimension of these 2 datapoints are displayed with exactly the same dimension.

I'd expect to see a map showing these 2 datapoints but one 4 times bigger than the other.

Any advice?


texts: Why does a .gdb I created with GDAL double in size after it's been modified in ArcMap?

texts:
I created a large (over 200 million features) .gdb file with Python/GDAL. I had to move it to my Windows machine to add and calculate a column in Arc, since I couldn't figure out how to do with GDAL. The only difference was the addition of the one long integer column, but when I move it back to my Linux VM, the .gdb is now double the size (244G vs 112G).
Why is that and is there anything I can do to change it?
Relatedly, a .gdb I had created earlier using the same data and process (but for a different part of the world) had almost 5 times fewer features (44.6 mil), but is 12.2 times smaller than new .gdb (20G vs. 244G). Granted, there are 9 fewer columns, but this still doesn't make sense to me. Another interesting thing is that this .gdb, when I created and calculated the column in Arc a few months ago, had SHAPE_Length and SHAPE_Area added as columns, which I assumed was an Arc thing. This time, however, those fields were not added.
Anyone has some insight as to why I may be seeing this behavior? I am just trying to understand.


texts:
I am trying to utilize GeoExplorer (made by OpenGeo) to browse data provided by GeoServer.

So far, I have downloaded source code from https://github.com/opengeo/GeoExplorer. Some instructions I found mention Ant and Git software to build app., but I can't make it work.

Could somebody suggest me how to set up GeoExplorer (non-Suite) on Windows? If you could ouline the required steps, that  would be great.
Thanks,


texts: What extra CS courses should a GIS Developer take?

texts:
I'm a CS major who is minoring in GIS to become a GIS Developer. I was wondering what later CS (or Software Engineering and Computer Engineering) courses would be useful to a GIS Developer?

For the major, I have to take CS1, CS2, Design Patterns (Software Enginering 101), Computer Organization (Comp Engineering 101), Data Structures/Algorithms, Operating Systems, and Database Management.

What more classes will be most useful in CS to a GIS Dev, 3 of the following?

I have options of Web Programming, Intro to Computer Graphics, Adv Computer Graphics, Computer Graphic Algorithms, Compiler Design, Networks, more Software Engineering courses, or more Computer Engineering?

My advisor doesn't know much about GIS, and I haven't taken any GIS courses yet so I don't know any GIS professors (next semester I will start).

I just need help future planning because if I do the wrong sequence, I could screw myself over prereq wise since some of these have weird offerings like fall of even years or spring of odd years etc.


texts:
Similar to Reclassify raster values of 9999 to NoData, leaving all other values alone, in ArcGIS Desktop?, I am attempting to convert numerical NoData values to actual NoData cells. The solution in the linked question doesn't work for me, however, as when I use the save method to save to TIFF the numerical values are reintroduced.

When I do this manually (using the GUI tools) I do not have the issue. When I use  arcpy.gp.SetNull_sa (over arcpy.sa) I do not have the issue, but using the gp approach doesn't work on one of our team member's machines for some reason. I have not found much documentation on the gp appr, either.

I can save the outputs as GRID and export using the context menu because that tool allows you to supply a NoData value, but the Raster To Other Format (Conversion) tool doesn't give that option. It appears to work, regardless, but does not let you name the output (only the workspace), which means there's another step to add to the script (renaming). This solution seems too convoluted to be the right one, and also takes twice as long as the other methods.

So, what are methods for saving the output from an arcpy.sa tool in TIFF format with NoData intact?

Environments:

arcpy.env.compression = "NONE"
arcpy.env.rasterStatistics = "STATISTICS 1 1"
arcpy.env.resamplingMethod = "NEAREST"
arcpy.env.snapRaster = inDEM
arcpy.env.tileSize = "128 128"
arcpy.env.pyramid = "PYRAMIDS -1 NEAREST DEFAULT 75 NO_SKIP"
arcpy.env.extent = inDEM
arcpy.env.nodata = "NONE"

This does not work:

con_ras = arcpy.sa.Con(arcpy.sa.IsNull(outBreachrast),outBreachrast,inDEM,""""Value" = 0""")
con_ras.save(tempDEM)

null_ras = arcpy.sa.SetNull(arcpy.sa.IsNull(inDEM),tempDEM,""""Value" = 1""")
null_ras.save(outDEM)

This does (but not on everyone's computers):

arcpy.gp.IsNull_sa(inDEM, isnullDEM)
arcpy.gp.SetNull_sa(isnullDEM, tempDEM, outDEM, """"Value" =1""")

This also does, but does not let you specify file names:

con_ras = arcpy.sa.Con(arcpy.sa.IsNull(outBreachrast),outBreachrast,inDEM,""""Value" = 0""")
con_ras.save(tempDEM)

null_ras = arcpy.sa.SetNull(arcpy.sa.IsNull(inDEM),arcpy.sa.Con(arcpy.sa.IsNull(outBreachrast),outBreachrast,inDEM,""""Value" = 0"""),""""Value" = 1""")

tiff_ras= arcpy.RasterToOtherFormat_conversion(null_ras, outDir, "TIFF")

This behaviour may only happen in ArcMap, and may not be an issue in ArcCatalog, which is puzzling.


texts: What can I do (print, publish…) with Esri's "free" maps accessible via their REST API?

texts:
I am mostly working with Python, and the matplotlib basemap toolkit has so far been sufficient for my mapping and plotting needs. However now I managed to get to an image size, where the included map backgrounds look rather bad.
Searching the web for possible high resolution backgrounds ESRI_Imagery_World_2D and similar things show up, e.g. here or here.
But I can't find a license for this data that clearly applies.
Can I use the above methods and outputs for a scientific publication?
The API website does not give me any license info. (Sorry, out of links…)
The "HUMAN READABLE SUMMARY OF THE TERMS OF USE " for ArcGISOnline seems related, but I think it applies to stuff in their web viewer, and not via their API. However it does say I may:

Include screen captures or a printed or plotted maps in the following ways

…

In academic publications (research journals, textbooks, etc)

Which would apply to my situation, but I am not sure if it also applies to stuff accessed via the REST API.


texts: Maptip containing an expression combining a CASE WHEN-clause and a hyperlink consisting partly of value from table column

texts:
I am using maptips in QGIS 3.4, but the below expression won't work.

What I am trying to do is to use a CASE WHEN-clause in combination with a HREF from a webserver where the last part of the URL comes from the table column.

Here is my code:

 [%CASE WHEN "pdf_link3" IS NULL THEN 'No other documents' ELSE 'Click here'
END %]

'pdf_link3' is a table column and contains the name of the pdf document, for instance 'lokalplan.pdf'.

Putting it like the above QGIS the CASE WHEN works, but the URL is loaded without the last part expected to be loaded from the attribute table.

Can anybody figure out how to make this work?


texts: ExtJS custom popup on OpenLayers 3 (and the terrifying event propagation problem :O)

texts:
I'm trying to make work a custom popup with a datagrid. The popup visually displays correctly, but the problem is how events propagate.

The popup is an ext panel (with a datagrid) who is rendered on a div. This div is set as the element of the ol.Overlay added to the map.

The event propagation problem

For the event propagation, the ol.Overlay have the 'stopEvent' option. If it is setted to true, the click events doesn't work on the popup. You can select elements on the grid, or select the 'Save' button but to press the button (or edit the grid element) you have to select the element and then press enter on the keyboard.

If 'stopEvent' option is setted to false, the popup events fires ok, but if you double-click the popup the map is zoomed, if you drag the popup the map is dragged too, etc...

You can view and test the code here:

http://jsfiddle.net/j8fLmju8/

Sorry I couldn't make the styles work on the fiddle it should look like this:

Another trouble with the fiddle (that works ok on my app) is when stopEvent is setted to true, if you select an element on the grid and press enter it always edit the first value, just ignore that.

Finally

I don't know what should be the next step. Try to use stopEvent=false and do something to priorize the panel events over the ol.Overlay? Use stopEvent=true and try to stop the map events when clicked on the panel? i don't even know if I can do anithing of this. Another option?

PD: I've tried the GeoExt 3 popup but it doesn't fit to my needs.


texts: Calculating distance of polygons from coastline using ArcGIS Desktop?

texts:
I have a series of polygons representing parcels of land. I would like to calculate the distance of each polygon from the coastline. I also have a polygon of the entire region of interest, i.e. an island.

What is the ArcGIS Desktop 10.1 procedure to achieve what I am after?


texts: How to use Morton Order(z order curve) in range search?

texts:
How to use Morton Order in range search?

From the wiki,

In the paragraph "Use with one-dimensional data structures for range searching",
it says:

  the range being queried (x = 2, ..., 3, y = 2, ..., 6) is indicated by
  the dotted rectangle. Its highest Z-value (MAX) is 45. In this
  example, the value F = 19 is encountered when searching a data
  structure in increasing Z-value direction. ......BIGMIN (36 in the
  example).....only search in the interval between BIGMIN and MAX....

My questions are:
1) why the F is 19? Why the F should not be 16?

2) How to get the BIGMIN?

3) Are there any web blogs demonstrate how to do the range search?


texts: Income data on a block level: how does a newbie do it?

texts:
I am a complete newbie at GIS data, and am trying to do the following:

Figure out census blocks in Miami-Dade (and several other counties) that have average household incomes > $75,000
Get the latitude-longitude boundaries of those blocks
Plot the blocks on a map

1 and 2 are the  most important. What is the best way to get this data, and then the lat-long boundaries? Quickest, easiest, cheapest are of course the top criteria.


texts: ArcGIS JS API 4.11, 4.12 Incorrect work GraphicsView2D (from esri/views/2d/layers/grapics/GraphicsView2D) hitTest method

texts:
ArcGIS JS API 4.11, 4.12
GraphicsView2D (from esri/views/2d/layers/grapics/GraphicsView2D) hitTest method return only one feature (graphic) with zero index.
Better to be return all graphics.

Is it issue or correct work?

UPD:
I opened source code https://github.com/Esri/jsapi-resources
Go to 4master branch
Open /views/2d/layers/graphics/GraphicsView2D.js

in hitTest method see implementation like

t.prototype.hitTest = function (e, t) {
        if (!this.view || !this.view.position) return l.resolve();
        var r = this.view.toMap(g.createScreenPoint(e, t));
        return this.searchFeatures(r).then(function (e) {
            return e && e.length ? e[0] : null
        })
    }

From all features from graphicsView return only one value (e[0] with zero index)

I created example https://jsfiddle.net/catcherholms/dmp1z387/143/ with many objects

when I clicked here in hitTest result I receive one feature and should be more than one


texts:
In pure Mapinfo SQL Update code I need to extract the last word of a string.

Ex

SUTHERLAND STREET
BURKE DEVELOPMENTAL ROAD
Street Not Specified

The following Update DCDB_DNRME Set road_type = Mid$(road_name,(InStr(1,road_name," ")+1),len(road_name)) gives me

STREET
DEVELOPMENTAL ROAD
Not Specified

I just need

STREET
ROAD
Not Specified

Is there an equivalent to rfind in python?


texts: How to get 'GetPropertyValue' result in JSON format in geoserver?

texts:
I have tried adding outputFormat = application/json but not working :( my URL looks like this
http://domain/geoserver/wfs?service=wfs&version=1.3.0&request=GetPropertyValue&typeNames=forest:division&valueReference=name&outputFormat=application/json

I also tried WFS version 2.0.0
http://domain/geoserver/wfs?service=wfs&version=2.0.0&request=GetPropertyValue&typeNames=forest:division&valueReference=name&format=application/json

I'm getting results as follows. :


texts: Finding points a certain distance from each other

texts: I would like to find the maximum amount of points within a specified domain that are a certain distance from each other. For example, how many points within an area are 1000ft from each other? I envision each point with a radius of the desired distance, but I don't know if I must start with a point or what.


texts:
I am trying to show all world borders. Either using OpenLayers or Leaflet. World borders are taken from here. This code does not work.

What I am doing wrong?

var Map = {

  map: null,
  layer: null,

  initializeMap: function() {
    map = new OpenLayers.Map('map');
    layer = new OpenLayers.Layer.OSM("World Map");
    map.addLayer(layer);
    map.zoomTo(2);
    var geojson_format = new OpenLayers.Format.GeoJSON();
    var vector_layer = new OpenLayers.Layer.Vector();
    map.addLayer(vector_layer);
    vector_layer.addFeatures(geojson_format.read(featurecollection));
   }

}


texts:
I'm new when it comes to GIS and PostgreSQL. I have searched everywhere for the cause of this problem but I can't find the solution. I was following the documentation here

First i created a table:

CREATE TABLE tableA (id integer, polygon geometry);

I couldn't use st_polygon because that was causing a "type not exist" error though GIS is fully updated.  (through Ubuntu repos)

Then I inserted a row. I couldn't do it like how it was described in the documentation because that was giving me a weird errors:

INSERT INTO tableA VALUES (1,
    st_polygon ('polygon ((40 40, 40 60, 60 60, 60 40, 40 40))', 0)
    );

ERROR: lwline_deserialize: attempt to deserialize a line which is really a Invalid type
    SQL state: XX000
    Context: SQL function "st_polygon" statement 1

INSERT INTO tableA VALUES (1,
    st_polygon ('polygon ((40 40, 40 60, 60 60, 60 40, 40 40))', 0)
    );

  ERROR: function polygon(unknown, integer) does not exist
SQL state: 42883
Hint: No function matches the given name and argument types. You might need to add explicit type casts.
Character: 44

In the end I was able to insert it like this:

INSERT INTO tableA VALUES (
1,
ST_PolyFromText ('polygon ((-30.8947 138.6522,
-34.8951 138.6522,
-34.8952 138.6512,
-34.8947 138.6512,
-30.8947 138.6522))', 1)
);

Now when I perform

SELECT * FROM tableA WHERE ST_OVERLAPS(ST_BUFFER(GeomFromText('POLYGON((-34.8947 138.6522,
-34.8951 138.6522,
-34.8952 138.6512,
-34.8947 138.6512,
-34.8947 138.6522 ))'),0), polygon);

It outputs the error:

ERROR: Operation on two geometries with different SRIDs
SQL state: XX000

Anybody got a clue what it could be now and how to fix it?


texts: How to Combine Two Image Overlaying Into One Image in Google Earth Engine?

texts: I have two images of different sizes. I can make them overlay each other using Map.addLayer() for each. But I want to combine them into one image as shown in illustration below. How to achieve this?


texts: Is it possible to allow ArcGIS Online users to set visibility range?

texts:
I'm not certain you are able to do this or not but is it possible to create a widget or something in an AGOL app to allow a user to set their own visibility range for a layer?

This is outside of me setting the visibility range for everyone as an administrator. I've had several requests from my users to allow them to change their own visibility ranges for certain layers.


texts: Publishing a big raster layer via GeoServer

texts:
I have a big raster layer (7GB) and I want to publish it so it can be accessible from my Geoportal using OpenLayers 4. I couldn't find a clear procedure to do that, I started by tiling the image but I don't know next steps.

gdal_retile.py -v -r bilinear -levels 4 -ps 2048 2048 -co "TILED=YES" -co "BLOCKXSIZE=256" -co "BLOCKYSIZE=256" -targetDir C:\Users\Downloads\Compressed\bmpyramid C:\Users\Downloads\Compressed\img.tif


texts:
I'm working in a GIS program with C++ and in this moment, i'm working in my raster visor.

I need reproject my Raster in ejecution time, but i dont want create a output file with this reprojected Raster. For this reason i can't use GDALWarp.

Does anyone have any idea how to reproject the raster without creating a new output file?

I think if I got the GeoTransform information from the reprojected Raster it would be enough


texts: How to split one modelbuilder process in several GUI windows?

texts:
I am working with ArcGis 10.2.

I would like to know how I can have the different steps of a modelbuilder's process displayed in successive windows once the previous main step is done? A bit like in a AV script.

Let's say, I have a process with 4 main steps (each has at least 6 to 10 Tools/steps) requiring each time different Inputs, Parameters and Outputs. What I will get as "GUI" tool will be a very long list of parameters to inform.

I would like to avoid this in order to have an easier tool to handle by "basic" ArcGis users.

I am quite Ok with the modelbuilder but sadly don't have skills in coding.

Many thanks

I have edited my post with a screenshot. I'd like to have the possibility to ask the user if the first Step (in the middle) needs to be done or not
- regarding the answer either run the second step (last box), after execution of the first step  or directly after the question " do you need to reclass your data?"

I hope this is helpful.
Many thanks


texts: Join attributes by location in GRASS GIS

texts:
I have spatially superimposed points in two layers from which I want to copy one column of attributes to the other layer. In QGIS you can do this with 'Join attributes by location' which worked great.

Now I want to write these steps in an R script using the R package rgrass7, which allows me to include GRASS GIS in R. Unfortunately I can't find a suitable equivalent to 'Join attributes by location' in GRASS GIS.

To give more background information:
I want to calculate the distance (as the crow flies) between locations and a river system. Therefore I calculated a network with v.net.connect and extracted the sections between the locations and the river with v.select.intersect. From this I got lines (of the sections), from which I could calculate the length with v.to.db, which I then saved in a column in the sections. Then I used v.to.points to create end nodes which got the length information. And now I want to copy exactly this length information into the attribute table of the locations. The places contain important information like a column with numbers which I need for further calculations and without I can't assign these points to other data later.

This may be a lot more complicated than it has to be, or 'thought around the corner', but in QGIS this method has worked.

A selected section between the location and the river with the length in the attribute table I want to transfer to sites:


texts: Matching ArcGIS and QGIS coordinate systems

texts:
I have the following info for a feature class in ArcGIS that is supposed to be in WGS 1984/UTM 16N projected coordinate system:

  Projection: Transverse_Mercator

  False_Easting: 500000.000000

  False_Northing: 0.000000

  Central_Meridian: -87.000000

  Scale_Factor: 0.999600

  Latitude_Of_Origin: 0.000000

  Linear Unit: Meter (1.000000)

  Geographic Coordinate System:
  GCS_WGS_1984

  Angular Unit: Degree
  (0.017453292519943299)

  Prime Meridian: Greenwich
  (0.000000000000000000)

  Datum: D_WGS_1984

  Spheroid: WGS_1984

  Semimajor Axis: 6378137.000000000000000000

  Semiminor Axis: 6356752.314245179300000000

  Inverse Flattening: 298.257223563000030000

I also have a vector shapefile imported into QGIS that is also supposed to be projected WGS 1984/UTM 16N, with the EPSG 32616 and ID 3100.

Since I don't know what the EPSG and ID means, I am not sure if the shapefile in QGIS and the feature class in ArcGIS are actually in the same coordinate system.

Can anyone help with this? And if they are not in the same system, what should I do to correct that?


texts: Trying to get more classes to classify landcover on TensorFlow + Earth Engine

texts:
I used this notebook to create a neural network to classify landcover on Landsat 8 imagery.
https://colab.research.google.com/github/google/earthengine-api/blob/master/python/examples/ipynb/TF_demo1_keras.ipynb#scrollTo=Ejxa1MQjEGv9
The notebook is an example made for a 3 class classification, but I need to make an analysis with 4 and 5 classes. I changed the number of classes in the defining variables section
# Your Earth Engine username.  This is used to import a classified image
# into your Earth Engine assets folder.
USER_NAME = 'username'

# Cloud Storage bucket into which training, testing and prediction
# datasets will be written.  You must be able to write into this bucket.
OUTPUT_BUCKET = 'your-bucket'

# Use Landsat 8 surface reflectance data for predictors.
L8SR = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
# Use these bands for prediction.
BANDS = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7']

# This is a trianing/testing dataset of points with known land cover labels.
LABEL_DATA = ee.FeatureCollection('projects/google/demo_landcover_labels')
# The labels, consecutive integer indices starting from zero, are stored in
# this property, set on each point.
LABEL = 'landcover'
# Number of label values, i.e. number of classes in the classification.
N_CLASSES = 4

And I changed the way the network writes it's predictions in the final image -my guess is that my error is here-
 Instantiate the writer.
writer = tf.io.TFRecordWriter(OUTPUT_IMAGE_FILE)

# Every patch-worth of predictions we'll dump an example into the output
# file with a single feature that holds our predictions. Since our predictions
# are already in the order of the exported data, the patches we create here
# will also be in the right order.
patch = [[], [], [], [], []]
cur_patch = 1
for prediction in predictions:
  patch[0].append(tf.argmax(prediction, 1))
  patch[1].append(prediction[0][0])
  patch[2].append(prediction[0][1])
  patch[3].append(prediction[0][2])
  patch[4].append(prediction[0][3])
  # Once we've seen a patches-worth of class_ids...
  if (len(patch[0]) == patch_width * patch_height):
    print('Done with patch ' + str(cur_patch) + ' of ' + str(patches) + '...')
    # Create an example
    example = tf.train.Example(
      features=tf.train.Features(
        feature={
          'prediction': tf.train.Feature(
              int64_list=tf.train.Int64List(
                  value=patch[0])),
          'aguaProb': tf.train.Feature(
              float_list=tf.train.FloatList(
                  value=patch[1])),
          'vegnativaProb': tf.train.Feature(
              float_list=tf.train.FloatList(
                  value=patch[2])),
          'soloUrbanorProb': tf.train.Feature(
              float_list=tf.train.FloatList(
                  value=patch[3])),
          'agriProb': tf.train.Feature(
              float_list=tf.train.FloatList(
                  value=patch[4])),
        }
      )
    )
    # Write the example to the file and clear our patch array so it's ready for
    # another batch of class ids
    writer.write(example.SerializeToString())
    patch = [[], [], [], [], []]
    cur_patch += 1

writer.close()

In this example I put 4 classes up for classification, but when the final image is generated, the fourth class is always empty, not a single pixel classified as it.
Anything I'm missing here?


texts: Constant Solving by ArcGIS Network Analyst?

texts:
I have created a transportation network dataset from NYC streets.  Upon building and creating the dataset, no significant errors were found in my attributes, connectivity, etc. However, when I run any analysis, whether it be finding drive time service areas or a simple route from point 1 to point 2, I get a constant 'solving...' message and ArcGIS does not respond. I have thoroughly looked through my set up for the dataset and there does not seem to be any obvious errors.

My network dataset is stored within a File .gdb and on a USB drive with ample space.

Does anyone have some troubleshooting ideas or possible solutions?


texts: How to create an extent indicator in QGIS 1.8 print composer?

texts: I'm looking to put an inset map for my main map frame onto my map layout in QGIS print composer. I would like it to indicate the extent and location of the main map on a larger regional map. This something I am familiar with accomplishing in ARCMap but I have yet to discover if QGIS has an equivalent function. The indicator should show the extent of the data in the main frame in another data frame.


texts:
I'm new to ArcPy and am creating a script to automate some interpolation work.

I have multiple GA Layers outputted by an Empirical Bayesian Kriging interpolation (using arcpy.EmpiricalBayesianKriging_ga) which are in the in_memory workspace. I want to convert these GA layers to rasters for further manipulation.

Currently, I'm doing the following:

 for lyr in layers:
    print("Converting " + lyr + " to raster.")
    name = lyr + "_raster"

    arcpy.GALayerToRasters_ga(lyr, name, "PREDICTION")

Effectively reading each layer in from the in_memory workspace and converting them one at a time. This is fine for a couple layers, but it takes forever on larger collections of data (dozens or hundreds of layers).

Since the layers are independent from one another, this seems like a good application of multiprocessing to me. However, because the GA layers are all in the in_memory workspace, I'm under the impression that I will have to write each one to disk, read it back in, and then perform the GALayerToRasters_ga operation.

I want to do something like the below, but it predictably returns with "The layer does not exist."

def ga_layer_to_rasters_converter(layer):

    arcpy.env.workspace = "in_memory"
    arcpy.env.overwriteOutput = True

    print("Converting " + layer + " to raster.")
    name = layer + "_raster"

    # Names cannot be longer than 13 digits
    if (len(name) > 13):
        name = name[:13]

    arcpy.GALayerToRasters_ga(layer, name, "PREDICTION")

p = Pool(4)
p.map(ga_layer_to_rasters_converter, layers)
p.close()
p.join()

Do I have no choice but to write the GA layers to disk only to immediately read them back in?

I was hoping for a faster way to convert these GA layers to rasters.


texts: Non-admin GeoServer users with advanced privileges

texts:
I have user, who has full access to only their workspace. How can I allow him to:

change their password
set password protection for reading/writing layers from his workspace
create own service rules for own workspace/layers


texts: has anyone used their Garmin, i.e. Garmin 60CSX to grab field data instead of a Trimble or something similar?  I'm considering buying a Trimble for Field Data collecting but am wondering if you are able to use a 60CSX Garming for the short-term.


texts: Cannot select a feature using ArcObjects because coordinate systems differ

texts:
I'm using ArcMap 10 and Bing Maps Road as base map. Then we have our adding which is programmed using C#. What I'm trying to do is to select a feature (A) that intersects another feature (B, from a different feature class). I can achieve this with ArcMap's Select By Location function.

But when I try to do the same with ArcObjects, I cannot find anything because feature A seems to have a different coordinate system than feature B! Bing Maps Road uses WGS 1984 and my feature classes should use EUREF_FIN_TM35_FIN.

An interesting thing is, that when I manually modify feature class' geometry and save it, then both features are using EUREF (as they should in the first place)!

If I right click the Layers and take Data Frame Properties, it shows "WGS_1984_Web_Mercator_Auxiliary_Sphere" (and this is when error happens). If I change it to EUREF_FIN_TM35_FIN then everything works.

Code sample:

IFeatureClass featureClass = ((IFeatureWorkspace)m_workspace).OpenFeatureClass(name) as IDataset;

// Here I can see the feature.Shape.SpatialRef changes
IFeature feature = featureClass.GetFeature(oid.Value);

Now I would like to know is this a bug in ArcObjects or am I doing something wrong? I'm quite new to GIS and ArcMap, so there might be some simple thing I don't get.

In a nutshell:

I want to select feature B that intersects feature A
When I try to do that in ArcMap, it works
In ArcObjects it doesn't, except when I modify a geometry and save
It has something to do with the base map's coordinate system

Edit:
Here's how to reproduce this issue:

Preparation
On the map and in same data frame (data frame spatial ref wgs)

layer1 (from fgdb) - spatial reference EUREF
layer2 (bing maps) - spatial reference wgs

Start editor
Using arcobjects open layer1 from workspace opened using workspacefactory

((IFeatureWorkspace)workspace).OpenFeatureClass("name of layer1);

Find feature from layer1
layer1 feature geometry has spatial reference WGS
use geometry to find (spatial query) other object FAILS
Edit the geometry using editor and save edits
repeat steps 2-3
layer1 feature geometry has spatial reference EUREF
use geometry to find (spatial query) other object SUCCEEDS
Stop editor & Repeat from 1. (always repeatable)


texts: Using a csv or similar to limit possible parameter inputs

texts:
Is it possible to use the contents of a single column csv file to limit the list of acceptable inputs for a tool in a python toolbox? I'm aware that I can use the code below to achieve my requirements, but the list I need often varies outside of ArcGIS and outside of my control and therefore the code would need regular updates.

options1 = [1, 2]
param0.filter.type = "Value List"
param0.filter.list = options1


texts: Geoserver - Sql View's parameters not works in WFS service

texts:
I've created sql view which will be calculate route between 2 points which are pass over service ... so here is Sql View configuration window:

When I try preview layer I got this result on init:

{"type":"FeatureCollection","totalFeatures":0,"features":[],"crs":null}

So I've tried add parameters to url and execute it ... but with same result.
Here is my url:

localhost:8080/geoserver/Maps/ows?service=WFS
&version=1.0.0
&request=GetFeature
&typeName=Bachelor:Route
&maxFeatures=50
&outputFormat=application%2Fjson
&viewparams=s:30623;t:6021

I tried execute query with parameters into pgAdmin, and works .. but here not :/ any ideas?

Bounding box of layer will be dynamically by entered points, so I add full range into bbox inputs like this:


texts: How to find embedded Coordinate Reference system in a Layer

texts:
I want to find CRS in an added shapefile.
Ex. In ArcGIS if we go to Layer Properties and Source we can find embedded CRS to a layer.

How to find in QGIS, embedded CRS in a layer/ Shapefile?

I think in QGIS, if we move to Layer properties and General, it doesn't show extract CRS of layer.

Because I Added shapefile with undefined CRS. In that case it shows EPSG:32244 - WGS 72 / UTM zone 44N. But I added undefined one????

I'm using qgis 1.7.4,
System is Windows 7

How to find in QGIS, embedded CRS in a layer/ Shapefile???


texts:
I have a directory full of geojson files that I need to convert to shapefiles for editing in QGIS.  As far as I have searched, there are no resources in GIS StackExchange that explain how to do this from a  geojson .

The output should be shapefiles, named exactly the same as the geojson, in a separate directory. Can anyone recommend the best method to do this programmatically in Python or in other software?


texts: PostGIS Raster Tile Size Affects ST_Intersect?

texts:
My understanding of tiles in PostGIS rasters had been that they optimize the speed of calculations, so I was confused to run into an example where the choice of tile size actually changed the results of an intersect.

I am trying to find the value of myraster at the points in mypoints_tbl, but found that some points did not intersect with the raster, even though they appeared to when viewed in ArcMap.

In addition, I found that different tile sizes produced different numbers of points that don't intersect.

For example:

$ raster2pgsql -s 26958 -d -I -C -M myraster.tif -t 100x100 myraster_tbl | psql -q
$ psql -c "select count(*) from mypoints_tbl left join myraster_tbl on st_intersects(rast, geom) where st_value(rast, geom) is null;"

 count
-------
 8

$ raster2pgsql -s 26958 -d -I -C -M myraster.tif -t 1000x1000 myraster2_tbl | psql -q
$ psql -c "select count(*) from mypoints left join myraster2_tbl on st_intersects(rast, geom) where st_value(rast, geom) is null;"

 count
-------
 0

$ raster2pgsql -s 26958 -d -I -C -M myraster.tif -t auto myraster3_tbl | psql -q
$ psql -c "select count(*) from mypoints left join myraster3_tbl on st_intersects(rast, geom) where st_value(rast, geom) is null;"

 count
-------
 9

I can get more values of count for other tile sizes.

What am I misunderstanding about tile size and st_intersect in this case? It's not clear to me at all why one should affect the other.


texts: How to reclassify a raster based on quantile and without no data in QGIS3?

texts:
Revisiting this topic, I managed to use the answer given by https://gis.stackexchange.com/users/45066/xunilk in version 2.18, but this code does not run in version QGis 3.16, as shown below:
from qgis.core import QgsProcessing
from qgis.core import QgsProcessingAlgorithm
from qgis.core import QgsProcessingMultiStepFeedback
from qgis.core import QgsProcessingParameterVectorLayer
from qgis.core import QgsProcessingParameterField
from qgis.core import QgsProcessingParameterNumber
from qgis.core import QgsProcessingParameterRasterLayer
from qgis.core import QgsProcessingParameterRasterDestination
from qgis.core import QgsProcessingParameterFeatureSink
from qgis.core import QgsProcessingParameterVectorDestination
import processing
from qgis.processing import alg
from qgis.core import QgsFeatureSink
import numpy as np
from osgeo import gdal, gdal_array, osr
@alg("rastertoquantiles", "Raster to quantiles", group="rastertools", group_label="Raster tools")
@alg.input(type=alg.SOURCE, name="INPUT", label="Input raster")
@alg.input(type=alg.SINK, name="OUTPUT", label="Output raster")
dataset = gdal.Open('F:\QGis3_Plugins\Tutorials about Python Programming\CartaPerigo_01.tif')
band = dataset.GetRasterBand(1)
nodata = band.GetNoDataValue()
array = dataset.ReadAsArray()
new_array = array
nan_array = array
nan_array[array == nodata] = np.nan
percentile_80 = np.nanpercentile(nan_array, 80)
percentile_60 = np.nanpercentile(nan_array, 60)
percentile_40 = np.nanpercentile(nan_array, 40)
percentile_20 = np.nanpercentile(nan_array, 20)
percentile_0 = np.nanpercentile(nan_array, 0)
for i, v in enumerate(new_array):
for j, element in enumerate(v):
if element <= percentile_0:
new_array[i,j] = 1
if element > percentile_0 and element <= percentile_20:
new_array[i,j] = 2
if element > percentile_20 and element <= percentile_40:
new_array[i,j] = 3
if element > percentile_40 and element <= percentile_60:
new_array[i,j] = 4
if element > percentile_60 and element <= percentile_80:
new_array[i,j] = 5
if element > percentile_80:
new_array[i,j] = 6
new_array[ new_array != new_array ] = nodata
geotransform = dataset.GetGeoTransform()
wkt = dataset.GetProjection()
Create gtif file
driver = gdal.GetDriverByName("GTiff")
output_file = 'F:\QGis3_Plugins\Tutorials about Python Programming\CartaPerigo_01_Quant.tif'
dst_ds = driver.Create(output_file,
band.XSize,
band.YSize,
1,
gdal.GDT_Int16)
#writting output raster
dst_ds.GetRasterBand(1).WriteArray( new_array )
#setting nodata value
dst_ds.GetRasterBand(1).SetNoDataValue(nodata)
#setting extension of output raster
top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
dst_ds.SetGeoTransform(geotransform)
setting spatial reference of output raster
srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)
dst_ds.SetProjection( srs.ExportToWkt() )
#Close output raster dataset
dataset = None
dst_ds = None


texts: i'd like to create menus in a base toolbar like arcHdro tools or spatial analyst . how can i create menus with arcengine 9.3 in my toolbar ( the toolbar shows in arcmap )?


texts: Scatter Plot between 2 images Google Earth Engine

texts:
I have 2 high-resolution images in Google Earth Engine. I want to plot a scatter between these 2 images. But I am getting some error Missing required arguments to function seriesByRegion(): reducer. I don't know how to correct this (or if this is due to the images).

First I imported geometry and my Image:

var RZSC = ee.Image("users/chandrakant/Max_RZSC_Layer_Trail_3"),
    geometry = /* color: #d6cbbb */ee.Geometry.Polygon(
        [[[-81.375, -56.125],
          [-34.625, -56.125],
          [-34.625, 12.625],
          [-81.375, 12.625]]]);

Then I visulatized the Image

var vizParams = {
  bands: ['b1'],
  min: 0.0,
  max: 1500.0,
  palette: ['blue', 'green', 'red']
};
Map.setCenter(6.746, 46.529, 10);
Map.addLayer(RZSC.clip(geometry), vizParams, 'Rootzone Storage Capacity');
Map.centerObject(RZSC);
print('RZSC Projection, crs, and crs_transform:', RZSC.projection());

Visualized my MODIS image

var MODIStc = ee.ImageCollection('MODIS/051/MOD44B')
                  .filter(ee.Filter.date('2000-01-01', '2017-12-01')).mean();
var percentTreeCover = MODIStc.select('Percent_Tree_Cover');
var percentTreeCoverVis = {
  min: 0.0,
  max: 100.0,
  palette: ['bbe029', '0a9501', '074b03'],
};
Map.setCenter(6.746, 46.529, 2);
Map.addLayer(percentTreeCover.clip(geometry), percentTreeCoverVis, 'Percent Tree Cover');
print('MOD44B Projection, crs, and crs_transform:', percentTreeCover.projection());
\\Here I added Water Mask
var waterOcc = ee.Image("JRC/GSW1_0/GlobalSurfaceWater").select('occurrence'),
    jrc_data0 = ee.Image("JRC/GSW1_0/Metadata").select('total_obs').lte(0),
    waterOccFilled = waterOcc.unmask(0).max(jrc_data0),
    waterMask = waterOccFilled.lt(50);

Here I reprojected the image (RZSC) at 250 m

\\Here I tried to increase the resolution from 25km to 250m
var RZSC_250m = RZSC.resample('bilinear').reproject({
  crs: RZSC.projection().crs(),
  scale: 250
});
Map.addLayer(RZSC_250m.clip(geometry).updateMask(waterMask), vizParams, 'RZSC_250m');
print('RZSC_250m Projection, crs, and crs_transform:', RZSC_250m.projection());

Tried a scatterplot here

\\Tried to plot scatterplot here (I want just the basic scatter plot between 2 variables)
var chart = ui.Chart.image.seriesByRegion(RZSC, MODIStc,
                ee.Reducer.mean(),250, 'nd').setChartType('ScatterChart');
print(chart);

The error I get:
Error generating chart: Collection.map, argument 'collection': Invalid type. Expected: FeatureCollection. Actual: Image<[Percent_Tree_Cover]>.
and Collection.map, argument 'collection': Invalid type. Expected: FeatureCollection. Actual: Image<[Percent_Tree_Cover]>.

Additional information:
Also, I don't know why, but the dimensions of MODIStc and RZSC_250m are different.
For MODIStc:

bands: List (1 element)
0: "Percent_Tree_Cover", double ∈ [0, 255], EPSG:4326, 48x73 px
id: Percent_Tree_Cover
crs: EPSG:4326
crs_transform: [1,0,0,0,1,0]
data_type: double ∈ [0, 255]
dimensions: [48,73]
origin: [-82,-59]

For RZSC_250m:

bands: List (1 element)
0: "b1", float, EPSG:4326, 20818x28371 px
id: b1
crs: EPSG:4326
crs_transform: List (6 elements)
data_type: float
dimensions: [20818,28371]
origin: [-36235,-6107]

Also, when I try to use reproject on MODIStc, the data looks smoothed and weird (a lot of details are lost).


texts: How to use last clicked x/y coordinates in a qgis action?

texts:
I created an action (layer properties > actions) that opens an online map.

For the action I use this expression to dynamically generate the url of the map.

[% ürl" %]/?op=flyTo&cam=[% x($geometry)%],[% y($geometry) %],500.00&look=[% x($geometry)%],[% y($geometry)%],25.00

This works: the online map opens at the x/y position of the center of the object that you click.

For example:
https://pointscene.com/scene/85300c83/share/8a7f7c68262ccc1fbebdc38031f68dc4/?op=flyTo&cam=93618.37,436479.01,500.00&look=93618.37,436479.01,25.00

Question:

Instead of opening the online map at the center of the object that I clicked, I would like to open the map at the exact mapcoördinates that I clicked.

I noticed these coördinates are shown in the info tool window, as 'derived' data. But I can't find them in the expression builder.

So: how do I use the last clicked coördinates in my action?


texts:
How Does Erdas calculates the ground coordinates for each pixel in the image

I have a conceptual question regarding photogrammetry, what is the formula used by Erdas to find the (X,Y,Z) on ground. Erdas requires inputs like ( focal length, H of the plane, image exteriors (x,y,z, omega, phi, kappa), GCPs, pixel size) then the result is an orthophoto that have accurate coordinates, but I am curious on how does this software calculates the ground coordinate for each pixel in the image !
Is there any formulas/references on this, can I calculate a coordinate for one pixel manually just to check that the software is working fine or just to prove it for myself.


texts: Disable tiling when Mapserver requests layers from external WMS server

texts:
We use Mapserver to create a composite of multiple layers from another WMS server. We noticed that a single request to Mapserver creates multiple requests to the "upstream" WMS server. Apparently Mapserver uses smaller tiles (that are based on the layer "SizeX" and "SizeY" config) instead of using the request bounding box it received.

E.g. incoming request to Mapserver:

And what Mapserver requests from the upstream WMS server (4 tiles):

Is there a way to disable this behavior in Mapcache? It should just use whatever bbox and resolution it got in its request for the upstream query.

We are using Mapserver 6.4.1-3~precise2 on Ubuntu but will upgrade fairly soon.

Our Mapserver config of the layer in question:

CONNECTIONTYPE UVRASTER
DATA "<GDAL_WMS>
          <Service name=\"WMS\">
              <Version>1.1.1</Version>
              <ServerUrl>http://host/wms?&amp;TIME=%TIME%&amp;PROJECTION=%2Binit%3Depsg%3A3857&amp;ENCODER=ENC1&amp;</ServerUrl>
              <SRS>EPSG:3857</SRS>
              <ImageFormat>image/png</ImageFormat>
              <Transparent>TRUE</Transparent>
              <Layers>L1</Layers>
          </Service>
          <DataWindow>
              <UpperLeftX>-20037508.34</UpperLeftX>
              <UpperLeftY>20037508.340</UpperLeftY>
              <LowerRightX>20037508.34</LowerRightX>
              <LowerRightY>-20037508.34</LowerRightY>
              <SizeX>40075</SizeX>
              <SizeY>40075</SizeY>
          </DataWindow>
          <BlockSizeX>256</BlockSizeX>
          <BlockSizeY>256</BlockSizeY>
          <Projection>EPSG:3857</Projection>
          <BandsCount>4</BandsCount>
          <DataType>Byte</DataType>
          <MaxConnections>4</MaxConnections>
          <OverviewCount>10</OverviewCount>
      </GDAL_WMS>"


texts:
I have a list of objects called for sample "A1_YYYYMMdd, B1_YYYYMMdd, A3B3_YYYYMMdd, etc..."

For me the characters before the "" identify the name of the block and the numbers after the "" identify the date.

I want create a new "object" containing only the new one for each block.
es.

a= 'A1B1_20160213','A1B1_20161001','B2_20160202','B2_20161010
run script
b='A1B1_20161001','B2_20161010'

I try to write this but doesn't work.

##########################  TEST 1
import string
delivery=['A1B1_20160213','A1B1_20161001','B2_20160202','B2_20161010','F3_20160203']
dely=string.split(str(delivery),"_")
l=len(dely)
##print l
list(dely)
print type(dely)
d=[]
print type(d)
for i in range(0,l,2):
    print i
    if i==0:
        d=dely[i]
        list(d)
        print type(d)
    else:
        print type(d)
        list(d)
        print type(d)
        d.append(dely[i])
        print d
print d

##########################  TEST 2
import string
delivery=['A1_20160213','A1_20161001','B2_20160202','B2_20161010','B2_201610113','A2B2_20160202','A2B2_20161010',\
'F3_20160203','F3_20160204','F3_20160205','F3_20160206']

delivery.sort
l=len(delivery)
for i in range(0,l-1,1):
    print i
    a=delivery[i]
    b=delivery[i+1]
    aa=a[:2]
    bb=b[:2]
    if aa==bb:
        if a<b:
            c=b
            print c
        elif a>b:
            c=a
            print c


texts:
I am trying to reproject the raster data to "EPSG:4326", but getting this error :
RuntimeError: Translating source or target SRS failed: EPSG:4326
and my code is:
filename = 'mosaic_band1_1.tif'
input_raster = gdal.Open(filename)
output_raster = 'mosaic_band1_1_rep.tif'
warp = gdal.Warp(output_raster,input_raster,dstSRS='EPSG:4326')
warp = None

GDAL version 2.4.1
Driver: GTiff/GeoTIFF
Files: mosaic_band1_2.tif
Size is 7200, 4800
Coordinate System is:
PROJCRS["unnamed",
    BASEGEOGCRS["Unknown datum based upon the custom spheroid",
        DATUM["Not_specified_based_on_custom_spheroid",
            ELLIPSOID["Custom spheroid",6371007.181,0,
                LENGTHUNIT["metre",1,
                    ID["EPSG",9001]]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]],
    CONVERSION["Sinusoidal",
        METHOD["Sinusoidal"],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["northing",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
Data axis to CRS axis mapping: 1,2
Origin = (3335851.558999999891967,6671703.117999999783933)
Pixel Size = (463.312716527916677,-463.312716527916507)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:

Upper Left  ( 3335851.559, 6671703.118) ( 60d 0' 0.00"E, 60d 0' 0.00"N)
Lower Left  ( 3335851.559, 4447802.079) ( 39d 9'43.99"E, 40d 0' 0.00"N)
Upper Right ( 6671703.118, 6671703.118) (120d 0' 0.00"E, 60d 0' 0.00"N)
Lower Right ( 6671703.118, 4447802.079) ( 78d19'27.97"E, 40d 0' 0.00"N)
Center      ( 5003777.339, 5559752.598) ( 70d 0'27.26"E, 50d 0' 0.00"N)
Band 1 Block=7200x1 Type=Int16, ColorInterp=Gray
  NoData Value=255
Band 2 Block=7200x1 Type=Int16, ColorInterp=Undefined
  NoData Value=255
Band 3 Block=7200x1 Type=Int16, ColorInterp=Undefined
  NoData Value=255


texts: Graduated symbol for percentage attribute in QGIS

texts: I want to use the Graduated Symbol style for a field which contains values in percentages. When I select the Graduated Symbol this field does not appear in the drop down box. Is there a way to make this field appear?


texts: I got a temperature raster (maximum temperature for January, global extent) and the value range seems to go beyond the valid range. The highest is 42.006 (which I assumed should be Celsius) but the lowest is -9999 (which corresponds to the ocean). The image is attached.


texts: Categorizing multiple values_multiple atributtes - QGIS 2.14

texts:
I have a problem with a polygon layer.

Its an equipment and services layer... I want to aggregate different types of services in similar legend, like this:

For example, these four :

"Ensino Especial"
"Ensino Básico"
"Creche/Jardim de Infância/ATL"
"Ensino Secundário"

I want to label all of them as "Equipamento de Educação" (Education Equipment).

Do I have to use "Ruled-Based"?
Do I write the RULE in every row?

By-hand, what I get is this:

Duplicated entries... When I switch to print composer...I get at least four "Equipamentos de Educação" :/

But I want one label for "Equipamentos de Educação" that conjugates all of these: "Creche/Jardim de Infância/ATL"; "Ensino Especial"; "Ensino Básico"; "Ensino Secundário"..."


texts: How to clip polygon such that only regions in 'line of sight' from a point remain

texts:
I am using ArcGIS.

My point represents a monitoring station and is located within a polygon representing a river catchment.

I have created a buffer of 300m representing the detection range of the monitoring station and I have clipped this to the polygon area using the intersection tool.

I want to further clip the buffer region so that it only contains the polygon area in 'line of sight' from the monitoring station.

The image shows the GIS output I have so far achieved, with the clipped buffer region in purple, but with areas marked in red that are not in line of sight of the monitoring station and should therefore be removed. Is there a tool to achieve this, perhaps with radiating lines from the monitoring station?


texts: Extracting data to table for two rasters with identical polygons (different values) using ArcGIS Pro

texts:
I have a raster image of watershed catchments where the raster value is the catchment ID. Using this raster and the "Zonal Statistics", I created a separate raster of land use where the raster value is land use. In each raster, there are identically shaped catchments. How do I generate a table of catchment and land use based on these images? I have ArcGIS Pro, it's just not clear what tool I would use.
Catchment---Land Use
020506.........4
95138..........12
etc.


texts: Grouping vector file column attributes under another column's attribute as a hierarchy group

texts:
I would like to output attributes in a vector file I have for Land use. There are three categories that are in hierarchical order. The columns represented in the attributes table are:
PRIMARY_20 > SECONDARY_ > TERTIARY_2
I would like to display the TERTIARY_2 column attributes but group them under their PRIMARY_20 category so it looks like this:
PRIMARY_20 (no symbology, just group heading)
(grouped below PRIMARY_20) TERTIARY_2 (with symbology displayed in the map layers panel and in the map legend (layout))
I have tried the method outlined here (I added the SECONDARY column just as a test): https://gis.stackexchange.com/questions/34428/dissolve-column-based-on-attribute-in-another-column[![enter image description here]2]2
and my own method trying to concatenate (see figure attached).
I can get the column attributes to join but cannot space them ( I have tried ' ' between concatenate) and tried to insert a new line before TERTIARY_2 (using '\n' function) to no avail.
I have attached figures to try and explain visually.


texts:
I have the following RasterLayer containing categorical values for each pixel. I am using the following code to plot it

par(xpd)=FALSE
plot(LC, legend=FALSE, col=c("lightcyan3","royalblue","forestgreen","coral3","papayawhip"),xaxt='n', yaxt='n')

par(xpd=TRUE)
legend("bottom", legend=c("Urban","Water","Forest","Agriculture", "Bare soil"),fill=c("lightcyan3","royalblue","forestgreen","coral3","papayawhip"),horiz = TRUE, inset=-0.175)

In total I have 5 classes but due to the configuration of the legend only 4 are swhon. I want all the elements of the legend to be closer. I have tried to used text.width=c(1,1,1,1) for example to set the distances but it is not doing the job.

I am using a similar approach as the one proposed here and was considering also to use the rasterVis package


texts: Consuming ArcGIS Online services from ArcGIS Desktop 10.0?

texts:
I work at a county government and was at one of our municipalities recently.
I was planning to show them how to consume our cached tile service from our most recent flyover that I have hosted on AGOL. I brought up AGOL from the ArcMap menu and searched for the relative tags. It showed my service in the results, but nothing happened when I clicked the "Add" option.  The tile service is shared to everyone so no sign in was necessary.

I did notice that they're still running ArcGIS 10.0.  Does anyone know if there is an issue with 10.0 not being able to consume certain services? I'm just trying to troubleshoot why it wasn't working for them. Their computer has Internet access.


texts: Which interpolation method is used in the close gaps module of Saga GIS?

texts:
I like to close gaps in a raster layer using the Module Close Gaps of Saga. I can not find which interpolation method is used by the module. Which method is used by the module and/or where can I find out?

http://www.saga-gis.org/saga_tool_doc/2.1.3/grid_tools_7.html


texts: Saving array of objects in PostGIS field?


texts: Why this pyqgis code cause crash dump in QGIS?

texts:
In ArcGIS there is gp.EucDirection to calculate Euclidean allocation, distance, and direction around a point, Is there equivalent function in pyqgis or gdal?

I wrote a script with help of gene and ndawson to create Euclidean direction raster in pyqgis, this script is ok when raster file is small, but when I want to use it in a large area (for example 10000*7000 pixel) it causes crash dump and close QGIS, Is there a solution for this problem?
link of script: Euclidean Direction


texts: Rotation of Text in a layout to be dependent on arbitrary Map

texts:
I have a Print composer with two map windows (MainMap and InsertMap).  I also have a bit of text which I want to rotate on the Composer and it to be dependent upon the rotation of MainMap what ever value it may take.

I believe I should be able to use rotation.  The question is how do I go about exposing the

@map_rotation

Variable associated with MainMap and using it in an expression associated with the text element in the composer and applying it to the rotation of the text element.

I am using QGIS 2.14.6 on a windows 7 box.

The use case is I am using a font symbol for my north arrow ('H' Esri North if you are interested...)


texts: Problem when creating database based on template_postgis

texts:
I have created a PostGIS template but when I try to create a database based on this template this is what happens:

Connecting to the template:

\c template_postgis;

You are now connected to database "template_postgis" as user "postgres".

\dt
               List of relations
  Schema  |      Name       | Type  |  Owner
----------+-----------------+-------+----------
 public   | spatial_ref_sys | table | postgres
 topology | layer           | table | postgres
 topology | topology        | table | postgres
(3 rows)

create database prueba with TEMPLATE template_postgis;
CREATE DATABASE

Connecting to database prueba:

\c prueba;
You are now connected to database "prueba" as user "postgres".

prueba=# \dt
              List of relations
 Schema |      Name       | Type  |  Owner
--------+-----------------+-------+----------
 public | spatial_ref_sys | table | postgres
(1 row)

You can see only the public schema was copied, no topology schema layer or topology were copied.

I had no errors when I created the template:

#createdb template_postgis
      #createlang plpgsql template_postgis
      #psql -d template_postgis -f /usr/pgsql-9.2/share/contrib/postgis-2.1/postgis.sql
      #psql -d template_postgis -f /usr/pgsql-9.2/share/contrib/postgis-2.1/postgis_comments.sql
      #psql -d template_postgis -f /usr/pgsql-9.2/share/contrib/postgis-2.1/spatial_ref_sys.sql
      #psql -d template_postgis -f /usr/pgsql-9.2/share/contrib/postgis-2.1/rtpostgis.sql
      #psql -d template_postgis -f /usr/pgsql-9.2/share/contrib/postgis-2.1/raster_comments.sql
      #psql -d template_postgis -f /usr/pgsql-9.2/share/contrib/postgis-2.1/topology.sql
      #psql -d template_postgis -f /usr/pgsql-9.2/share/contrib/postgis-2.1/topology_comments.sql

Could you please help me?


texts: ArcGIS JavaScript API Graphic doesn't Display

texts:
I am using map.graphic layer to present a query result, which is a featureset. I assigned symbol and template to the graphics and added them to map.graphics. However, when I apply the function, it doesn't display the graphics. I added another single graphic to the map and it presents well. Can anybody help me with the issue?
Here is the webpage:
http://www.public.asu.edu/~jzhan300/HCIC/HCICInteractiveMap.html
When you hit the select, and then click the apply it will execute the function.
Here is the query and display part of code:

var whereExpression = "";
myQueryTask = new QueryTask("http://arcgis-cabhpgis-2026970499.us-west-1.elb.amazonaws.com/arcgis/rest/services/HCIC/160524_interactive_map_jzz_v02/MapServer/0");
myQuery = new Query();
myQuery.returnGeomery = true;
myQuery.outFields=["*"];
on(dom.byId("applyServiceSelection"),"click",executeQuery);
function executeQuery(){
    console.log("executequery started");
    whereExpression = "";
    for(i=0;i<serviceList.length;i++){
        if(i==0) {
            whereExpression += ("ServiceProvided LIKE '%" + serviceList[i] + "%'");
        }
        else{
            whereExpression += (" AND ServiceProvided LIKE '%" + serviceList[i] + "%'");
        }
    }
    myQuery.where = whereExpression;
    console.log("where statement completed in executequery function");
    myQueryTask.execute(myQuery,showResults);
    console.log("myquerytask completed in executequery function");
}
var symbolPoint = new SimpleMarkerSymbol(SimpleMarkerSymbol.STYLE_X, 10,
    new SimpleLineSymbol(SimpleLineSymbol.STYLE_SOLID, new Color([255, 0, 0]), 1),
    new Color([0, 255, 0, 0.25]));
function showResults (results) {
    console.log("showresults started");
    console.log("total graphics in graphics layer: " +mainMap.graphics.graphics.length);
    //mainMap.graphics.clear();
    dom.byId("info").innerHTML = results.features.length;
    for (var i=0; i<results.features.length; i++){
        var aGraphic = results.features[i];
        aGraphic.setSymbol(symbolPoint);
        aGraphic.setInfoTemplate(BHInfoTemplate);
        mainMap.graphics.add(aGraphic);
        console.log("one feature launched, total graphics in graphics layer: " +mainMap.graphics.graphics.length);
        console.log(aGraphic.attributes);
    }
    console.log("showresults finished");
}
mainMap.on("load", function () {
    var point = new Point([-112,34.3]);
    var pointGraphic = new Graphic(point, symbolPoint);
    mainMap.graphics.add(pointGraphic);
    console.log("one feature launched, total graphics in graphics layer: " +mainMap.graphics.graphics.length);
});

You can check the full code from the link.


texts: Add / collate QGIS layer to existing PostGIS table?

texts:
With the database manager in QGIS, it is possible to export a QGIS layer to a PostgreSQL database, by creating a new table in the database. For instance, if I have a cities layer in my QGIS project, I can export it to a new table cities in my PostgreSQL database.
The tool provides the option of replacing the destination table in the database, if this table already exists (see screenshot below, sorry about the French locale).
Now, my question is: is there any way to "concatenate" the data of my QGIS layer cities to the existing citiestable of my PostgreSQL database, instead of erasing the existing table and replacing it by a new one?
(assuming of course, that the existing database table and the QGIS layers have the same CRS and the same geometry type)


texts: Is the postgis layer option in qgis 2.12 Lyon two way? Was the spit plugin built in?

texts:
I am trying to install Postgresql with Postgis on a remotely accessed terminal only installation of Ubuntu 14.04 that I want to use to connect to QGIS 2.12 and hopefully have a two way transfer between editing data and committing it to my database and grabbing it to edit.

If so, is there any known compatibility issues if I use Postgresql 9.5 and Postgis 2.1 given that QGIS 2.12 is only a month old?

Essentially, I would like to be able to edit data in QGIS and have a Postgis database connected to Geoserver for display.


texts: GROUP BY on geography in postgres 9.4.0 returns wrong/different results

texts:
table definition:

CREATE TABLE search (
  id bigserial NOT NULL,
  location geography(Point,4326) NOT NULL,
  location_aggregated geography(Point,4326) NOT NULL,
  created timestamp with time zone NOT NULL,
  CONSTRAINT search_pkey PRIMARY KEY (id)
)

The table was filled with 1699012 rows containing random geography data around some points in a regional area for the month february. If that is important: i wrote a small java program to fill in data using jdbc and the postgis jdbc lib for the insert with org.postgis.PGgeometry class for the location. location_aggregated is the result of ST_SnapToGrid using the value of location and 0.005.

I'm using this statement:

SELECT (location_aggregated), count(*)
FROM search
where created between '2015-02-21 00:00:00' and '2015-02-28 23:59:59'
group by (location_aggregated)
having count(*) > 1
order by 2 desc;

The result is then 104 rows with all having count of 2. But there are much more results for the same geography point.

If i'm testing that with e.g. ST_GeogFromText('SRID=4326;POINT(11.05 49.47)') i get 70 hits in the specified created range.

If i use ST_AsText(location_aggregated) in the GROUP BY instead, or just if i add location_aggregated = ST_GeogFromText('SRID=4326;POINT(11.05 49.47)') to the WHERE clause the result is very different. Using where it returns 1 row with the correct count of 70, using the function in GROUP BY i get 39010 rows including the one with the count of 70 mentioned before.

This looks like the GROUP BY on a geography (or just this special definition) doesn't read all rows to aggregate or the value just isn't read with the same precision for both ways of access. Checking with ST_AsText and also converting to a varchar/string does return the same values for the 70 rows of the example point.

Is it a bug or do i have to use a function on the geography column for grouping?


texts:
Now I have a lot of files that contain XY coordinates, like the picture, but this picture just contains one day. However, I want to integrate all the files I have.

Secondly, I want to integrate all point files. This is the main problem I have. The "Merge" or "Append" tool did not work for me at all or I did not know how to use it correctly, so I can not integrate all point files.
Actually, this is the model I would like to do but it always pops up the error. It says "You start with the invalid character".


texts: How to convert data from a .gdb into a shapefile without ArcMap?

texts:
I have a .gdb folder with LOCK Files, FREELIST Files, ATX Files, GDBINDEXES Files, GDBTABLE Files, GDBTABLEX Files and, SPX Files.

I want to convert the parcel data within this folder into a shapefile so that I can load it into Arc GIS online. If I still had my ArcGIS desktop version this would be a simple task however, I no longer have access to this version of the software.

Any ideas on how to get this parcel data into a shapefile? I haven't tried downloading QGIS yet so I'm not sure if this software would have that capability as ArcGIS 10.1 did


texts: Why there is no "create network dataset" option when I right click my shapefile?

texts: I wanted to use Network Analyst toolbox to do a vehicle routing problem in ArcGIS. When I right click the related shapefile which includes links of the network, there is no "create network dataset" option there. I created a network before but that was not the one I wanted to use. Now that I want to create a new network dataset, it does not let me do that. when I chooses "a new vehicle routing problem" from Network Analyst toolbar, it seems that it is using the information from the previous network dataset that I created. I am confused. I wonder if I should remove the previous network dataset so it let me create new network datasets. I even removed the previous network dataset but it seems that the toolbox is still considering that.


texts: How to translate ESRI cursor example to pyQgis?

texts:
For example, how to translate this code into pyqgis?

    rows = gp.UpdateCursor("D:/St_Johns/data.mdb/roads")
    row = rows.Next()

    while row:
       row.buffer_distance = row.road_type * 100
       rows.UpdateRow(row)
       print row
       row = rows.Next()

errata:
buffer_distance is a calculate column
road_type is a calumn with calculate value

Follow your advice I type this code in to python windows:

>>> from qgis.core import QgsVectorLayer, QgsFeature
>>> layer = QgsVectorLayer(r"D:\fold", "boundingBoxes.shp", "ogr")
>>> road_type_index = layer.fieldNameIndex("road_type")
>>> buffer_distance_index = layer.fieldNameIndex("buffer_dis")
>>> layer.select(layer.pendingAllAttributesList())
>>> layer.startEditing()
True
>>> for feature in layer:
...     newvalue = feature.attributeMap()[road_type_index].toInt()[0] * 100
...     feature.changeAttribute(buffer_distance_index, newvalue)
...     layer.updateFeature(feature)
...

next when I press "enter" the qgis (1.7.4) crash and exit


texts: Grouping points with distance using ArcGIS Desktop?

texts:
I have a point layer and I want to create polygons arround points which are in distance of 10 meters.

I did create a grid vector of 10m but it doesn't group as well as they have to be.

I did create a buffer arround my points in arcgis, as you know arcgis buffer tool gives you the posiibility to choose a dissolve type yet in the merged buffers I can have also a point which is in a distance of 20m.

I'm wondering if there is a possibility, using ArcGIS Desktop, to select from the same layer the features within a distance of  10m  and create a polygon around.


texts: Calculate certain records (Calculate Field)

texts:
I have records that need to be updated.
For this I select all records that are false and then process them with Calculate Field. The other records should not be changed.
The output has NULL in each record. I don't understand where the error is.
I also tried instead of return ags else: continue.
#select records for updating
arcpy.management.SelectLayerByAttribute("Strassenverzeichnis", "NEW_SELECTION", "AGS_NR_KOMMUNAL = '03154003' Or AGS_NR_KOMMUNAL = '03154010' Or AGS_NR_KOMMUNAL = '03154011' Or AGS_NR_KOMMUNAL = '03154020' Or AGS_NR_KOMMUNAL = '03154023' Or AGS_NR_KOMMUNAL = '03154401' Or AGS_NR_KOMMUNAL = '03156003' Or AGS_NR_KOMMUNAL = '03156012' Or AGS_NR_KOMMUNAL = '03156013' Or AGS_NR_KOMMUNAL = '03156016' Or AGS_NR_KOMMUNAL = '03156017' Or AGS_NR_KOMMUNAL = '03158015' Or AGS_NR_KOMMUNAL = '03158024' Or AGS_NR_KOMMUNAL = '03158026' Or AGS_NR_KOMMUNAL = '03158029' Or AGS_NR_KOMMUNAL = '15085010'", None)

#locale variables
inTable = "Straßenverzeichnis"
fieldName = "AGS_NR_KOMMUNAL"
expression = "upd(!AGS_NR_KOMMUNAL!)"
codeblock = """
def auf(liste):
    neu_liste = []
    for i in liste:
        neu_liste.append(upd(i))
    return neu_liste

def upd(ags):
    if ags == "03154003":
        return "03154028"
    elif ags ==  "03154010":
        return "03154028"
    elif ags ==  "03154011":
        return "03154027"
    elif ags ==  "03154020":
        return "03154027"
    elif ags ==  "03154023":
        return "03154027"
    elif ags ==  "03154401":
        return "03154016"
    elif ags ==  "03156003":
        return "03159004"
    elif ags ==  "03156012":
        return "03159036"
    elif ags ==  "03156013":
        return "03159036"
    elif ags ==  "03156016":
        return "03159036"
    elif ags ==  "03156017":
        return "03159002"
    elif ags ==  "03158015":
        return "03158039"
    elif ags == "03158024":
        return "03158040"
    elif ags ==  "03158026":
        return "03158039"
    elif ags ==  "03158029":
        return "03158040"
    elif ags ==  "15085010":
        return "15085330"
    return ags
"""
arcpy.CalculateField_management(inTable, fieldName, expression, "PYTHON3", codeblock)


texts: is there a way to save a temporary file as a permanent raster layer

texts: I am trying to build a raster table for a temporary layer i created from an image difference on ArcGIS image analysis. Is there a way to save the temporary layer from this image analysis as a permanent .img file/hdf file or something along those lines ?


texts: Is it possible to publish a ArcGIS .mxd file in OpenLayers?

texts: It is possible to publish a .mxd in OpenLayers? I am using ArcGIS 10.0.


texts: Putting frame around legend using QgsComposerLegend?

texts:
I can already see this being a simple answer but I'm stumped with this. I want to create a frame around a legend created by a PyQGIS script

Code snippet is:

legend = QgsComposerLegend(c)
c.addComposerLegend(legend)
legend.model().setLayerSet(map_settings.layerSet())
legend.setTitle('Clay %')
legend.setItemPosition(10,200, 0, 0)
legend.updateLegend()
c.addItem(legend)

I've looked through the API expecting to find something along the lines of legend.enableFrame() but am at a loss.

Any ideas?


texts: Install GRASS GIS extension r.pi including r.pi.prox and implement the Addon in QGIS

texts:
I would like to install the GRASS GIS addon r.pi within QGIS,

I only have access to basic GRASS commands

I installed QGIS using the osgeo4w-setup-x86_64.exe advanced installation using the http://osgeo4w-oslandia.com

Downloadsite and I installed GRASS GIS 7.4.1 where I was able to install four r.pi commands (r.pi.cor.mw, r.pi.grow, r.pi.lm, r.pi.rectangle) using following command:

g.extension extension=r.pi -s

but I need to compile the r.pi.prox (Calculation of Proximity and Modified Proximity Index for patches in buffer region).

The command:

g.extension extension=r.pi.prox -s

returns:

  ERROR: Extension  not found

And I don't know how to implement the extensions (Addons) in QGIS.


texts: ArcPy replacing network mapped drive with UNC

texts:
I am writing a script in Python to look in the current MXD for paths which are on a mapped network drive instead of the UNC path to the network.
However I am getting a syntax error. I am not very versed in Python.
How can I fix this script?
import arcpy
mxd= arcpy.mapping.MapDocument
for lyr in arcpy.mapping.ListLayers(mxd):
    if lyr.supports("DATASOURCE"):
        if r"company shared data" in lyr.datasource.lower():
            lyr.findAndReplaceWorkspacePath(r"S:",r"\\network.local\company shared data")
mxd.saveACopy(r"\\network.local\company shared data\Scientific Data\ArcGIS\MXD")
del mxd


texts: If I fill a field through a dropdown menu in QGIS, will the content still be there in after importing the shapefile into a different GIS programme?

texts:
Is the content of the field in the attribute table turned into "normal" text even if entered through a dropdown menu once I saved the shapefile in QGIS?
So that I may open it with a different GIS programme and it will work fine or will it be empty, because I don´t have the dropdown menu in there?


texts:
I would like to take the data generated in the doySeries chart and use it as an array of values. Is this at all possible? What I am interested in is using the mean over the years to create a list of cumulative precipitation values that I can then use in another chart as a reference point to identify anomalous days in later years. Below is the script I have so far:

var polygon = ee.Geometry.Polygon([[-92.56227555674951,40.858223918075275],
                                   [-92.56227555674951,40.85666590660416],
                                   [-92.56227555674951,40.85666590660416],
                                   [-92.55961480540674,40.85828883442468]]);

var startDateEarly = ee.Date('1981-01-01');
var endDateEarly = ee.Date('1983-12-31');

var chirpsColl = ee.ImageCollection('UCSB-CHG/CHIRPS/DAILY')
                .filterDate(startDateEarly, endDateEarly)
                .filterBounds(polygon);

var chirpsMean = ui.Chart.image.doySeries({
  imageCollection: chirpsColl,
  region: polygon,
  scale: 1,
  yearReducer: ee.Reducer.mean(),
  //seriesProperty: 'label'
});
print(chirpsMean);

on a side note, is the scale a buffer? If so, if I set it to 0, will it only reduce within my selected region?


texts: PostgreSQL PostGIS: Spatial join but keep all features that don't intersect

texts:
Spatial joining a point table to a polygon table, creating a new output table. So far I can only get the output table to contain features that satisfy the spatial join. However, I want the output to contain ALL features, even the ones that don't intersect. I also want to choose which columns are pulled through into the final output (in the case below, 5 columns from 'topo' and only 1 from 'adb')

Here's my code as it stands at the moment:

DROP TABLE IF EXISTS public.wolv_topobuild_uprns_all;
CREATE TABLE public.wolv_topobuild_uprns_all AS
SELECT
  topo.id,
  topo.geom,
  topo.ogc_fid,
  topo.theme,
  topo.class,
  adb.uprn
FROM
  public.wolv_adb_slim as adb,
  public.wolv_topobuild as topo
WHERE
  ST_WITHIN(adb.geom, topo.geom)

From my research it seems I want to mix this spatial join with a full outer join, but being able to specify the columns I want in the output table.

I want the resulting output table to contain all of the original polygons but with the chosen metadata from the points where they fall inside a polygon. There will be some points that lie outside of the polygons... we can ignore them from the final output.


texts:
I have a shapefile with many polygons that are bigger than my working extent. I want to rasterize only fields that are fully contained in my working extent.

How can I do this?


texts: Adding special characters (e.g. $) to double value in Identify window?

texts:
I don't know what I want is possible in ArcMap? I have a field whose type is "Double". The value in it is accurate to 2 decimal places. When I use Identify to see the detail information of a feature, I wonder whether there is a way to add a $ in front of the double value?

What I see in Identify is as follows:

What I intend to see in the CONTRACT_SUM field is $24222568, is this possible?


texts: Expression to select features by comparing attribute values

texts:
I'm working in QGIS and I have a point layer with two fields and I'm trying to compare the accuracy of the distances in feet. One field highlights the old calculated distances and another field highlights the new.

How do I select all the features in this point layer that have a plus or minus change of less than 50 (feet)?

Both fields are already in feet unit so I do not need any statements including conversions.

Does anyone know an expression to do this selection?


texts: How to use QgsRasterFileWriter to save and reproject raster?

texts:
I have to save georeferenced png images as GeoTiff and reproject them from different coordinate systems (UTMs ang Google) to given UTM. I want to do this with pyton code and tried to use QgsRasterFileWriter.

I face two problems here:

How to save file with 4 bands (transparency band)
How to reproject file at the same time

For now my code does the first thing and looks like this:

# Load Raster
fileName = myRaster.png
baseName = myRaster
rlayer = QgsRasterLayer(fileName, baseName)

crsSrc =  rlayer.crs()  # may be any
crsDest =  QgsCoordinateReferenceSystem(4326)  # WGS84 or UTM Zone

renderer = rlayer.renderer()
provider = rlayer.dataProvider()

pipe = QgsRasterPipe()
pipe.set(provider.clone())
pipe.set(renderer.clone())

file_writer = QgsRasterFileWriter(baseName + '_TEST.tif')
file_writer.Mode(1)

file_writer.writeRaster(pipe, provider.xSize(), provider.ySize(), provider.extent(), provider.crs())

I tried with using QgsRasterProjector but never got it right.


texts: Converting TopoJSON data to GeoJSON using JavaScript in browser?

texts:
There are many articles about how convert GeoJSON to TopoJSON using JavaScript, but I could not find how to convert TopoJSON to GeoJSON in the browser.

Response:

Using http://d3js.org/topojson.v1.min.js, TopoJSON API’s
feature() method we can convert a high-level TopoJSON object back into a GeoJSON structure:

geojson = topojson.feature(topojson,topojson.objects[topologyKey]);


texts: Getting pixel value of GDAL raster under OGR point without NumPy?

texts:
I am working on a computational model of the abundance of wild pollinators across a landscape.  The model itself is complete, and I am now struggling with a post-processing step.

I have my GDAL pollinator supply raster that looks something like this (lighter colors mean higher pollinator visitation to a pixel):

And I have an OGR shapefile of points representing sample locations on the landscape:

I'm trying to run some analysis on the pixels under these points, but to do so, I need to be able to extract the value of a pixel under a point.

Is it possible to extract the value of a pixel under a point using only OGR and GDAL through Python?  I would prefer to avoid reading the entire raster into memory through ReadAsArray(), as my output rasters are very, very large (too large to fit into memory).

I noticed this post, which is similar, but requires a command-line call.


texts: Turn off Radiometric Enhancement ("DRA") when downloading Worldview-2 imagery

texts:
I am trying to calculate the TOA reflectance for worldview-2 images, but every time I download the WV-2 8 band imagery from DigitalGlobe, the metadata states that the image has been radiometrically enhanced and the absolute calibration factor is -9999, which is obviously not correct. Does anyone know how to set radiometricEnhancement to "off"? I can not figure out how to do this.

radiometricLevel = "Corrected";
radiometricEnhancement = "DRA";
bitsPerPixel = 8;
compressionType = "None";
BEGIN_GROUP = BAND_C
    ...
    absCalFactor = -9.990000e+02;


texts: Rectangles Ovals Digitizing Plugin error in QGIS

texts:
I am using Rectangles Ovals Digitizing plugin in QGIS 2.18.2 for macOS. What I noticed is when using the Rectangle by extent tool it creates two layers so I have duplicate layers or features with duplicate geometries.
Is this a bug or I have to change something on shapefile layer properties?

Note: To remove duplicate I am currently using MMQGIS> Delete Duplicate Geometries. But I want to remove the error while creating the boxes.
Also any other plugin which can create Bounding boxes/ Rectangles


texts: How to view "quick" statistics for a table in ArcGIS Pro

texts:
In ArcMap i can view statistics describing the values in numeric columns of the attribute table by right-clicking the field name and then click Statistics:

I cant find statistics like this in ArcGIS Pro, instead it will open the geoprocessing tool "Summary Statistics". How can i display quick statistics for a field in ArcGIS Pro?


texts: How to digitize polygons from georeferenced raster file

texts:
I used QGIS to georeference a raster file of a custom US map. The map has features which are not strictly defined by state borders. I would like to digitize the map by creating contiguous borders for each feature. I tried that a few different way, but none of them really worked:

traced the borders with polylines, and digitizing with the Polygonizer plugin. The plugin kept crashing.
used Conversion > Raster > Vectorize (Raster to Vector). That one worked, but the result was many overlapping polygons for each feature, each with slightly different borders.

I have somewhere between 50-100 different raster images to process, which is why any automated solution is highly preferred. I'm new to GIS software, but could easily pre-process the images if needed, either manually or with a custom script (using Python/PIL, ImageMagick, Photoshop batch processing or similar tools).

Thanks for any pointers!

Original

Georeferenced image


texts: How to calculate area in sq m - troubleshooting projection issues in a global dataset in QGIS

texts:
I am using QGIS version 2.12.3.

I have 3 polygon shapefiles: A, B and the intersect of A and B, which we can call C.

All I want is to add a field to the attributes table of each shapefile giving the area (in square m or square km) of each feature in that layer.  A, B and C each contain upwards of 10,000 individual features (polygons of protected areas), distributed over every continent.

I have used a global projection (WGS 84/NSIDC EASE Grid-Global EPSG:3975, have also tried 3410), which is listed under 'Projected Coordinate Systems' so I assume it is projected and not geographical, and set this as the coordinate system for the project and each one of the layers.

When I use the field calculation tool and the $area function, the resulting field values are definitely not in square m or square km and bear no resemblance to area as I understand it.

I have run it several times having fiddled with settings that I'm fairly sure I don't understand, and I've ended up with either really tiny values (e.g. 0.0000002345) or random sets of values including negative numbers.

Can anyone identify where I've gone wrong and suggest how I can go about completing this fairly simple task?


texts: Counting number of pixels using ArcGIS Spatial Analyst?

texts:
There is a raster data which contains specific areas separated by layers. My purpose is to count the number of the pixels in those areas.

Is it possible to do it in the ArcGIS Desktop 10.5.1 with a Spatial Analyst license?


texts: How can I use gdal warp on a .sid to create an .tif of the same resolution?

texts:
I’m trying to convert a .sid to .tif in OSGeo and was left this script by a former co-worker. I’m not very savvy in scripting, so I’m unsure of my translation.
First, I set the directory. The conversion command includes "-tr  ": set output file resolution (in target georeferenced units). Does “target georeferenced units” here mean the units of the source data? The .sid I am trying to convert is in degrees.
If I wanted the output to be the same resolution as the source, wouldn’t the following numbers be 1 and 1? Or is it asking what the size (in the units specified in the source data) of each pixel should be in the output?
If it’s the latter, then if the source data was in degrees, with 30 cm resolution, this script:
for /r %a in (*.sid) do gdalwarp -tr .000005 .000005 %a %~dnpa.vrt && for /r %b in (*.vrt) do gdal_translate -co COMPRESS=JPEG -co TILES=YES %b %~dpnb.tif && del %a && del %b

…would mean each pixel in the output is .000005 degrees?
How can I calculate what the -tr resolution needs to be to keep the output as the same resolution as the input?


texts:
I converted a .tif file into a tileset using the gdal2tiles plugin on QGIS. The resulting files are a bunch of .png files, .html files and a .xml file. I zipped these and tried uploading them to Carto but got the error that the file types were not supported.

What am I doing wrong? How does one upload these tiles?


texts: Use identify tool on a specific coordinate

texts:
Just a simple question. I would like to use the identify tool for a specific coordinate to get the pixel value. Unfortunately it's not possible to place the cursor that precisely. I used the go to xy tool to set a point but even that doesn't help.

Anybody knows how to solve this?


texts: Tables in ArcMap 10.0 - are they layers or tableviews

texts: Is a table in ArcMap a layer or a tableview.  I want to be able to feed a table into an UpdateCursor via a toolbox script.  If my  table has a definition query on it how do I check for that with python?


texts: Leaflet control window is not a function

texts:
I trying to add a leaflet control window to my map, but i get an error in the console of browser. The error is in console is:

Uncaught ReferenceError: L is not defined
at L.Control.Window.js:1
Uncaught TypeError: L.control.window is not a function
at visor.js:16

I copied the folder of plugin to my public html, and added the url to css and js file of plugin in my index.html.

My code is:

var win = L.control.window(mapa, {title: 'Hello world!', content: 'This is my first control window.'})
    .show();