texts: Contours generated from merged TIFs huge - much bigger than individually

texts:
I'm generating contours from four 5x5 degree SRTM DEMs in southeastern Australia. I'm using gdalwarp to reproject the DEM, then gdal_contour to generate contours.

First, I generated contours for each block separately. The shapefiles were 9, 45, 73 and 270 MB. Then I merged the source TIFFs, and generated a single contour file. It's 2.0GB (that is, more than 4x the size of the combined individual contours).

Why is this? Here are the two scripts:

Individual contours:

for f in srtm_66_19 srtm_65_19 srtm_66_20 srtm_65_20
do
gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3785 -r bilinear $f.tif $f-3785.tif

# contours
gdal_contour -a elev -i 20 $f-3785.tif $f-3785-contour.shp

done

Merged contours:

gdal_merge.py srtm_66_19.tif srtm_65_19.tif srtm_66_20.tif srtm_65_20.tif -o srtm.tif
f=srtm
gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3785 -r bilinear $f.tif $f-3785.tif
gdal_contour -a elev -i 20 $f-3785.tif $f-3785-contour.shp

And the complete file listing:

-rw-rw-r-- 1 ubuntu ubuntu  20M Sep  9 17:36 new/srtm-3785-contour.dbf
-rw-rw-r-- 1 ubuntu ubuntu  399 Sep  9 17:08 new/srtm-3785-contour.prj
-rw-rw-r-- 1 ubuntu ubuntu 2.0G Sep  9 17:36 new/srtm-3785-contour.shp
-rw-rw-r-- 1 ubuntu ubuntu 7.6M Sep  9 17:36 new/srtm-3785-contour.shx
-rw-rw-r-- 1 ubuntu ubuntu 404M Sep  9 17:07 new/srtm-3785.tif
-rw-rw-r-- 1 ubuntu ubuntu 8.9M Sep  9 13:34 srtm_65_19-3785-contour.shp
-rw-rw-r-- 1 ubuntu ubuntu  68M Sep  9 12:28 srtm_65_19-3785.tif
-rw-rw-r-- 1 ubuntu ubuntu  69M Sep  9 12:28 srtm_65_19.tif
-rw-rw-r-- 1 ubuntu ubuntu  73M Sep  9 13:38 srtm_65_20-3785-contour.shp
-rw-rw-r-- 1 ubuntu ubuntu  67M Sep  9 12:28 srtm_65_20-3785.tif
-rw-rw-r-- 1 ubuntu ubuntu  69M Sep  9 12:28 srtm_65_20.tif
-rw-rw-r-- 1 ubuntu ubuntu  45M Sep  9 12:28 srtm_66_19-3785-contour.shp
-rw-rw-r-- 1 ubuntu ubuntu  68M Sep  9 12:28 srtm_66_19-3785.tif
-rw-rw-r-- 1 ubuntu ubuntu  69M Sep  9 12:28 srtm_66_19.tif
-rw-rw-r-- 1 ubuntu ubuntu 270M Sep  9 13:37 srtm_66_20-3785-contour.shp
-rw-rw-r-- 1 ubuntu ubuntu  67M Sep  9 12:28 srtm_66_20-3785.tif
-rw-rw-r-- 1 ubuntu ubuntu  69M Sep  9 12:28 srtm_66_20.tif

Why so big? Is there anything I can do about this? I'd prefer to work with merged files as they're easier to manage in TileMill.


texts: PostGIS nest ST_ClusterKMeans on results of ST_ClusterKMeans

texts:
I am using ST_ClusterKMeans in PostGIS to perform a spatial clustering.
CREATE TABLE kmeans10 AS
WITH s AS
(SELECT p.*, ST_ClusterKMeans(p.geom, 10) OVER() cid
 FROM p)
SELECT ST_UNION(s.geom) OVER (partition by s.cid ORDER BY s.cid) geom, s.cid,
FROM s

The initial result looks like this:

Based on the initial result, I'd like to divide the cluster with the maximum number of points into another 5 sub-clusters.
WITH s AS
(SELECT * FROM kmeans10
WHERE cid in
(SELECT cid
FROM kmeans10
GROUP BY cid
ORDER BY COUNT(cid) DESC
LIMIT 1))

SELECT s.geom, ST_ClusterKMeans(s.geom, 5) OVER () s_cid
FROM s
ORDER BY s_cid DESC

However, I do not get the 5 sub-clusters as desired but this notice:
kmeans_init: there are at least xxx duplicate inputs, number of output clusters may be less than you requested

What's wrong with my query? Is there any better solution to nest the spatial clustering action?
UPDATE ON 2021-05-10
I refine the initial kmeans10 query as per the advice from @geozelot, by using GROUP BY instead of ST_Union. Then the consecutive query works properly.
CREATE TABLE kmeans10 AS
WITH s AS
(SELECT p.*, ST_ClusterKMeans(p.geom, 10) OVER() cid
 FROM p)
SELECT s.geom, s.cid
FROM s
GROUP BY s.cid, s.geom


texts: PyQt Push Button randomly stops working in PyQGIS plugin

texts:
I keep running into a very weird issue with clicked.connect not working in various plugins I am working on. The error happens randomly in QGIS 3.x versions and it seems to happen on some machines only.
It is very hard to research this online and I am not sure if it is linked to QGIS or to PyQt or to a combination thereof.
What happens: I click a push button but the method that has been connected to it just does not trigger (even though it has triggered many, many times before).
The only fix I found so far (after hours of trying the most random things): I intentionally insert an error that does not break the plugin on reload but only once the GUI is opened and the handlers are connected.
Workflow:
My button (btnAddShapefile) stops working (i.e. the method addShapefileToQgis is just not called when clicking the button)
self.btnAddShapefile.clicked.connect(self.addShapefileToQgis)

Now I change the method to a method that does not exist, e.g.:
self.btnAddShapefile.clicked.connect(self.HELLO)

Now I reload the plugin and open it again. The plugin then crashes, because:
AttributeError: 'UploadData' object has no attribute 'HELLO'

(UploadData is the name of the class)
Now I fix the error, i.e. I revert to the method that actually exists:
self.btnAddShapefile.clicked.connect(self.addShapefileToQgis)

Now I reload the plugin, open it again, click the button, and the method is called without any issues.
I am absolutely clueless why this happens in the first place and why an intentional error that breaks the plugin actually helps me fix it.


texts: How to get JSON data from cql_filter WMS layer with OpenLayers 4

texts:
I am trying to get the JSON feature data which I have filtered using cql_filter.

<script>
      var wmsSource = new ol.source.TileWMS({
        url: 'http://localhost:9090/geoserver/ilfs/wms',
        params: {'LAYERS': 'ilfs:dharani_samplevillage',
            'TILED': true, 'format': 'image/png',
            'VERSION': '1.1.1', 'SRS':'EPSG:4326',
                "TRANSPARENT": true,
                "ISBASELAYER": true,
                "RATIO": 1,
                "TILED": true,
                "CQL_FILTER": "bhucode=11223344444555 and kid=243"
                //kid = 237, 236, 78, 241, 243
        },
        serverType: 'geoserver',
        transition: 0
        //crossOrigin: 'anonymous'
      });

      var wmsLayer = new ol.layer.Tile({
          source: wmsSource,
          opacity: 0.7,
          visible: true
      });

      var view = new ol.View({
          projection: 'EPSG:4326',
          //Longiture first, then latitude
          //center: ol.proj.transform([76.63612, 16.94924], 'EPSG:4326', 'EPSG:3857'),
          //center: [76.63612,16.94924],
          center: [0, 0],
          zoom: 18
      });

      var map = new ol.Map({
          layers: [wmsLayer],
          target: 'map',
          view: view,
          renderer: 'canvas'
      });<script>

With the above code displaying the filtered map, but I need to retrieve feature data with the format of JSON, and that should be in load into browser console.


texts: Eliminate irrelevant points based on time and distance

texts:
I'm collecting raw GPS points from android device for tracking, now when i'm creating a linestring based on the points i collected there are some points which makes no sense (due to GPS fluctuations). I have created a query to select only those coordinates which are satisfying a condition i.e (1.4 meter in 1 sec) which is average human walking distance. But still some points are zigzag. I'm not an expert in gis and also PostGIS

Here is the query i used

SELECT a.*
FROM location_history a
JOIN location_history b
  ON (a.id+1 = b.id)
where a.socket_client_id=187 and DATE(a.date_time)='2019-04-17'
  and ST_DistanceSphere(ST_SetSrid(ST_MakePoint(a.lng, a.lat), 4326),
                        ST_SetSrid(ST_MakePoint(b.lng, b.lat), 4326))
      <= ((EXTRACT(EPOCH FROM (b.date_time - a.date_time)))*1.4)
order by a.date_time;

Here are some points which are zigzag


texts: modify raster metadata in qgis to change raster coordinates

texts:
I have two rasters of the same place, however they differ 3metres between them. So I want to change the metadata of one of them, to match the coordinates of the other one. (they both use the same projection and coordinate system, its just a raster geolocation issue)

I've read that this can be done using "gdal_translate -a_ullr X-coordinate Y-coordinate" but the command is failing.

Anyone knows an easy/fast way to change metadata coordinates in QGIS?


texts: How is device-responsive layer selection working for you?

texts:
I ran into problems getting layer selection to work with data categories that I didn't document well in this process, but I would just like to show my results here in case people are interested in adding responsive navigation buttons to their maps. Here is my result: http://jsbin.com/sihayino/4/edit (Ultimately, I'd like to add a custom base map—still figuring out how to keep the delicate balance I have here while adding that icing to my cake).

Has anyone a better method for incorporating layer toggling for viewing map layers on mobile devices/tablets? If not, it's no big deal. I'm pretty ok with this.


texts:
I have two file geodatabases with identical schema that have a feature class and four related tables. I need to copy features from one to the other with the information from all of the related tables.

Does anyone have a work flow for this?

All of the standard ways of moving data (append, simple data loader, copy) do not accomplish what I am trying to do.


texts: Merge adjacent parallel lines and add count as attribute in QGIS?

texts:
I have a layer containing multiple lines representing flown tracks of airplanes. My goal is to create maps to analyse the flow of traffic in certain airspaces.

I would like to merge all parallel, or near parallel line, segments that are within a certain distance of each other (say, 1 nautical mile).

The new, merged, lines should contain an attribute containing the number of original lines that have been merged. For example, if 5 parallel lines are merged to one, the new merged line should have a count attribute of 5, which would allow me to style it accordingly.

Consider the following schematic example:

Input:

Three individual lines, two of which have a common parallel segment.

Desired output:

The common part of the green and orange lines have been merged to a single line (red). The new red line will have a count attribute of 2.

Note: in this example, the red line will be an average of the two original line segments. The remaining green and orange line segments do not have to connect exactly to the new red line. In other words, the topology does not have to be intact, since the result will only be used for visual representation.

Basically, I am trying to create something like this:

(source)

How would I go about achieving a result like this in QGIS? Bear in mind the actual dataset is much larger than my example above, so the solution should be scalable.

I came across this tutorial which does something very similar to what I am trying to do. However, I would like to be able to do it using QGIS.


texts: Making PDF exported from ArcMap have pop-up when clicking on feature?

texts: When converting your map into a PDF from ArcView 10.1, is there a way to make the PDF Map to have a pop-up when you click on a tract?


texts: Python labeling code - remove empty values

texts:
I'm a GIS technician who is (slowly) learning python. I'm using it to generate custom labels in ArcGIS Desktop 10.4. I have a point layer with 10 fields but no point uses all 10. I'm trying to display only the fields that has a value.

I think part of the problem is that my table has  values, "none", or "0" listed in some of the fields. I've tried replacing every "empty" value with a null but I don't think python/Arc likes this.

I've attached pics of what I have and what I want.

def FindLabel ([ELEV_OBS], [INVERT_D_1], [INVERT_D_2], [INVERT_D_3], [INVERT_D_4], [INVERT_D_5], [INVERT_D_6], [INVERT_D_7], [INVERT_D_8], [INVERT_D_9], [INVERT_D_10] ) :
    return ([ELEV_OBS] + "\n" +
            [INVERT_D_1] + [INVERT_D_2] + "\n" +
            [INVERT_D_3] + [INVERT_D_4] + "\n" +
            [INVERT_D_5] + [INVERT_D_6] + "\n" +
            [INVERT_D_7] + [INVERT_D_8] + "\n" +
            [INVERT_D_9] + [INVERT_D_10] )


texts: Creating fontset to use it later as Character Marker Symbol in ArcGIS Desktop?

texts:
Using Inkscape to create fonts as indicated in Creating new custom vector marker symbol for ArcMap? does not works well for me, as the drawn symbol is not aligned correctly in the Symbol Property Editor, and I don't want to set offset for each of my symbols.

What is the best software (free, if possible) to use for creating these fontsets in a correct way?

I have followed the aformentioned tutorial. The problem is, that in Inkscape the image seems centered, but the SVG font editor automatically adds some additional spaces around it (called bearings) and that makes the image to misalign to the center. Have you successfully created a centered glyph somehow?


texts: Check for intersection of raster and vector in R

texts:
I have a rasterstack and a vector-geometry (sf-object). The raster's values are either ǸA or some integer. I would like to check if the raster (in the best case just as a binary response) has any intersection with the shape in the points that are non-NA in the raster. And I don't want to use the bounding box of the raster.
I thought about polygonzing the raster where the values are not NA. But this seems a lot of effort and space. I'm not asking for any specific code, but maybe an approach that could be used in R


texts: Any open source alternative for ESRI maps for office

texts:
I am looking for an alternative for ESRI maps for office. I know it is very good product for analyzing data within MS office environment.

My requirement:

-I have few stores data contains location details.

-Trying to analyzing the trend on the map itself in the excel document +Statistical analysis

-Is there any tool/app in the Open office(or MS Office & going through EuroOffice tool but it looks like its paid extension)?

-Researching on Google.


texts: Create lines from a Blob field in ArcMap/Catalog or FME (or any other way)

texts:
Is there some way to create lines from geometry stored in a Blob field? It's a SDE-database line feature class that has crashed in some way, but I can still get the raw data and the feature class works in a third party application. No, I can't recover the database from a backup, because it turns out our IT department only stores backups 7 days.

The said third party application reads it's data from a view made of the data table and geometry table. This view I have exported to a regular file geodatabase (to be on the safe side with experiments). So, I need to somehow extract the geometry data from the Blob field, so I can recreate the lines. I'm having a hard time finding anything like this on Google.

There's also a field named ENTITY, and one named NUMOFPTS. I'm guessing these have got something to do with the Blob field (named POINTS).


texts: Applying symbology to a layer with python, How to only have symbology in TOC that the sub layer has

texts:
i am trying to apply a symbology to multiply layers that i have created from one feature class.

The original feature class is used to create several sub layers see
layer loop

this involves using arcpy.MakeFeatureLayer_management

Before using arcpy.SaveToLayerFile_management

i would like to apply some symbology to the layers that have been created.  I thought if i took the original feature class and applied some the desired symbology and then saved this as a layer file. This could be used when using the Apply Symbology From Layer

The layers were made and then the symbology was applied.

But when drag the layers into ARC MAP , they all have the identicaly symbology showing in the TOC.

The actually polygons are coloured up correctly.

Example _ Original feature class symbology

1-100  = ORANGE

101 - 200 = RED

201- 300 = BLUE

301 - 400 = GREEN

401 -500 = PURPLE

Created layer  value

10

200

300

In the TOC the symbology contains all of the items within the the template

How the layer is showing in the TOC

1-100  = ORANGE

101 - 200 = RED

201- 300 = BLUE

301 - 400 = GREEN

401 -500 = PURPLE

Desired Output

1-100  = ORANGE

101 - 200 = RED

201- 300 = BLUE

Is it for the TOC to show the desired output. My worry is that someone might look at the symbology and assume there are data points. When actually none exist


texts:
I downloaded the AWS Earth search notebook from the intake-stac examples and whithout changing anything. Just run the cells...
But I've got everytime the error: KeyError: 'open_stac_item_collection'
Here is my code:
%matplotlib inline
import intake
import satsearch

bbox = [35.48, -3.24, 35.58, -3.14] # (min lon, min lat, max lon, max lat)
dates = '2020-07-01/2020-08-15'

URL='https://earth-search.aws.element84.com/v0'
results = satsearch.Search.search(url=URL,
                                  collections=['sentinel-s2-l2a-cogs'], # note collection='sentinel-s2-l2a-cogs' doesn't work
                                  datetime=dates,
                                  bbox=bbox,
                                  sort=['<datetime'])

items = results.items()

catalog = intake.open_stac_item_collection(items)

This is the result/error:
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
~/opt/anaconda3/envs/openspacedata_app/lib/python3.7/site-packages/intake/__init__.py in __getattr__(attr)
     60     try:
---> 61         return gl[attr]
     62     except KeyError:

KeyError: 'open_stac_item_collection'

During handling of the above exception, another exception occurred:

AttributeError                            Traceback (most recent call last)
<ipython-input-29-45f007afa1be> in <module>
----> 1 catalog = intake.open_stac_item_collection(items)

~/opt/anaconda3/envs/openspacedata_app/lib/python3.7/site-packages/intake/__init__.py in __getattr__(attr)
     61         return gl[attr]
     62     except KeyError:
---> 63         raise AttributeError(attr)
     64
     65

AttributeError: open_stac_item_collection

Update 1
After the solution of Kadir and replace intake_stac.catalog.StacItemCollection with intake.open_stac_item_collection, I've got the following error:
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-61-0321aa237297> in <module>
----> 1 catalog = intake_stac.catalog.StacItemCollection(items)

~/opt/anaconda3/envs/openspacedata_app/lib/python3.7/site-packages/intake_stac/catalog.py in __init__(self, stac_obj, **kwargs)
     61             self._stac_obj = self._stac_cls.from_file(stac_obj)
     62         else:
---> 63             raise ValueError('Expected %s instance, got: %s' % (self._stac_cls, type(stac_obj)))
     64
     65         metadata = self._get_metadata(**kwargs.pop('metadata', {}))

ValueError: Expected <class 'pystac.catalog.Catalog'> instance, got: <class 'satstac.itemcollection.ItemCollection'>

Update 2
Now this works. But then, I'll get an error on the next line:
item = catalog['S2A_36MYB_20200804_0_L2A']

Full error:
KeyError                                  Traceback (most recent call last)
<ipython-input-75-60497ae0cb6f> in <module>
----> 1 item = catalog['S2A_36MYB_20200804_0_L2A']

~/opt/anaconda3/envs/openspacedata_app/lib/python3.7/site-packages/intake/catalog/base.py in __getitem__(self, key)
    415                 out = self[part]
    416             return out()
--> 417         raise KeyError(key)
    418
    419     def discover(self):

KeyError: 'S2A_36MYB_20200804_0_L2A'


texts: Source for high resolution satellite images

texts:
I am looking for high resolution imagery of the Nicobar Islands.
Ideally 3 images (2004, 2005 and 2014/15).
My budget is 1500 USD.
Which image providers would you advice?
Until now I have only worked with free (Landsat) imagery therefore I do not know where to ask.


texts: Converting image pixels to a list so that I can use .map() function

texts:
There is an image that is a 2-dimension matrix of the pixel values, and I want to do some operations for each pixel.

// For each pixel i (i = 1,..., N) in the image, the pixel value is P_i;
// For the other pixels j (j != i) in the image, P_j = P_j - P_i;
// For a given value q, if (P_j < q) P_j = 1; else P_j = 0;
// Finally, P_i = sum(P_j);

//The code provided by JonasV is a version aimed at only one pixel i.
//However, what I want is that every pixel in the image (total number is N) do the same calculation.
//Therefore, it seems need a list to contain these pixels firstly.
//because the value of pixel i is a variable, not a constant.

var i = ee.Image(5)  //i indicates the value of one pixel i in the image
var j = ee.Image(2)  //j indicates the other pixels j (j!=i) in the image

var j_2 = j.subtract(i)

var q = ee.Number(3)

var j_3 = j_2.lt(q)

Map.addLayer(j_3)

Therefore, I think I need to arrange these pixels into a list, in order to use .map() function.

How can I achieve this goal?


texts: Avoiding or reducing invalid geometries when geoprocessing in ArcGIS Desktop?

texts:
I use geoprocessing with large polygon datasets.

I use Multipart_to_Singelpart and Repair_Geometry to remove all invalid geometries, prior using the tools like Union or Intersect.
Regardless the correct geometries in the input data the result will contain invalid geometries. So, the invalid geometries appear during geoprocessing with ArcGIS.

These invalid geometries hamper further processing (ArcGIS freezes), so I desperately put a Multipart_to_Singelpart and Repair_Geometry in my script after every step of geoprocessing.

But is there a best practise or setting to avoid or reduce invalid geometries when geoprocessing in ArcGIS?

Using a projected CS (EPSG:4647) with ArcGIS 10.2.2 and ArcGIS 10.6


texts: R - Easy way to reproject an extent object

texts:
I have an extent object that I would like to re-project to lat and lon.

>outer_extent
    outer_extent
    class       : Extent
    xmin        : -5559753
    xmax        : -4447753
    ymin        : -4447852
    ymax        : -3335852


texts: Model only works one time and then fails, Error 000732, "Precondition fails"

texts:
I have a huge ol' model that's built to copy a bunch of data from one GDB into a different GDB and then delete two of the feature classes. If I build the model from scratch, it'll run fine the first time, but then encounters an error every subsequent time. The only way around this error is to manually re-input the data to be deleted every time. After the first time the model is run, the next time it will technically run but it skips the delete feature class step, even if it's a precondition for other steps. If I open the step, there's a

"Error:
"000732: <value>: Dataset <value> does not exist or is not supported."

notice next to the data input. The data DOES exist, and I can manually re-select it from the folder.

Does anyone have an idea what could be causing this issue? I need the entire model to run automatically without user input, so I need it to be able to run every time without re-inputting the data to be deleted.

I've tried turning off background processing, I've tried running the model with or without the edit window open.

EDIT:
I decided to see if I could run literally just the delete as a simple, separate model and it has the same issue. I've boiled it down to just this:

It will run but nothing will happen, with this message:

Executing (Delete Media Filter SMPs): Delete O:\GIS\ArcGIS_Pro\Planning\workspace.gdb\gswiMediaFilter #
Start Time: Thu Feb 28 16:42:20 2019
The process did not execute because the precondition is false.
Succeeded at Thu Feb 28 16:42:20 2019 (Elapsed Time: 0.00 seconds)


texts: How to get raster layer metadata in python for QGIS 3.6?

texts:
I am loading imagery as raster layer through a URL. I need to store its metadata as JSON string in db for future reference.

I tried getting the metadata via layer.metadata() and then tried to get json.dump() for that object but an empty string was returned.

When I use the layer.htmlMetadata() function I do get all the information but I need all the details as JSON.


texts:
So I've been assigned to clean up someone else's model.

Here's the situation.  The model is supposed to iterate over a few hundred folders containing feature classes and convert the contents into shapefiles.  The model fails if a file called "arc" does not exist in the folder, so I need to write some python to check if it is there, and if it does not, then skip the current iteration and continue onto the next folder.

My problem is that with the way modelbuilder is set up, there is no way to skip an iteration within an iterator tool.  Any ideas?  I'm halfway thinking of just rewriting the whole tool in python.

Here's a screenshot of the iterator

It has been suggested that I use calculate value to skip an iteration.  The code should be fairly simple, however how do I use the true/false output to skip an iteration?


texts: Getting MongoDB working with GeoServer 2

texts:
I'm trying to use MongoDB with GeoServer but it is actually not working at all.
Someone once posted this: Getting MongoDB Working with GeoServer, but since he could solve his problem without any help there's not much info about how I can do it too.

I got the MongoDB GeoTools plugin (what is the correct .jar file and on which folder should I put it?). I'm trying to make the "new data source" option appear like this guy did.

If someone could help me in a "tutorial like" way I'ld be very happy :)


texts: Problem with batch mode Saga Raster Calculator

texts:
I have to do a couple of raster multiplications for which the batch mode would save me a lot of time. But the only thing I get are rasters that contain nothing but nan values. If I multiply the rasters one by one without batch mode, it works well. the batch mode also works with the same rasters doing other operations like additions or the gt() command. The Grids product tool shows the same behavior.

Sometimes I also get the alert Oooops! The following output layers could not be open which disappears after restarting QGIS.

To me it seems like that is a bug, is there another easy way to do those operations automatically? Let me know if you need the text from the Log File.

Edit: I know there is another Raster Calculator in GRASS, but to use this one i'd have to import all layers to GRASS which would take me a lot of time as well I fear.


texts:
I am trying to calculate the mean raster from a large RasterStack, with the following (using package raster):

raster_mean <- calc(raster_stack, conditonalmean)

However, I do not want a standard mean, but I wish to add some conditions to the computation. I have written the following function (where three conditions appear, one is to remove null values and the others are numerical):

conditonalmean <- function(x, na.rm){
v = var(x,na.rm = na.rm)
mean(v[v>=-10 & v<=10])

}

Which however returns the following error:

Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) :

cannot use this function

Any suggestion?


texts: Using raster2pgsql in Python console of QGIS?

texts:
I have to upload GeoTiFF files into a PostGIS/PostgreSQL database and I am new to raster2pgsql.

How do I run a raster2pgsql command in the Python console?

import psycopg2
import subprocess

db_name = 'enter_qgiscloud_db'
db_host = 'db.qgiscloud.com'
db_user = 'enter_qgiscloud_user'
db_password = 'enter_qgiscloud_pw'

conn = psycopg2.connect( "dbname={0} host={1} user={2} password={3}".format( db_name, db_host, db_user, db_password ) )
cursor = conn.cursor()

cmds = 'raster2pgsql -s 3857 -p -Y -I -C -M C:\qgis_cloud_data\ -F -t auto .tif | psql'
subprocess.call(cmds, shell=True)

cursor.execute(cmds)
conn.commit()


texts: How to calculate distance in meter between QgsPoint objects in PyQgis?

texts:
I tried to calculate distance between QgsPoint objects in qgis, i used QgsPoint StartPoint.sqrDist(EndPoint).

but i don't find an expected result, i also tried

import math
math.sqrt(StartPoint.sqrDist(EndPoint))

but still getting a very high result.
does somebody know a way how to calculate distance in meter between 2 Qgspoint objects.


texts: Where are default layout templates stored for ArcGIS Pro?

texts:
Previously in ArcMap 10x you could save custom layout MXDs to C:\Program Files (x86)\ArcGIS\Desktop10.3\MapTemplates or C:\...\AppData\Roaming\ESRI\Desktop10.3\ArcMap\Templates. However I can't seem to find a similar directory for ArcGIS Pro, does anyone know if there is a directory for these layouts I think they are .pagx files?

I looked in C:\Program Files\ArcGIS\Pro\Resources and C:\...\AppData\Roaming\ESRI\ArcGISPro but couldn't find anything.

My goal is to have custom layouts show in the 'New Layout' window (shown below). I guess its possible Pro creates these on the fly when you click. I realize I can just use the 'Import a layout file' option but it would be cool to have them in the window. Any ideas?


texts: Why is the unique ID field in Distance matrix empty for some inputs?

texts:
I have some coordinates in a .TSV file I import using the delimited text layer button. I have a different set of coordinates (actually, one point) that I import in the same way.

What I've tried so far is the following:

I upload the coordinates (test, in pink) using the delimited text layer button.
I upload the localized coordinate using the same method above (test_Sheet(1)).
I wanted to try to accomplish this by going vector-->Analysis Tools-->Distance Matrix.

However for one of those, the one with multiple coordinates, it's not giving me any scroll box for to input the unique ID field. If I switch things around so it's the localized coordinate, it's totally fine no matter what.

Because of this I can't run analysis.
Can anyone figure out what my issue is?


texts: Choosing right scale for global aggregations in GEE

texts:
What is a correct scale to aggregate global datasets like VIIRS Nighttime Lights in GEE, a scale that will retain maximum (original) amount of information when aggregating over any place in the world (any latitude)?
For example, VIIRS Nighttime Day/Night Band Composites have resolution of 15 arc seconds, which corresponds to ~500 m at the Equator but up to ~320 m in high latitudes. Does it mean that setting scale: 500 will keep original amount of information everywhere?
result = VIIRScoll.reduceRegions({
  collection: areasInAnyLatitudes,
  reducer: ee.Reducer.mean(),
  scale: 500, // ~15 arc sec on equator but ~9 arc sec in 49 degrees north latitude.
});


texts: How do I get my Data Driven Pages Inset Map to center the grid featured?

texts: I'm doing data driven pages manually. I had this working out just fine but I must have changed the settings and can't figure out how to get the locator map to follow the grid featured.


texts: Extract one geometry type from GML file in GeoPandas

texts:
I try to use GML data from Landesvermessung Brandenburg in GeoPandas. It generally works to load files like this:
df = gpd.read_file(r'/vsizip//vsicurl/https://data.geobasis-bb.de/geobasis/daten/3d_gebaeude/lod2_gml/lod2_33379-5805.zip/lod2_33379-5805_geb.gml')
This gives the following warning, which I tend to ignore:
C:\Daten\Anaconda\lib\site-packages\geopandas\geodataframe.py:422: RuntimeWarning: Sequential read of iterator was interrupted. Resetting iterator. This can negatively impact the performance.
  for feature in features_lst:

But one file in the bulk does contain one odd polygon, where linestrings are expected:
df = gpd.read_file(r'/vsizip//vsicurl/https://data.geobasis-bb.de/geobasis/daten/3d_gebaeude/lod2_gml/lod2_33379-5806.zip/lod2_33379-5806_geb.gml')
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Daten\Anaconda\lib\site-packages\geopandas\io\file.py", line 133, in _read_file
    columns = list(features.schema["properties"])
  File "C:\Daten\Anaconda\lib\site-packages\fiona\collection.py", line 208, in schema
    self._schema = self.session.get_schema()
  File "fiona/ogrext.pyx", line 719, in fiona.ogrext.Session.get_schema
  File "fiona/_geometry.pyx", line 81, in fiona._geometry.normalize_geometry_type_code
fiona.errors.UnsupportedGeometryTypeError: 2147483648

I only want linestrings (MULTILINESTRING Z) and want to ignore other geometry types. Is there a way in GeoPandas?


texts:
Using ArcGIS Desktop 10.1, I have 2 Feature Datasets in a geodatabase. One dataset contains bus routes and the other contains bus stops for each route. I want to create a script that will snap the bus stops to the vertex of the relevant bus route. There are over 100 routes and stops. They are stored in each dataset in order. They are named like this S102, S104 etc. for stops and B102, B104 for routes. I'm unsure about how to do a for loop that will get a stop and then get the right bus route

I've started writing a script. What I'm trying to do is get the last 3 characters of the stops ID (for S102 that would be 102) and then search for B + those characters in the routes dataset. I haven't gotten to the searching bit yet but I'm getting a Parsing error SyntaxError: invalid syntax (line 20) which is arcpy.env.workspace = RoutesD

import arcpy

StopsD = "C:\College Stuff\Dissertation\POWSCAR2011\Practice.gdb\Stops_temp"
RoutesD = "C:\College Stuff\Dissertation\POWSCAR2011\Practice.gdb\Routes_temp"

arcpy.env.workspace = StopsD
arcpy.env.overwriteOutput = True
StopsList = arcpy.ListFeatureClasses()

for Stops in StopsList:
        Stopname = Stops[1:4]
    Stopname2 = "B" + Stopname
        StopnameEx2 = "\"" + Stopname2 + "\""

    Routes = str(RoutesD) + "//" + "B"+ Stopname

    arcpy.env.workspace = RoutesD
    RoutesList = arcpy.ListFeatureClasses(Stopname2)
    for routes in RoutesList
        print routes

arcpy.env.workspace = StopsD


texts:
I have created an sde table in PostgreSQL 9.5 named parcels from ArcCatalog version 9.5.

ID   NAME  TYPE_ID
1    x     1
2    y     1
3    z     2

I created parcel types as Domain in ArcCatalog. So I can see my data in ArcMap like this:

ID   NAME  TYPE_ID
1    x     Type-1
2    y     Type-1
3    z     Type-2

So I will publish this table on ArcGIS Server Enterprise 10.5. How will users see my TYPE of parcels, when they get feature info. As numeric or domain description value?


texts:
Is there a WFS for OSM data?
Is there a WFS for OSM data?


texts: Issue with 'ESRI.ArcGIS.Client.Geometry.Geometry' in WCF service

texts:
I am having this issue with using Geomtry that is returned from a WCF service. As part of the service, I have added
a service reference for the map service that I am using in this project. I am trying to return 'ESRI.ArcGIS.Client.Geometry.Geometry as a return type but when I assign the Geomtry to my result, it comes back with this error: 'Geometry' is an ambiguous reference between 'GIS.Services.Data.BoundaryMapService.Geometry' and 'ESRI.ArcGIS.Client.Geometry.Geometry'. I am not sure what needs to be done to resolve this issue. Appreciate any feedback. Thanks Jay

ESRI.ArcGIS.Client.Geometry.Geometry resultGeometry = null;

                for (int i = 0; i < attributes.Values.Count(); i++)
                    if (attributes.Values[i] is PolygonN)
                    {
                        PolygonN featureExtent = attributes.Values[i] as PolygonN;
                        resultGeometry = featureExtent; // Error here...

                    }


texts:
After installation of the Midvatten plug in (QGIS 2.12 Lyon), I tried to get confidence with the plugin following the tutorial. When trying to plot stratigraphy data, a QGIS error shows up: "the stratigraphy plot failed, check Midvatten settings and your data". I'm going to link an image of the error

Is there anyone having some experience with this plug in?


texts: Reverse geocoding multiple Lat\Long points in single go with QGIS?

texts:
I have a file which is having Lat\Longs for thousand of records, I need address for all these records to make it useful, I got a way to reverse geocode these records one by one but this process is time taking.

Can anyone suggest, how to reverse geocode multiple records in QGIS?


texts:
I have a found a lot of questions about converting shapefiles to GeoJSON, but how do you convert a GeoJSON to a Shapefile?
I've seen  Converting GeoJSON to Shapefile using ogr2ogr? that uses ogr2ogr, but they look like they are commands for a shell.
Is there a pure Python script to do it?
I've also found this solution that runs without any errors, but the output is still a GeoJSON.
Perhaps explaining how to write as a .shp?
import geojson
import subprocess
import urllib.request as ur

url = 'http://ig3is.grid.unep.ch/istsos/wa/istsos/services/ghg/procedures/operations/geojson?epsg=3857'
response = ur.urlopen(url)
data = geojson.loads(response.read())

with open('D:/Scripts/Stand19North.geojson', 'w') as f:
    geojson.dump(data, f)

args = ['ogr2ogr', '-f', 'ESRI Shapefile', 'D:/Scripts/converted_shp.shp',
        'Stand19North.geojson']
subprocess.Popen(args)


texts: How to interpret the red trend map of time series on Google Earth Engine?

texts:
I'm starting to play with Google Earth Engine. Following https://docs.google.com/document/d/1mNIRB90jwLuASO1JYas1kuOXCLbOoy1Z4NlV1qIXM10/edit# I managet to put together this https://code.earthengine.google.com/625b11fca42309bd7f37159434cd21cf, which generates a time series plot of NDVI. I kind of understand it until the plot.

After that it sets independents and dependent variables. And computes the trend which is added to the map.

May I ask for help on how to interpret this? What red means?


texts:
I have a bunch of geolocated points that I am looking to cluster on a google map using hexagons. These clusters will be stored in a db for quick retrieval and will be reduced to a single point and radius. I am using the http://www.geocodezip.com/v3_MW_example_eshapes.html library to draw the hexagons but because of the mercator projection the hexs overlap at higher latitudes.

Here is a jsfiddle attempt with the overlap http://jsfiddle.net/kworthin/Lkuqykog/2


texts: Follow-up query on "QGIS expression for label text to display vertex coordinates of polygons"

texts: Just a follow-up query about the existing thread QGIS expression for label text to display vertex coordinates of polygons. Is there an expression in QGIS to collect all the coordinates from the geometry generator and add to a virtual field of the polygon?


texts: Select interaction shown visually outside layer extent

texts:
I'm running with openlayers 4.3.1

When I add ol.interaction.Select to ol.Map, the most recent Select event is shown visually outside the layer extent

An example where a Select event is triggered when drawing a Point, and the location is shown outside the extent (but not drawn):

https://codepen.io/anon/pen/ayGqdN?editors=1010

My problem is that I cannot figure out how to avoid it to be shown.


texts: Counting number of overlayed points in line feature using ArcGIS Desktop?

texts:
I created a polyline representing a road using GPS points.

Then I have another set of data of random GPS co-ordinates which i have to overlay on that polyline (Road).

I am trying to calculate how many random points occur in each 1 km section of that line(road).

Can anybody give me a solution to this problem?


texts: Change vertical alignment of label in ArcGIS (currently option in grayed out)

texts:
I am trying to place my labels below my polygons using ArcGIS 10.0 (and centered). I have tried using

Properties... Labels

Choosing Symbol... from the interface above

On the General tab, I want to modify Vertical alignment  but options are grayed out and I can't modify anything! Even though it says "bottom", it shows up as centered.

UPDATE: as suggested by the comment, I tried to play around with the placement properties and modified Label weightand Feature weight but it just makes for really odd placement that changes depending on zoom and each feature.

Facts about my settings:

shapefile is a polygon
I have tried with and without a halo
I am using Arial font (and have tried others)
... run out of ideas!!!


texts: Cartesian world sized projection WKT/EPSG code?

texts:
I am after the EPSG code or Well Known Text (WKT) of a projection that puts the centre of the earth a 0,0,0 and each point is a cartesian XYZ distance from that.

My existing data is in WGS84 lat/long.

Does such a projection already exist? If not, any advice on do this projection? (I am very new to GIS)


texts:
I'm not very skilled within the GIS world, however, I would like to know if it's possible to use this map:

http://wfs-kbhkort.kk.dk/k101/wms?service=WMS&version=1.1.0&request=GetMap&layers=k101:theme-startkort&styles=&bbox=12.451709828167994,55.63150896601625,12.64714535243575,55.72503797474024&width=512&height=495&srs=EPSG:4326&format=application/openlayers

within leaflet?

I've looked at the API for leaflet and can see that it is possible to use custom maps/tiles with the tileLayer option. However, I cannot get it to work.

Can someone point me in the right direction?


texts: Anaconda (Python36) error after installing ArcGIS: No module named 'arcpy'

texts:
I have a Anaconda with Python 3.6, and I installed the ArcGIS with the command

conda install -c esri arcgis

but when I went to test it presented the error:

  No module named 'arcpy'.

I saw the tutorials and I didn't find the file Desktop10.3.pth in the path site-packages had a file called Desktop10.3.pth. I didn't find in anywhere.


texts: Create a folium map with a pop up window appearing when clicking a shapefile feature

texts:
I'm trying to create a folium map that will have the following features:

A layer will be created from a shapefile including multiple polygon features.
When clicking on a single feature a pop up window will appear that will either show a link directing to a locally saved image file or will view the image. The image paths are saved in a specific column of the shapefile (one image per feature).

How is this possible (if it is even possible at all)?


texts: Filtering earth engine image collection by image ID

texts:
I'm trying to filter an image collection using not bands within each image of the collection, but by pattern matching in the image ID. Using:

var dataset_2099 = ee.ImageCollection('NASA/NEX-DCP30')
                  .filter(ee.Filter.date('2099-01-01', '2099-12-31'));

we can create an image collection with forecasted climate data. If I, suppose, want to access just the RCP 2.5 data, I'm in a bind because there doesn't seem to be a clear-cut way to filter the image collection based on individual image names, which begin with "rcpXX" where XX may be 26, 45, 60, or 85 depending on the RCP scenario being mapped.

Image ID's can be printed using the following code, where you can see the pattern that I'm describing:

// Image ID listing code by Tyler Erickson
// https://gis.stackexchange.com/a/281489/67264
// Create a list of image objects.
var imageList = dataset_2099.toList(5000);
print('imageList', imageList);
// Extract the ID of each image object.
var id_list = imageList.map(function(item) {
  return ee.Image(item).id();
});
print('id_list', id_list);

The problem is I don't just want to get the names of the images, I want to subset the image collection based on those names.


texts: Leaflet reads GeoJSON [x,y] as [y,x]. How can I correct this?

texts:
Hey I am in the middle of a project and I suddenly realized that my layer is projected somewhere else on the globe. After some experimentation I've realized that this is because the program translates lng for lat. I really don't understand why is that, but at the same time I need to solve the issue ASAP, as time is running out.

Case: [37.994, 24.75906] is being read as [24.75906, 37.994]

I have used the e.latlng method to acquire the data so it really puzzles me that the GeoJSON reads coordinates in reverse.

Edit: OK it turns out that this is a known issue of conflicting standards. Now the issue remains: How can I quickly resolve this since I already have hundreds of JSON [y,x] couples? Last resort for me will be regex but I would like to avoid it.

Do you have any ideas on how to solve quickly the issue? I have hundreds of data like that.

Note: if you know this to be a duplicate please leave a comment and I'll delete the question.


texts: Managing ID's between parent and child forms/tables in QGIS

texts:
I would like to have a line ("ll")and a point-layer ("pl") in QGIS, in this case for example each of them representing the point- or line-geometry of maintenance tasks in the field.
The third layer ("ml") in the project is a geometryless table, which contains all maintenance-tasks (autoincrement id) regardless their geometry-type.
Each "pl"/"ll"-layer has a foreign-key-field for the id of "ml"-table. All are PostgreSQL tables.
Till now the tables were filled with data independently from each other.
I'd like to collect the data using parent/child-forms in the future to avoid renundant data.
To achive this I set up a relation between geometry tables and geometryless-maintenance-tasks-table.
Each of the geometry-layers ("ll"/"pl") is associated to the same "ml"-table.
As a result I got a parent / child form, which shows all existing maintenance-tasks as expected for all former collected data.
If I digitize a new "ll"/"pl"-parent-feature and fill data in the "ml"—child-subform the connection between child and parent is lost for the new data.
Now it is possible to read parent-values (i.e. id's)(https://github.com/qgis/QGIS/pull/34827 - from QGIS 3.14. on) and store them in my child-table. In my case this does not help because I have two parent-layers, which will have the same ID in both tables. So referencing to the layer-id would give wrong results.
So I have to work with foreign-keys ("ml"-id) in the parent-table. I have to write the ID of the child-table-row as a foreign-key to my parent-table. When digitizing the maintenance-tasks I can not know the correct ID's because the child-id is only shown as nextval('….':regclass) before saving the data to database.
Does anybody know how to solve the problem or is there a different approach ?


texts: Mapping from one plane on the other plane despite of masking regions

texts:
have a set of data given here, where in the first and second columns there are the sky coordinates (ra,dec) e.g. right ascension and declination, respectively and in the third and forth, the coordinates in a Cartesian system (x,y), e.g. positions on the CCD image.

I need to make a two-dimensional interpolation surface using coordinates x and y and another using Ra and Dec. The problem is the existence of masked regions, as shown in the figure above. That is what I so far tried and didn't give the right answer:

from scipy.interpolate import griddata
import numpy as np
import matplotlib.pyplot as plt

data = np.loadtxt('test.asc')

ra = data[:,0]
dec = data[:,1]
Xpos = data[:,2]
Ypos = data[:,3]

xi = np.linspace(Xpos.min(), Xpos.max(), 1000)
yi = np.linspace(Ypos.min(), Ypos.max(), 1000)
xi, yi = np.meshgrid(xi, yi, copy=False)
ra_int = griddata(data[:,2:4], ra, (xi.flatten(), yi.flatten()),
                  method='cubic')
dec_int = griddata(data[:,2:4], dec, (xi.flatten(), yi.flatten()),
                   method='cubic')

Using griddata fails and return just NaN values. Is there any way to do this interpolation in order to estimate the values of Ra and Dec from a given x and y coordinates even in the masked regions (map from x and y to ra and dec)?


texts: How to find nearest nodes of polyline in QGIS using Python

texts:
I have two line layers assume A and B and each layer contains multiple features.
I have selected one feature in layer A, by using that I need to find nearest nodes of a features in layer B .

I tried this

for rfeat in LayerA.selectedFeatures():
    rgeom = rfeat.geometry()
    geom_buffer = rgeom.buffer(50, -1)
    for feats in LayerB.getFeatures():
        sgeom = feats.geometry()
        if sgeom.intersects(geom_buffer ) :
            print(feats["id"])

but I am getting all features in layerB.


texts: GRASS GIS i.landsat.toar hangs & holds CPU at 100%

texts:
I use python scripting to automate processing of several Landsat 8 scenes.

I use a VPS (4 core, 8GB memory, SSD) running Ubuntu 14.04 with GRASS 7.0.1

Line in script for i.landsat.toar:

grass.run_command('i.landsat.toar', input='b', output='rb', metfile='/home/ben/map/working/%s_MTL.txt' % scene, sensor='oli', method='dos4', pixel='1000')

Usually this works flawlessly but on occasion (1 in 100) when calculating dark pixel the script will just hang endlessly in some sort of loop.

I added a timer function so that if it hangs for more than 15 min it will raise error and move onto the next scene. That works for now but even when I restart the server and download the Landsat scenes again I still get the same error. I also run GRASS 7 on a Windows 2008 server but do not have the same problem with exact same scene.

Curiously even when it works it still leaves the command open. When opening htop I often have 3 or 4 i.landsat.toar commands still running taking 100% CPU of each core.

Is this a known problem? I had a quick search of known bugs and couldnt find anything.

Update: still no solution to problem. The same scenes keep causing this problem.


texts: Latitude above and below the surface of the Earth

texts:
I know that Geodetic Latitude is measured with respect to the normal at a point on the surface of the reference ellipsoid. But what about points above and below the surface? Do they follow a hyperbolic path? (See the graphic I created.) Or do they follow a straight line?

Wikipedia says, "[Ellipsoidal] coordinates are the natural choice in models of the gravity field for a uniform distribution of mass bounded by the reference ellipsoid."

Latitude should follow gravity if possible, shouldn't it?


texts: Plotting 1M+ intensity reading as custom shapes with Google Maps

texts:
I currently have a database (Mongo, atm) holding roughly 1M points of data corresponding lat/lon/intensity values that I'd like to plot over a region using the Google Maps API.

Ideally, I'd like to have full control over the representation of the points, as the intent is to allow the end-user to control color/opacity/shape through a series of controls. I also need the ability to specify a color range which maps back to the intensity of a given point.

That said, I've investigated using Google Fusion Tables, but for a heavily customized app I found that display styles limiting. I've also looked a bit into marker clustering, as well as a number of external plugins, but none seem to provide the option to generate the marker type dynamically.

Any help/advice/direction is greatly appreciated.

Thanks!


texts: Convert shapefile to KML: KML polygons don't match shapefile

texts:
I've edited my shapefile quite a bit using the Cut and Merge tool (through the Editor toolbox, not the geoprocessing tool). It looks like this:

For some reason, when I export the shapefile to a KML file, it seems to "remember" the old polygons, as though the merge didn't take. The outlines should look identical to the above picture, but they don't:

Is this something to be concerned about? None of my other shapefiles are behaving this way. Before I created the KML file I made sure to save my edits and clear my selected features. I also ran the Repair Geometry tool, and validated my topology (which uses the Must Not Have Gaps and Must Not Overlap rules).

I use ArcGIS for Desktop 10.6.1 on Windows 7.


texts: Can't find the option to clip or union with only selected features in QGIS (2.18.11)

texts:
I just downloaded the new version of QGIS (2.18.11) and I can't find the option to clip or union with only selected features

Typically, there was a selection where you could choose only the selected layers.


texts: I use ArcGIS 10.4. and Win 8.0. .I have a problem with Agent Analyst in a model-level action. I did all steps of the Exercise 2f in the book.The model start running, the Repast output suddenly stopped working, and an error shows " refresh has stopped working", "check online for the solution...".The updateDisplay action in Agent Analyst does not update my ArcMap display. When I deleted the agent's updateDisplay, I didn't receive any Error.The attached file is the print screen of the Error.


texts:
I am writing a application which allow users to save a address into database (db) and in the next step it shows all addresses on a map. But in the saving step I need automatically calculate a GPS coordinates and save them into db too. But I don't know how calculate them because I don't want to show map in saving step so I can't use a Google Maps API.

Note: the Geocoding API may only be used in conjunction with a Google map;
geocoding results without displaying them on a map is prohibited.

Is there a similar service without the restriction?


texts: Adding Map Service Geoprocessing Widget Results Web App

texts:
I have a Web Appbuilder Application using several geoprocessing widgets that have extremely varied result values depending on what the user inputs. Therefore, I would like to symbolize the results with a Map Service, since that is flexible enough to apply unique symbology to the varying range of values.

However, I would like the user to be able to run the tools multiple times without the results getting replaced, as well as automatically enabling the pop-up, as you can do when you are viewing the results with AGOL symblology for the output. The Map service results do not give the option to add the layer to map, does not have the functionality to automatically enable the pop-up, and does not give the option to export results.

I feel stuck between choosing accurate/flexible symbology by viewing the results with a Map Service, and being able to configure a pop-up, download, and add multiple results to the web app session without it. I have looked at the below link, but have not been able to get my application to work.

http://server.arcgis.com/en/server/latest/publish-services/windows/using-result-map-service-in-web-applications.htm

Do I put the code inside each geoprocessing widget? at the bottom? how can I keep multiple Map Service results at once?


texts: QGIS polygonize issue: GDAL DLL load failed

texts: I am using QGIS Desktop 3.10.0 and I have some issues with the polygonize tool. When I run the polygonize tool, a warning shows up saying that there is an “Import Error: DLL load failed” (see Figure 2). I also tried to use this tool in QGIS 3.4 (stable version), but it also did not work. How can I fix this issue, please?


texts:
It is possible to disable the Geosearch Esri controls in leaflet after the first geocode search? When I run Chrome console and I select the search geocode icon this is the code:

<input class="geocoder-control-input leaflet-bar">

I have tried to disable it using

$('.geocoder-control-input leaflet-bar').attr('disabled', true);

but it doesn't work.


texts: Projection for Calculating Distances in Continental United States?

texts:
TL;DR Question: I want to calculate distance from point to edge of nearest polygon for 2,000 points in Continental United States. What is the best projection to use?

Full Question: I have two datasets. First, I have long/lat data for the point of origin of ~2,000 wildfires in the United States. Second, I have a shapefile for Census Designated Places (created by merging shapefiles for individual states from the US Census' website: https://www.census.gov/geo/maps-data/data/cbf/cbf_place.html).

I want to calculate the distance from each of these fires to the edge of the nearest Census Designated Place. Typically, I believe that I would do this using "Generate Near Table" in the ArcToolbox. However, the Census Designated Places shapefile uses a geographic coordinate system (GCS_North_American_1983). So I think I need to change the coordinate system to one that allows me to correctly estimate distance. Unfortunately, I don't know which is the best projection to use.

Would USA Contiguous Equidistant Conic be a safe choice for calculating distances for so many points spread across such a large area?


texts: Random points inside polygon QGIS: what is the unit of minimum distance?

texts: What is the unit of "minimum distance" in the "Random points inside polygon" function in QGIS?


texts: Importing large OSM data to QGis failed (using QuickOSM, SpatialLite layer or add vector layer)

texts:
I downloaded quite large OpenStreetMap dataset from https://www.geofabrik.de/ in osm (or pbf) format. Namely Czech Republic. Geofabrik provides ESRI SHP files too, but I unfortunatelly need more attributes (are not included in prebuilded shp files).

Question is, how to import OSM to the QGIS?
I know 3 ways, but each of them has issues or doesn't work:

1) Vector->OpenStreetMap->Import Topology from XML...->Export Topology to SpatialLite

This works, but:

relations are not imported; some attributes missing (when compared
with OSM map on the web)
polygon holes are not imported correctly (courtyards of houses are
filled)

2) Add Vector Layer...

Simple, but it imports nothing! After openning .osm or .pbf file (using UTF-8 encoding), QGis shows dialog, where number of elements is ZERO (points, lines, polygons, ...)

3) Import OSM File using QuickOSM plugin

Plugin is computing many hours, but at the end, there is no new layer, or sometimes layers named OsmFile, but empty.

I have newest QGis (2.18.12) and newest QuickOSM plugin (sorry, I don't see it's version in gui). My PC is Intel Core i7, 16GB RAM, SSD disc, Windows 10 Pro.


texts: How can I correct a .pgw file for a PNG image?

texts:
I have a PNG file and I need its .pgw file to be set with NAD 83 as the projection.  I have created a .pgw file:

799
0
0
-799
-122.562
37.1946

However when I open the PNG in QGIS or ArcMap 10.2 it is not being placed in the right location.  I have looked around on how to form a .pgw file and I think I am doing it right, but something is obviously missing.  My PNG file has an associated KML which is where I got the coordinates for my .pgw file.  My PNG and.pwg files have the same file name as well.


texts: Interpolating Z values from XYZ point dataset to alternative XY coordinates in QGIS

texts:
I have a dataset with X1,Y1,Z1 values (x coordinate, y coordinate, elevation -> the red dots in the image below). I also have another dataset with alternative X2,Y2 coordinates (-> the blue dots in the image below), but these data do not include elevations Z2.

How can I interpolate Z1 values to X2,Y2 coordinates, for example with IDW interpolation in QGIS?

I would like them to remain as point coordinates and not transform them to a raster format.


texts: QGIS 'Processing_xxx' folder in C:/User/AppData/Local/Temp directory isn't deleted on application exit

texts:
I wrote a Processing algorithm as Python script.
The script chains many processing.run() calls; the outputs of most of them are written to disk using the QgsProcessing.TEMPORARY_OUTPUT syntax for the 'OUTPUT' value of the argument dict.
As a result QGIS saves these all temporary outputs to the 'processing_xxx' folder created in "C:/User/AppData/Local/Temp", each one in a new separate folder. The process is quite long and the final processing_ folder can be around 10Go.
This 'processing_xxx' folder is supposed to be temporary and discarded when QGIS is closed, except it isn't.
I have read many threads on this subject and it seems the issue was signaled (and apprently addressed?), here's the most relevant link I could find so far : https://github.com/qgis/QGIS/issues/33641, several gis.stackexchange entries point to it. It says the issue was fixed in the 3.12 version with Windows 10.
I tried the 3.14 and 3.16 releases and the folder still isn't discarded on exit on my computer; I have windows 10 installed.
I can manually delete the files by going to the 'processing_xxx' folder after QGIS was closed.
Am I doing something wrong? Is anybody else facing the same issue? Would cleaning these temporary files by adding a few lines at the end of the script be a viable solution?


texts: Creating a .bat file to run multiple FME scripts

texts:
I have the task of creating a batch file that will run FME scripts that update data layers to AGOL. I am completely unfamiliar with .bat and have not found anything online that really summarizes the commands to do this.

Does anyone here have experience with creating a bat file that will run FME scripts in succession of each other?

I would like the file to run whether or not the FME files ran successfully or not, however it would be a bonus if I could somehow see if any of the FME scripts failed.


texts: Converting KML To Layer using ArcPy losing objects?

texts:
I am trying to a convert kml file to a layer using python.
This is my code:

import arcpy
arcpy.KMLToLayer_conversion('C:\kml\\result.kml', 'C:\kml', 'result')

But the layer has lost a lot of objects as compared to the layer converted by Conversion tools -> kmltolayer

I attach 2 screenshots:

What could be the problem?


texts: Null output from Transform function of Spatialite

texts:
I am using Spatialite, an extension of SQlite and I have some problems running the transform function.

According to the Spatialite functions page (http://www.gaia-gis.it/gaia-sins/spatialite-sql-4.2.0.html), the transform function should take in a geometry object and an integer of the SRID that I want to transform an object to, and then output a geometry.

I did the following query, but I'm always given a NULL output.

SELECT Transform ((SELECT Geometry FROM healthierdining), 3414)

I have checked that Geometry FROM healthierdining is indeed of type geometry, and yes, I have PROJ.4 version Rel. 4.6.1 installed. What could possibly be wrong with this code?


texts:
I've been trying to automate the merge vector layers via PyQGIS but I can't seem to get it right. This is my code:

import os
import processing

loc = r"C:\Users\loc"
os.chdir(loc)

shapelist = os.listdir(os.getcwd())

for file in shapelist:
    if "R4A.shp" in file:
        R4A = print(file)

Parameter1 = "QgsProcessingParameterMultipleLayers"
Parameter2 = "QgsProcessingParameterCrs"
Parameter3 = "QgsProcessingParameterFeatureSink"

#parameters = "Parameter1: file1.shp; file2.shp; file3.shp; file4.shp, Parameter2: 'EPSG:4326', Parameter3: r"C:\Users\destination"

processing.run("qgis:mergevectorlayers", parameters)

I always get this error:
SyntaxError: invalid syntax


texts: How to Update Control When Layers are Changed

texts:
I have a control that allows the user to switch between satellite/street tiles. I want the text to update to show the next tile text when appropriate. Currently I am using the Control.render hook but this runs whenever the view changes. I just want a function to run when the control is added to a map and whenever the tile layers are modified. Is there an event I could listen to for that?

Here is the constructor of my current code:

constructor(options: any) {

    const button = document.createElement("button");
    const element = document.createElement("div");
    element.appendChild(button);

    super({
      ...options,
      element,
      target: options.target,
      render: () => {
        button.innerHTML = this.getNextLayer().getProperties().label;
      }
    });

    button.addEventListener("click", this.switchTileSet.bind(this));
}


texts: Gdal raster to polygon in QGIS3 (Python3)

texts:
What shall I install to be able to run the GDAL algoritms in python3 in qgis3?

In Windows if I install gdal-python version 2.4.0-3 I get the gdal_polygonize.bat file, but it required python2.7 to be installed, is this a bug or I'm I missing something?


texts: QGIS 2.99 projection for unprojected x y- coordinates

texts: I vaguely remember that I have read that QGIS 2.99 has special support for unprojected x-y-coordinates. (In qgis 2.x I used epsg 4326 for this). Where can I find "no-projection" "unprojectes" or something similar in QGIS 2.99?


texts: Cannot login GeoServer after choosing to install as service?

texts:
I work on windows7. While installing GeoServer, I selected the 'install as a service' option. As I know when I installed GeoServer as a service, geoserver is already running. But I nagative to Web Administration Interface by url:

  localhost:8089/geoserver/web/

I can't log in by username and password. It have an announcement "Invalid username/password combination". Do anyone know this problem and fix it. Help me, please!


texts:
I search different online PPP (PRECISE POINT POSITIONING) services, but I am new in geodesy, so I can find only a few services like:

CSRS-PPP (Canada)

GPS Analysis and Positioning Software (GAPS)

magicGNSS

Are there any other PPP services?


texts: arcpy.SummarizeWithin_analysis error: error in executing tool

texts:
I am trying to identify and sum the lengths of line segments within a polygon boundary using the arcpy tool:

arcpy.SummarizeWithin_analysis('my_gdb.DBO.gridzones', 'my_gdb.DBO.lines', 'my_gdb.lines_within', 'KEEP_ALL', [['Shape.STLength()', 'SUM']], 'ADD_SHAPE_SUM', 'METERS', 'state', 'ADD_MIN_MAJ', 'ADD_PERCENT', 'my_gdb.DBO.inspection_percent')

but when I run this in a stand-alone Python script, I get this error:

arcpy.SummarizeWithin_analysis(in_polygons, in_sum_features, out_feature_class, keep_all_polygons, sum_fields, sum_shape, shape_unit, group_field, addMinMaj, add_group_percent, out_group_table)
  File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\analysis.py", line 2265, in SummarizeWithin
    raise e
  File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\analysis.py", line 2262, in SummarizeWithin
    retval = convertArcObjectToPythonObject(gp.SummarizeWithin_analysis(*gp_fixargs((in_polygons, in_sum_features, out_feature_class, keep_all_polygons, sum_fields, sum_shape, shape_unit, group_field, add_min_maj, add_group_percent, out_group_table), True)))
  File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\geoprocessing\_base.py", line 496, in <lambda>
    return lambda *args: val(*gp_fixargs(args, True))
arcgisscripting.ExecuteError: ERROR 999999: Something unexpected caused the tool to fail. Please refer to the error help for potential solutions, or contact Esri Technical Support http://esriurl.com/support.
Failed to execute (SummarizeWithin).

and when I run it from within ArcGIS Pro (with parenthesis around the feature class names), I get this error:

Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\analysis.py", line 2265, in SummarizeWithin
    raise e
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\analysis.py", line 2262, in SummarizeWithin
    retval = convertArcObjectToPythonObject(gp.SummarizeWithin_analysis(*gp_fixargs((in_polygons, in_sum_features, out_feature_class, keep_all_polygons, sum_fields, sum_shape, shape_unit, group_field, add_min_maj, add_group_percent, out_group_table), True)))
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\geoprocessing\_base.py", line 496, in <lambda>
    return lambda *args: val(*gp_fixargs(args, True))
RuntimeError: Object: Error in executing tool

I am following the documentation verbatim, but cannot get this tool to work, and am not getting a concise error message to troubleshoot this issue.


texts: Getting the center and boundaries of cell of raster object in R

texts:
Raster visualization

Raster information
class      : RasterLayer
dimensions : 81808, 37369, 3057083152  (nrow, ncol, ncell)
resolution : 0.0001843661, 0.0001843661  (x, y)
extent     : 119.7351, 126.6247, 4.554176, 19.6368  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs
source     : /path/to/raster.tif
names      : raster
values     : 0, 255  (min, max)

Information derived from R/RStudio
Objective

Above is the zoomed-in visualization of the same raster. Given some LatLng values with 7-to-8-decimal places (white marker), I think there should be a way to find the center and boundaries locations (black markers) of a certain cell/pixel.
I tried the functions cellFromXY() and xyFromCell() but got no luck because they were just returning the same values like so:
raster.data <- raster( "/path/to/raster.tif" )
given.point <- c( 123.17987768 , 12.09548919 )
cell.target <- cellFromXY( raster.data , given.point )
benchmark.point <- xyFromCell( raster.data , cell.target )

Both the given.point and benchmark.point returned the same values and can not be used this way to get the center values. Which puts me to the question...
How to get the center and boundaries of a cell of a raster object using R?


texts: how to get the projection from a vectorlayer in qgis?

texts:
I am trying to set the projection on a raster to match that of a vector point layer. Thus I need to find out what is the projection of a given layer, to use it in the GDAL.Dataset.SetProjection() so that I can create the GEOTIFF with the appropriate projection.

How do I do That? (in Python)


texts: I want to make a raster showing the average values of 12 rasters. I use Raster Calculator and simply add all of them and divide them by 12 but the result is not correct! The max is far less than what it should be! I already used ArcGIS to do the same job and in ArcGIS I use Mean function and it works as expected and correctly. I'm really confused what is wrong! What is the simple function for mean in QGIS?


texts: Automatic bus route map without overlaps

texts: I'm working with drawing maps of all bus routes in a defined area (municipality) in QGIS. I have a problem about showing two routes with the geometry (the same road) parallel. I know it's possible to do it manually (line offset). Is it possible to do automatically?


texts: Osm2pgql failed to output expired tiles for "delete" operation when running in cron

texts:
I am using osm2pgsql to:

Append .osc files into gis rendering database (based on mod_tile / Mapnik stack) to "delete" some features

Generate expired tiles list via -e -o options

Fed the expired tiles list to render_expired to re-render the tiles

I run all the steps in an script via crontab as the user myself.
The weird thing is:

a. When I run the script manually via login shell (bash), the expired list comes out as expected.
b. However, after run in the crontab job, for osm changes "delete", the .osc file is generated property, but the expired list is always "empty"!

Anyone met with this kind of problem before?
The keys here lie in:

osm2pgsql run by crontab in a bash script
.osc file is not empty
Only "delete" changes works wrong
pgsql database log shows "delete" statement works properly (so as for login & other operations)
I hacked osm2pgsql code, by the logging, m_dirty_tiles is not empty when pasring .osc xml


texts: Reclassify a raster in R based on two rasters with the output containing over 1000 bands

texts:
I am trying to convert the outputs of a model into a raster based on both a temperature raster and a seasonality raster. I have over 2100 bands so ifelse statements aren't going to cut it.

An example would be that I had my model outputs in a dataframe - with the range of temperature values in column 1, the seasonality in column 2 and the outputs in column 3:

a<-rep(1:100,10)
b<-sort(rep(1:10,100))
c<-rnorm(1000,5,3)
df<-data.frame(a,b,c)

And I have a raster of temperature

t = raster(matrix(1:100,10,10))

And a raster of seasonality

s = raster(matrix(rep(1:10,10),10,10))

I need to write a new raster which is based on the values of s and t, but consists of the model outputs of column 3 of the data frame.


texts: Read shapefiles in Python and render on Google Maps?

texts: Is possible to read shapefiles in Python and render them on Google Maps? After research I'm really confused. Should I convert shapefiles to KML?


texts: How to compute the conditional function and output binary image in GEE?

texts:
I am now using Google Earth Engine to process the Landsat 8 images. Now I have calculated NDVI, EVI, mNDWI and generated some maps. I am now want to generate a function based on those indices values, output binary numbers (0 or 1) and layout binary image. But I get stuck, and I hope anybody could help me.
I want to justify the water body, there's a formula:

water = 1 EVI < 0.1 and (mNDWI > EVI or mNDWI > NDVI)

water = 0 For other cases

The GEE guideline says the traditional if statement is not recommended and yes I tried it doesn't work, I used the way they suggested:

var water = ee.Algorithms.If(evi < 0.1 && (mndvi>evi || mndwi > ndvi), 1, 0);
var waterViz = {min: 0, max: 1, palette: ['00FFFF', '0000FF']};
Map.addLayer(water, waterViz, 'water', false);

And it says:"Cannot add an object of type ComputedObject to the map."

I have seen a similar question here, but that user used ArcGIS. I am new to JavaScript.


texts: Combining multiple feature classes into single feature class using ArcGIS for Desktop?

texts:
I am using ArcGIS 10.2.2 for Desktop?

I would like to combine multiple feature classes (same geometry) into a single feature class. The catch is that all of the feature classes have identical fields within them except one, which is unique to that feature class.

Is there a way I can combine these features classes and preserve all of the identical fields with the data correlated to them, while combining the unique field from all of the classes into a new field that just uses the value of the unique field?

I have tried both the append and merge processes. Append unfortunately does not bring any attribution if its not the same field name and merge brought everything, creating a new field for every unique field.


texts: Creating Longitude and Latitude Lines in Google Maps

texts:
I'm working on an independent high school project where i'm trying to build a simple differential equation model for diffusion. I'm trying to overlay the model on top of a map to visualize spread, so for example, oil in water.

Right now I have a rectangular box that have the coordinates of ( 31.0195,-102.8072) NW, (19.2019, -102.8072) SW, (31.0195, -75.3963) NE, and (19.2019, -75.3963) SE. I am trying to build a 2D array of lat/long coordinates spaced evenly. For example, an array of lat/long coordinates that are 500m apart.

Is there a formula/algorithm where I can convert these coordinates to a distance, break the distance up into 500m chunks, and then get lat/long coordinates along the way (i.e. so I have lat/long coordinates every 500m)?

Thanks


texts:
I want to use GeoServer plus OpenLayers in a Liferay portlet. My main question is how to do this? I mean how to include GeoServer in Liferay, how to open GeoServer's admin panel in Liferay.

Next question its how to use default GeoServer's CSS. When I create simple HTML page I use:

<link href="http://localhost:8080/geoserver/openlayers/theme/default/style.css" type="text/css"/>

But in Liferay all CSS must be included in file liferay-portlet.xml can someone say how to add this CSS in this file?

Next question: how to use WMS? In HTML I use this:

layer = new OpenLayers.Layer.WMS( "OpenLayers WMS",
    "http://localhost:8080/geoserver/wms", {layers: 'cite:filedata', transparent: true, format: "image/png"});

I must change this link to WMS to use this in the portlet?

And final question: I try add OpenLayers.js in layfray-portlet.xml but after deploying my portlet I don't see any map controls on the page just an empty DIV.

This is part of liferay-portlet.xml

<portlet>
 <portlet-name>MyGreeting</portlet-name>
 <icon>/icon.png</icon>
 <instanceable>true</instanceable>
 <header-portlet-javascript>/js/jquery.min.js</header-portlet-javascript>
 <header-portlet-javascript>/js/OpenLayers.js</header-portlet-javascript>
 <footer-portlet-css>/css/main.css</footer-portlet-css>
 <footer-portlet-javascript>/js/main.js</footer-portlet-javascript>
 <css-class-wrapper>MyGreeting-portlet</css-class-wrapper>
</portlet>

TThis is what I get:

I forget to mention one thing. I see an error in Firebug with this text:

localhost:8080/poller/receive error 404 Not Found

Anyone have ideas?


texts: Creating new vector layers / save output in QGIS scripting?

texts:
I need to create a script in QGIS which stores outputs. In this case, official documentation didn't help to me.

When I use geoprocessing (processing.runalg()) for vector layers. In some cases, output is selection of the input vector (eg. select by location) and in other cases, output is a vector layers (eg. split by lines with lines).

How do I save both types of outputs as:

new memory layer;
new shapefile saved on disk?


texts: Import Illustrator/Inkscape svg into qgis

texts:
I tried to import lines from svg into qgis 3.x by using dxf format. The problem is that the scale of dxf is of course not the right one (a few centimeters whereas my map is kilometers long).

I tried to look at affine transformation but it is very annoying to use and doesn't really quite work as expected.

Do you have any advice on how to import the drawings into qgis and how to scale them properly?

Or only on the last step, on the scaling and translating issue?


texts: How to get all the column number and row number of DTED file

texts:
I am not very familiar with GDAL. I use this coding to get all the 3601*3601 elevation value into an vector array.The issues that it only show elevation point and not the row and column. How can I get the number of column and row that align with that elevation point?
int _tmain(int argc, _TCHAR* argv[])
{
    const char* pszFilename = "C:\\ElevationSample.DT2";
    double adfGeoTransform[6];
    float buffEleveationValue[1];

     int rows = 3601;
     int columns = 3601;

    std::vector<float> buffElevation(rows * columns);

    int row;
    int col;

    GDALRasterBand *poBand;

    GDALDataset  *poDataset;
    GDALAllRegister();
    poDataset = (GDALDataset *) GDALOpen( pszFilename, GA_ReadOnly );

    printf( "Driver: %s/%s\n",
        poDataset->GetDriver()->GetDescription(),
        poDataset->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME ) );
    printf( "Size is %dx%dx%d\n",
        poDataset->GetRasterXSize(),poDataset->GetRasterYSize(),
        poDataset->GetRasterCount() );
    if( poDataset->GetProjectionRef()  != NULL )
        printf( "Projection is `%s'\n", poDataset->GetProjectionRef() );
    if( poDataset->GetGeoTransform( adfGeoTransform ) == CE_None )
    {
        printf( "Origin = (%.6f,%.6f)\n",
                adfGeoTransform[0], adfGeoTransform[3] );
        printf( "Pixel Size = (%.6f,%.6f)\n",
                adfGeoTransform[1], adfGeoTransform[5] );
    }

    poDataset->GetGeoTransform(adfGeoTransform);
    poBand = poDataset->GetRasterBand(1);

    poBand->RasterIO(GF_Read, 0, 0, columns,rows, buffElevation.data(),columns,rows, GDT_Float32, 0, 0); // read the cell value

            getchar();

            return 0;
    /*  }
    }*/

}


texts: A faster polygon intersection tool than rgeos gIntersection in R?

texts:
I'm looking for a R function able to retrieve the polygon intersecting between two SpatialPolygonsDataFrame.
I need to use R and the rgeos::gIntersection() is too slow in comparison with ArcGis...


texts:
I'm tempted to propose porting some windows forms based arcmap tools to WPF.

Before doing that though, I'd like to hear from others who have done something similar, and what sorts of issues were encountered.

Compared to Windows Forms, how stable is WPF in Arcmap?


texts: GeoServer GeoTIFF not rendering properly in QGIS and GeoServer Previews

texts:
I have a WMS raster ortho image that I have published as a GeoTIFF on two different GeoServers, a development server and a test server. On the dev server the image can be previewed in GeoServer just fine. In the test server the preview is invisible except for the black 'noData' around the perimeter. Both instances of the raster were published the same way in two different GeoServers.

When I open both the WMS layers in QGIS, the image is displayed incorrectly and only as a single band instead of the four bands (RGB + alpha).

What's happening?

This is how it looks displayed in QGIS but also looks that way in GeoServer.

QGIS only recongizes a single band

EDIT: It appears that I can see the layers when I zoom in to certain extents using both GeoServer and QGIS but not at the default extent and not at other extents when the WMS is brought in to QGIS/previewed in GeoServer.


texts: Getting 3D shapefile into 3D view using QGIS

texts:
Despite trying and following several guides/tutorials (and CRS types) I am unable to get a 3D .shp into QGIS 3d view - all i can get is a flat plane with no z values and .shp projected in 2D.
The same .shp open as 3d polylines and polygons in various other applications.
To be clear; i am not looking to extrude or use dem/tin/raster elevations base, but simply render the 3D .shp as-they-are...
Is it even possible to bring 3D .shp (that render unproblematically in ESRI and simple .shp viewers) into QGIS 3d as stand-alone data?
Do I need to import some kind of base elevation data first?
A problem for me is that all the guides i can find are primarily about terrain and elevation maps, and then extruding ... instead i have fully 3d ,shp that just need to be rendered, and i can not find clear information about that.


texts: Performing Raster Reclassification in QGIS?

texts:
Is there a way to perform a reclassify of a raster image in QGIS?

In ArcGIS you can do that in Spatial analyst Tools→Reclass→Reclassify

Any ideas?


texts: Why is this SRTM TIFF not visible in GRASS GIS 6.4.3?

texts:
I want to use interpolated SRTM data for Germany to generate an STL terrain map for use in a CFD application.

My SRTM data comes from the OpenDEM project. SRTM data normally has voids, but this data has been assembled from multiple tiles and the voids filled with other DEM data and through interpolation.

The DEM is in the form of a TIFF, which I am presuming is a GeoTIFF since I've yet to see any other GIS TIFF formats.

When I import the TIFF into GRASS, everything seems to work at first:

(Sat Aug 31 23:22:22 2013)
r.in.gdal input=/home/sfbosch/Documents/Quelle/OpenFOAM/GIS/srtm_germany_dtm/srtm_germany_dtm.tif output=srtm_germany_dtm
Projection of input dataset and current location appear to match
r.in.gdal complete. Raster map <srtm_germany_dtm> created.
(Sat Aug 31 23:22:28 2013) Command finished (5 sec)

The moment I try to make the layer visible in the Layer Manager, though, I get an error:

Command 'd.rast -o map=srtm_germany_dtm@PERMANENT' failed
Details: region for current mapset line 4: <north:
371149.387386>
run "g.region"

The image does not appear in the GIS Map Display window. My first thought was, I've imported the data incorrectly, or there's something wrong with it. But when I use QGIS with the GRASS plugin to access the very same layer and set Contrast Enhancement to "Stretch to MinMax", there it is (though I would like to know why the top half is so dark):

Why does it work in one, and not the other? Unfortunately, it's not good enough that it works in QGIS -- I need to get this working in GRASS because the Python script I am using to generate the STL file depends on it.

One clue might be that, when I move the cursor over the (blank) Map Display window in GRASS after rendering the layer and getting the error message, I can't make any sense of the coordinates. They seem to be in some mish-mash of UTM and LatLon, and they are not consistent. The numbers jump all over the place.

If I get the region, here's what comes out:

(Sat Aug 31 23:49:02 2013)
g.region -p
projection: 3 (Latitude-Longitude)
zone:       0
datum:      wgs84
ellipsoid:  wgs84
north:      56:00:01.5N
south:      0:59:58.5S
west:       0
east:       17E
nsres:      1
ewres:      1
rows:       57
cols:       17
cells:      969
(Sat Aug 31 23:49:02 2013) Command finished (0 sec)

I'm using GRASS GIS 6.4.3, QGIS 1.8.0, GDAL 1.9.2 on Gentoo Linux 3.7.9-gentoo x86_64.


texts: Accuracy/ Positional error in converting GIS data from one format to another

texts:
Is there a chance of error creeping in to the positional accuracy of features/objects when converting from one GIS format to another.

Say from SHP >> DWG >> SHP

If yes, by what amount?


texts: Applying visualisation parameters to all images in image collection when exporting to drive

texts:
I want to apply an image visualisation to all the images from a filtered image collection, and export the images to Google Drive using the GEE Javascript API.

I can do the export fine, but I just end up with binary black/white images that aren't useful, so I need to add a {min:-20, max:1} visualisation to each image in the collection either before or during the export job.

Here is the code I'm using to try and achieve that, adapted from the answers to Export each image from a collection in Google Earth Engine.

// Import batch tools from fitoprincipe
var batch = require('users/fitoprincipe/geetools:batch');

// Define a rectangular image region for filtering the collection and exporting images
var geometry =
    /* color: #d63000 */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[-73.85413391074309, -47.41730172316857],
          [-73.85413391074309, -47.47117019964287],
          [-73.72641784628996, -47.47117019964287],
          [-73.72641784628996, -47.41730172316857]]], null, false);

// Load the Sentinel-1 Image Collection and filter by properties
var IW_H = ee.ImageCollection('COPERNICUS/S1_GRD')
        .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
        .filter(ee.Filter.eq('instrumentMode', 'IW'))
        .filterBounds(geometry)
        .filterDate('2017-07-01', '2017-07-20')
        .select('VV');

// Filter to get images from different look angles
var DescCollection = IW_H.filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'));

// Define visualisation parameters for the exported Sentinel-1 images
var visParams = {min: -20, max: 1};

// Export all images in collection to test folder in my drive
batch.Download.ImageCollection.toDrive(DescCollection.visualize(visParams), 'test',
                {scale: 10,
                 region: geometry,
                 fileFormat: 'GeoTIFF',
                 formatOptions: {
                  cloudOptimized: true
                  }
                });


texts: Convert a tab (Mapinfo) file to ESRI Shapefile with ogr2ogr

texts:
I want to convert a tab file created with MapInfo to an ESRI Shapefile (.shp) using ogr2ogr in command line.

I did the following

ogr2ogr -f 'ESRI Shapefile' "file.shp" "file.tab"

but it gives an error:

Unable to open datasource `file.tab' with the following drivers.

I use gdalv2.4

Do I have to specify anything else?


texts: How do i install ESRI File GDB API in Ubuntu 16.04 so QGIS 2.8 can see it?

texts:
Can someone tell me how to get QGIS 2.8 in Ubuntu to see the ESRI File GDB API?

I am trying to install the ESRI File GDB API (https://github.com/Esri/file-geodatabase-api/tree/master/FileGDB_API_1.5) but i seem to be unable to do so.  I know that QGIS has the open fileGDB, but that won't let me edit.  I am not finding a lot of instructions out there besides the readme provided by the API

I followed the instructions in the API to make the samples, but these result in errors such as:

../../lib/libFileGDBAPI.so: undefined reference to `GetFileAttributes(wchar_t const*)'
../../lib/libFileGDBAPI.so: undefined reference to `IID_IGeometry2'
../../lib/libFileGDBAPI.so: undefined reference to `FindNextFile(void*, _WIN32_FIND_DATA*)'

I am running Ubuntu 16.04.4 on a virtual box vm. The GDAL version for this machine is 1.11.3, released 2015/09/16. My system buildout is as follows:

    QGIS version    2.8.6-Wien  QGIS code revision  exported
    Compiled against Qt 4.8.7   Running against Qt  4.8.7
    Compiled against GDAL/OGR   1.11.3  Running against GDAL/OGR    1.11.3
    Compiled against GEOS   3.5.0-CAPI-1.9.0    Running against GEOS    3.5.0-CAPI-1.9.0 r4084
    PostgreSQL Client Version   9.5.13  SpatiaLite Version  4.3.0a
    QWT Version 6.1.2   PROJ.4 Version  492
    QScintilla2 Version 2.9.1


texts: Finding layer file used to add layers to MXD in ArcPy?

texts:
I work with layer that located in G:\PROJECTS\bneyAyish\gis\lyr\New File Geodatabase.gdb\polygon, and I created a layer file and located it in G:\PROJECTS\bneyAyish\gis\lyr under the name "indisLabel":

I added the layer file "indisLabel" to some of the MXD files and after that i used this code in order to see if this code will find all the maps that have this specific layer file:

import arcpy
from arcpy import mapping as m
from os import path, walk

root_directory = r"G:\PROJECTS\bneyAyish\gis"
path_to_find = r"G:\PROJECTS\bneyAyish\gis\lyr\indisLabel.lyr"
def FindMaps(root_directory, path_to_find):
    maps = []
    for root, dirnames, filenames in walk(root_directory):
        for fname in [f for f in filenames if f.endswith(".mxd")]:
            mxdPath = path.join(root, fname)
            if not path.isfile(mxdPath):
                print "didn't found"
                continue

            mxd = m.MapDocument(mxdPath)
            df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
            for df in m.ListDataFrames(mxd):
                for lyr in m.ListLayers(mxd, data_frame=df):
                    if lyr.supports("DATASOURCE"):
                        if lyr.dataSource == path_to_find:
                            print(mxdPath)
                            maps.append(mxdPath)
    return maps
FindMaps(root_directory, path_to_find)

And it doing nothing:

>>> ================================ RESTART ================================
>>>
>>>


texts: Use restriction on locating facilities and route finding

texts:
On a project I came to the following problem in ArcGIS 10.6.1:
I have a road network which consists of paths, local roads, national roads and highways. I have also a bunch of possible facility locations. Some of the facilities are located near a highway. Since the facilities cannot have a direct access to highways, I want to prevent the facilities to be located on highways.
I found this question already here: Restricting portions of road network from analysis in ArcGIS Network Analyst?
I want to apply @Alex Tereshenkov idea of not locating the facility near those restricted routes, but to still let them run through them. I added in my road network layer a coloumn called RESTRICTED (Type: Integer) with the following idea behind:
0 = no restriction
1 = prohibited to locate facilities on this road

I fail however to define in the network dataset that the facilities are prohibited to be located on highways.
What am I missing here?
The very same problem to tell the network dataset that some roads are unilateral while others are bilateral.


texts: Is "puntal" an accepted and defined geometric term?

texts:
There are definitions for the geometric qualities of polygons and lines:

  polygonal: Having many angles; hence characteristic of a polygon.

Ref: http://en.wiktionary.org/wiki/polygonal

  lineal: of, relating to, or consisting of lines; linear.

Ref: http://oxforddictionaries.com/us/definition/american_english/lineal

However, I have been unable to find a scholarly reference for "puntal", that is, having the geometric qualities of a point.

The only solid references to "puntal" (in English) I have found are in the discussions and source code of projects like GEOS and JTS Topology Suite, where there might be classes or interfaces corresponding to Polygonal, Lineal, and Puntal.

Is there a canonical definition of the geometric sense of this term?


texts: After seeing the $600 price tag for Adobe Illustrator, I am especially interested in learning open source solutions for finishing maps created in a GIS.  What open source, or commonly available, software do you use to add the finishing touches (e.g. arrows, textures, colors, graphics, text, etc) to maps created in Arc, QGIS, or your choice of GIS software? Examples welcome.


texts: Auto-updating map service via REST in ArcGIS for Server?

texts:
I have one map service which is updating per day so I need to manually overwrite the service in ArcGIS 10.1 for Server.

In previous versions 10.0 its very easy..just update the .mxd then REST will be auto updated.

Is there any way to auto-update the map service & REST  in ArcGIS 10.1 for Server

I tried the ADMIN API but its not the proper way ... also there are no settings for service definition file in 10.1.

Any help/ideas/workflow?

The problem is resolved by adding arcgis account.


texts: How to visualize pgRouting layer in leaflet?

texts:
I have completed all the queries of the of pgRouting tutorial by Underdark where visualization was done with the help of RT SQL Layer Plugin. As I am using Windows OS I have no access of that tool. That is why I was unable to visualize the resultant layer of Shortest_path. If there any way to visualize this resultant layer through leaflet API it will be very helpful.

Thanks in advance.


texts: Identifying unknown coordinates/coordinate system

texts:
I have a set of data for the country of Chile which has geospatial information in the form of 'XX0000X' where X is a letter and 0 is a number. I have not been able to figure out what these coordinates are in my googling.

Could anyone please point me in the right direction?

Example:

GEOREF

GT5346D
GL9366D
GK7748C
HD3342D
IB7926D
IB7926D
ZQ6108B
IB7932D
GL9366D
GL9366D
YG8782A
GB1982B
YG8374B
DE1302A
GL9368C
ZM5324B


texts: Installing ArcGIS Web Adaptor without Windows Authentication?

texts: Is is possible to install ArcGIS Web Adaptor without Windows Authentication? I am installing ArcGIS Web Adaptor for Server and Portal. I need to enable ISS components and I realized that I do not have Windows Authentication. I am using Windows 10 Home edition. How should I go about installing Web Adaptor as the installer will not let me proceed without enabling all the ISS components requirements.


texts: Geodatabase Featureclass not editable in Python

texts:
I have a geodatabase featureclass that I need to add a field to. I can open this file in ArcMap and make all the edits I need, adding fields, calculating those fields, etc., but when I try to replicate the results in Python using the Add Field tool, I get this error message:

ERROR 000499: table is not editable
Failed to execute (AddField).

The file is stored in a file geodatabase, and I've made sure no other Arc programs are accessing it, so locks shouldn't be an issue as far as I can tell.
Why would I be able to edit the file in one program but not the other? How can I fix this?
It seems like it also works to bring up the python window within ArcMap/Catalog. The problem only exists when I use the Python Shell outside of Arc.


texts:
I have trail data in shapefiles.  Each trail record includes a "difficulty" field, and I would like to use this field to display a standard ski slope difficulty icon next to each trail.  I have these icons as SVG files.

This seems like a labeling matter to me.  I have, for example, use trail length fields to display the length next to the trail.  And of course I have the trail name displayed as well.  I would like to display difficulty, not as the numerical values of "1", "2", "3", or "4", but as a different SVG icon based on whether the trail's difficulty field is set to 1, 2, 3, or 4.

But I can't find a way to do this in QGIS 2.18, in which "labeling" seems focused on text data and it's formatting.


texts: How to remove layers that are not in visible df extent

texts:
I am working on ArcMap 10.3 and I have a map with 50 layers in it.

I am trying to check if layers are in a data frame.

Unlike How to check if layers in a dataframe with Arcpy, I want to define the area to the current data frame.

My code is:

import arcpy,os,sys,string
import arcpy.mapping
from arcpy import env

env.workspace = r"C:\Project"
for mxdname in arcpy.ListFiles('*.mxd'):
     print mxdname
     mxd = arcpy.mapping.MapDocument(r"C:\Project\\" + mxdname)
     df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
     data_frame = df
     for lyr in arcpy.mapping.ListLayers(mxd,"",df):
          if lyr in data_frame == False:
               arcpy.mapping.RemoveLayer(df, lyr)
               print lyr
     mxd.save()
del  mxd

and I get:

 >>>
Project.mxd
>>>

but no layers have been removed.


texts: Drawing order of lines and polygons (line between polygon)

texts:
I have 2 layers. One with forests and lakes (polygon) and one with rivers (line). Now I want the features to draw in that order:

The river have to be UNDER the LAKES and the forest have to be UNDER the river and OVER the lakes. Tricky?

Because I have to export as vector-pdf, I can not use the feature blending modes. (or is there a hidden way to export blending modes to vector?
I need the borderlines of the lake and the river.
I can not physically change the layers.

...or is there a other way to draw the river flow into the lake?


texts: around filter: at least x meters away from anything

texts: With the around filter I can find things close to anything for example: way(around.golf_course:250)["natural"="beach"]; How can I form the opposite? I need things who are at least minimum x meters away from anything.


texts: Why doesn't my interpolation result in qgis look good?

texts:
My goal is to use some RTK GPS data from combines or sprayers to generate an interpolated DEM of fields for agricultural purposes. I used QGIS for that. So far I opened the data and did a simple TIN interpolation with the interpolation tool in the raster drop-down. As input, I took the point layer and as the attribute the elevation. Somehow the result always looked like that:

Does anyone know why it doesn't work? I did the same steps as in the qgis tutorial about interpolating point data.

here is a part of the sample point attribute table:


texts: QGIS 2.0 Plug Repository Connection Failed

texts:
I'm trying to connect to the QGIS Official Plugin Repository but the status within the Plugin Manager window is always"unavailable" for the URL - http://plugins.qgis.org/plugins/plugins.xml.

I'm running QGIS 2.0 within Windows 7.

Could the problem be that QGIS is unable to detect the proxy settings used on the corporate network?


texts: Pyproj transform returning incorrect values (EPSG 3875 to EPSG 4326)

texts:
I am trying to convert from EPSG 3875 to EPSG 4326 and have this code:

from pyproj import Proj, transform

def convert_coordinates(input_format, output_format, x, y):

    input_format = "+init=epsg:" +str(input_format)
    output_format = "+init=epsg:" +str(output_format)

    inProj = Proj(input_format)
    outProj = Proj(output_format)

    return transform(inProj,outProj,x,y)

x, y = 4152144, 219039
x2,y2 = convert_coordinates("3875", "4326", x, y)
print x2, y2

This returns 112.264404642 154.392043378 which is out of bounds for EPSG 4326. But ArcGIS returns 37.2993441706517, 1.96727415792103 which is in bounds and looks better.

Is there anything obviously wrong with this code?