texts:
I'm using GeoServer to provide maps generated from raster data (satellite imagery from Sentinel) to a web application via WMS layers.
Each raster consists of a file of about 1.8-2 gb. To save space on disk I tried to apply a gdal_translate using some compression algorithms.

I tried a simple gdal_translate -of GTiff -co "COMPRESS=JPEG" infile.tif outfile.tif.
Input GTiff file is obtained by extracting RGB band from the L2A source. Using JPEG compression each file reduces the dimensions to less than 500 mb.
Testing the web application (built using OpenLayers) I noticed a marked deterioration of performance in terms of speed and responsiveness.
Is the problem due to the computation required to GeoServer to process a compressed file? Is there a better solution among the compression methods provided by GDAL or the best possible performance requires a non-compressed raster?


texts: ArcPy publishing to image server failing at UploadServiceDefinition_server

texts:
We are trying to automatically publish datasets into our image server.

We have an environment with a Portal enterprise server, an image server, a raster server and a machine running ArcGIS Pro. I have tried this running on both the image server and the ArcGIS Pro machine.

We are following this documentation example, slightly modified so it only takes a single file rather than all the files in a directory, and gets the parameters from other places.

It runs through with out any problems to the arcpy.UploadServiceDefinition_server call. At which point it fails with the dreaded ERROR 999999

If I try and upload the service definition 'manually' from ArcGIS Pro it works OK. I am using the same .ags file in both cases.

My best guess is that my script is not logged into the portal correctly. I've attempted to use arcpy.SignInToPortal but it seems to have made no difference to the error I am getting.

Does this process only work if run inside of ArcGIS Pro and not the Esri Python environment?

Is there a step that I'm missing when running outside of ArcGIS Pro?


texts: I have multiple polygon and line shapefiles on which I would like to apply a translation so I was wondering if it is possible to batch process them (maybe with graphic modeler ?) or do I have to change them one by one using affine transformation tool in QGIS 2.2?


texts: How to make the LoCoH toolbox work in ArcMap 10.2

texts:
I'm trying to estimate home range with the Local Convex  method. I cannot make the LoCoH toolbox work in ArcMap 10.2. I get following error:

Executing: alocoh "C:\Users\Linn\Desktop\Fra Bjørn\Export_Output.shp" 1 3 1 100 # #
Start Time: Mon Nov 10 13:56:12 2014
Running script alocoh...

Traceback (most recent call last):
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\ArcToolbox\Toolboxes\locoh.py", line 14, in
    R.Init('R')
  File "", line 2, in Init
com_error: (-2147352567, 'Undtagelse opstod.', (0, None, None, None, 0, -2147221485), None)

Failed to execute (alocoh).
Failed at Mon Nov 10 13:56:12 2014 (Elapsed Time: 0,12 seconds)

Could it just be that my version of ArcGIS is to new?


texts:
I'm using QGIS 1.8.0 and as CRS EPSG:6372 a projection for Mexico, defined by the Instituto Nacional de Estadística y Geografía (INEGI). I have the following questions:

If I measure the minimum distance between points projected with that CRS, what are the units? meters?!
I'm also doing some Vector Grids with maps projected with that CRS, when QGIS asks for the parameters for the grids (as in the attached image) are the units meters? So from the image I will be dividing the map into 10000 x 10000 m. squares, right?


texts: Exporting label rendering options based on zoom level correctly in QGIS

texts:
I am making a webmap using QGIS3.
However, I set the text-label to be invisible from a certain zoom-level onward, and visible from another zoom-level. It displays correctly locally in QGIS edit screen but it does not work after I export it to web. (using qgis2web plugin)
Goal: to make text disappear at 1:1800 scale(zoom) level and above. (text appear below 1:1800 when zooming out)
This is what I mean:

Result after exporting is that "Lot number and land area" will not even show up at all.
I try both OpenLayers and Leaflet.
(this is the example with no label display setting based on zoom level)
--> export outcome: https://landing-at-ven.web.app/kelson_output_1/index.html
--> Ideally, I want to hide the appearance of text label totally when I zoom out to a certain level.
Bad impression:

I wonder is there anyway I can modify the output code to modify and realize what I want to achieve or am I getting some settings wrong?


texts: Kriging on GRASS+R - cannot allocate vector of size 146.5 Mb

texts:
I'm trying to do a groundwater elevation map with an Ordinary Kriging on R, but always resulting on the lack of memory

  cannot allocate vector of size 146.5 Mb

allocation (windows xp 32bit). Should I increase memory size and how, or should I reduce the grid size and how (kind of a newbie in R)?

following this example http://casoilresource.lawr.ucdavis.edu/drupal/node/438


texts: raster tile - find the correct tile for a given coordinate

texts:
I have a large raster map which I loaded into postgis with 1000x1000 meter tiles... the corresponding table has about 200 rows, each row representing a tile. I want to get the raster value for a specific coordinate, but I don't know how to do that most efficiently.

I did something like:

select ST_Value(rast, ST_Transform(ST_GeomFromText('POINT(lat lon)',4326),mysrid)) from rasterdata;

but it gives me a lot of:

Attempting to get pixel value with out of range raster coordinates:...

messages because finally only one tile contains the coordinate that I'm searching for.


texts: Creation of buffer area of 20mtrs upon click event on the base layer

texts:
I have a base layer upon which i wish to import a kml with certain lat lons and create a buffer region of 20mtrs around any point upon click event.

i got google api for the same but want openlayers api to execute this on my base layer.


texts: Performing viewshed analysis in ArcGIS Desktop?

texts:
I want to run a viewshed analysis on a site I am interested in. I have the relevant elevation raster (taken from the OS UK, which is in the ASCII format) and a point feature class from the lat/lon coordinates of the particular site (I have also converted both into the same coordinate system). However when I try to run a viewshed analysis in ArcGIS 10.1, it doesn't work and there is an error message saying that 'there is no valid observation points present'.

I was wondering whether anyone could tell me what I am doing wrong and how to correctly create observation points and run the viewshed analysis. I have looked at some books on the subject and have successfully followed their instructions on the data provided by them but I seem to be having trouble when I am using my own data! Any information would be massively appreciated!

I have used the Geographic Coordinate system for both: OSGB_1936. I have taken the DTM from OS UK website (which is called Panorama and it is in ASCII format). First I added both to ArcMap and created a Hillshade. Hillshade did work but there was a yellow exclamation mark and obviously I got the error message I mentioned earlier when I tried to do the viewshed.

Today I changed the coordinate system of the DTM to the British National Grid (that’s under Projected coordinate systems) but left the XY point the same: Geographic Coordinate system OSGB_1936. When I ran the Hillshade with this combination it worked properly i.e. there was no yellow exclamation mark.

Then I ran a viewshed and it seems like this time it worked but the only thing I did different was change the coordinate system for the DTM. I don’t understand why that would make such a difference and I was hoping whether you would be able to let me know if I have done the right thing here (I did try to attach printscreen shots to make the steps I have done clearer but it didn't work - so, sorry about that). I wanted to add 15m and run another viewshed. So I set a OFFSETA (which appeared on the attribute table of my XY point) and another viewshed was generated - indicating that the area of visibility had increased which seems to suggest that it worked. But I was wondering whether you might be able to tell me why changing the coordinate system may have had such an impact and whether I have missed any crucial steps along the way?


texts: Finding all ways that are partly within a bounding box

texts:
I'm trying to write a query that will return all OSM ways which are contained within a bounding box of min/max gps coordinates.

I get these cooridinates from the greater london bounding box on http://wiki.openstreetmap.org/wiki/Bounding_Box

From googling I have this query, but it doesn't return any results:

SELECT
    planet_osm_roads.*
FROM planet_osm_roads
WHERE ST_Contains(
    ST_SetSRID(
        ST_MakeBox2D(
            ST_Point(0.3340155, 51.2867602), ST_Point(-0.5103751, 51.6918741)
        ),
        4326
    ),
    planet_osm_roads.way
)

Any idea what I'm doing wrong?

Update

So I've swapped the points around, tested that the polygon created is what I expect by getting the geometry as geojson and visualising it on geojson.io - it looks fine, however I still get 0 results and it happens basically instantly (which seems weird). The table is definitely fully populated! I also changed to use ST_Intersects but that hasn't made it return anything either

SELECT
    planet_osm_roads.*
FROM planet_osm_roads
WHERE ST_Intersects(
    planet_osm_roads.way,
    ST_SetSRID(
        ST_MakeBox2D(
    ST_Point(0.3340155, 51.2867602),
            ST_Point(-0.5103751, 51.6918741)
        ),
        4326
    )
)

It gives me this GeoJSON when I SELECT ST_AsGeoJSON

{
  "type": "Polygon",
  "coordinates": [
    [
      [-0.5103751, 51.2867602],
      [-0.5103751, 51.6918741],
      [0.3340155, 51.6918741],
      [0.3340155, 51.2867602],
      [-0.5103751, 51.2867602]
    ]
  ]
}

Which as I said, looks correct on geojson.io :)

Update #2

I've copied a query from another SO question that gets the closest point on a road or something, with a point in London, it returns some results roughly as I expect

SELECT
 osm_id,
 highway,
 name,
 ref,
 ST_X(ST_ClosestPoint(ST_Transform(r.way, 4326), point.geom)),
 ST_Y(ST_ClosestPoint(ST_Transform(r.way, 4326), point.geom)),
 ST_Distance_Sphere(ST_ClosestPoint(ST_Transform(r.way, 4326), point.geom), point.geom)
FROM
  planet_osm_roads r,
  (Select ST_SetSRID(ST_MakePoint(-0.1367229, 51.5377733), 4326) as geom) point
ORDER BY 7 ASC
LIMIT 10;

returns

41383228;"trunk";"Camden Street";"A400";-0.137183850111055;51.5373791549638;54.195918203
253353304;"primary";"Royal College Street";"A5202";-0.13568782069246;51.538616921758;118.002625807
17513372;"primary";"Royal College Street";"A5202";-0.134746843460286;51.5374616146768;140.996741591
166369355;"primary";"Royal College Street";"A5202";-0.136110935219222;51.5391336013284;157.068830304
253352352;"primary";"Royal College Street";"A5202";-0.136578148998492;51.5397000983131;214.484260552
6247747;"primary";"Pratt Street";"A503";-0.140021211819466;51.5376420920654;228.588127843
211950233;"primary";"Bayham Street";"A503";-0.140021211819466;51.5376420920654;228.588127843
253350830;"primary";"St. Pancras Way";"A5202";-0.135177855133607;51.5396274637255;232.22039167
11266300;"primary";"Royal College Street";"A5202";-0.136734905015571;51.5398841985372;234.722999647
253350834;"primary";"St Pancras Way";"A5202";-0.135399101773244;51.5397738652997;240.557396864

The point is actually closest to Mandela Street which runs next to Camden Street, but I'll worry about that after I have my initial query working.

My data came from download england osm data and running this command:

osm2pgsql --create --database england_osm -C 6144 ~/Downloads/england-latest.osm


texts: How to use ITransform2D on dimension features (arcobjects)

texts:
I am trying to move a dimension feature inside a personal geodatabase using ArcObjects VBA. I was able to move other types of features using similar code as shown below. However it is not working for dimension features. Can someone help me on this?

Sub AdjustDimensionFeature()
Dim pMxDoc As IMxDocument
Set pMxDoc = ThisDocument

Dim pmap As IMap
Set pmap = pMxDoc.FocusMap

Dim pFtrLyr As IFeatureLayer
Set pFtrLyr = pmap.Layer(0)

Dim pFC As IFeatureClass
Set pFC = pFtrLyr.FeatureClass

Dim pFtr As IFeature
Set pFtr = pFC.GetFeature(8362)

Dim pTransform2D As ITransform2D
Dim pPolygon As IPolygon
Set pPolygon = pFtr.ShapeCopy

Set pTransform2D = pPolygon

pTransform2D.Move 20, 20

Set pFtr.Shape = pPolygon

pFtr.Store

End Sub


texts: Efficient way to handle NA values from analysis in R with Raster

texts:
I am not 100% sure if GIS is the correct stackexchange for this, but here we go.

A problem I have been bested by:
I have a raster layer of approximately 240 million cells.  It could be a DEM or anything similar.  I also have a series of sampling-site polygons of irregular shape located throughout the extent of the DEM but only occupy about 0.1% of the space.  I am using statistical tests to compare the distribution of values (elevations in this example) between the sample-site locations and the overall DEM (considered background values).

I have used the MASK() function in the Raster package of R to extract the the DEM values where they intersect non-Null values of the sample-site raster.  It returns a raster that is 99% NA and the remainder are the values at the sample sites.

From here, I draw roughly 50K random samples from the DEM and I test them against the distribution of values at sampling-site locations.  However, I have only found the most inefficient ways to remove or ignore the NA values (eg. x <- na.omit(y))

*edited to add:
What is an efficient way to extract the values from the sampling-site raster into a vector that contains no NA values?

I aware of the na.rm flag fro some functions in the raster package, but I have not used it to success.
I have successfully implemented this general process by turning the sampling-site raster into cell centroids within sampling-sites and extracting values, but that has extra steps.
I apologize for being so verbose, but I would appreciate any and all ideas.

edit*  this is the most effective way I have found is site_smpl_freq <- as.data.frame(freq(site_smpl, useNA='no', progress='text')  but it still takes a few minutes to run

edit* to add a code example.

    ###sample background data and put in data frame
    rand_smpl <- sampleRandom(raster("DEM"), collapse =''))),50000)
    rand_smpl_freq <- as.data.frame(count(rand_smpl))

    ### then a bit that modifies the table if values exceed a certain threshold

    ### extract raster values at the location of sampling-sites
    sites_extract <- mask(DEM, sampling_sites_raster)

    ### Here is where I create the freq table using the useNA option
    ### It works, but it takes a long time and seems inefficient
    site_smpl_freq <- as.data.frame(freq(sites_extract, useNA='no', progress='text'))

    ### then I modify the table as above if values exceed a certain threshold
    ### This is followed by some manipulations that add columns for cumulative counts and percentages, then to ggplot()


texts: Weighted Overlay Decimal Values for Evaluation Scale

texts: I would like to do a weighted overlay of different raster datasets. I would like my evaluation scale to be from 1 to 10 with 0.1 as the incremental value. This way for example a criteria could have the value 1.2. So far ArcGIS does not let me enter decimal points. Any ideas how I can solve this?


texts: How to append text string information from CAD annotation to polygon feature class attribute by location?

texts:
I have a CAD file with polygons. They are signed with an annotation.

Is it possible to assign text values from the annotations to polygon feature class attribute.

As I understand only the location connection between the text and the polygon.


texts: Bay Area Synthetic Network node/line coordinates

texts:
I just started working with the Bay Area Synthetic Network and was wondering if this dataset has location (lat/lon) information. I have visualized the files with QGIS, but when it comes to exporting the data, I have not found a solution to obtain exact locations.

I understand this data is synthetic, so I was wondering if it has any real data for coordinates or the node location is with respect to other nodes.


texts: Problem running shell command with PyQGIS

texts:
I am building a QGIS plugin with Python (3.7) that is supposed to run a shell command.
This command launches an FME Script with some parameters
os.system("fme D:\script.fmw --param1 D:\file.txt --param2 value")
When i run this command on my python plugin main script, there is no effect.
But when i run the same command in a shell, it works. When i run it from a regular python script, it works too.
I also tried using subprocess module to get more information on my result. It raises a compilation error with code 1.


texts: How to create multiple graduated layers using different attributes/columns?

texts: I have a shapefile for which I'd like to create multiple graduated layers for each attribute/column. At present, I manually duplicate the layer and then change the column from the drop down list, select a color ramp and update the classes. As I have numerous layers to create, how do I go about automating these steps? Can I do this using a python script?


texts: GeoPackage in QGIS: what are .gpkg-shm and .gpkg-wal files?

texts: In QGIS, if I open some GeoPackage file like mylayers.gpkg, it will create two files named mylayers.gpkg-shm and mylayers.gpkg-wal in the same folder. After I close QGIS, these files disappear. I understand that they are temporary files, but I wonder what their purpose is.


texts:
I am using QGIS 3.16.3.
I am trying to edit an existing thematic shapefile layer composed of contiguous polygons and polygons in polygons. When I draw a new polygon within or crossing existing polygons or part sharing boundary; to get the properties correct, I need to 'clip' to remove existing attributes within the new polygon. Clip doesn't seem to work. I tried the Plugin Clipper but this did not work. Clip is an essential tool for editing polygons in polygons in ArcGis. I cannot seem to edit a polygon shapefile in QGis without it.
Is Clip in later versions of QGIS or is it only on a wish list?


texts: Standalone GeoWebCache drop some parts of features in tiles

texts:
I have a local GeoServer with standalone GeoWebCache to serve vector tiles in PBF file format. The layer displays country borders in the world but the problem is in different tiles some parts of polygons missed and have not shown. You can see the difference between the two following images.

The geoserver logs these line but I cannot find any layer group named web:

2020-05-05 17:14:22,133 DEBUG [filter.GeoServerSecurityContextPersistenceFilter$1] - SecurityContextHolder now cleared, as request processing completed
2020-05-05 17:14:36,404 DEBUG [geoserver.wps] - Removing statuses matching [[[ NOT [ completionTime IS NULL ] ] AND [ completionTime Before 2020-05-05T17:09Z ]] AND [[ NOT [ lastUpdated IS NULL ] ] AND [ lastUpdated Before 2020-05-05T17:09Z ]]]
2020-05-05 17:15:15,606 DEBUG [security.IncludeQueryStringAntPathRequestMatcher] - Checking match of request : 'Path: /web/wicket/bookmarkable/org.geoserver.web.admin.logpage, QueryString: filter=false'; against '/web/'
2020-05-05 17:15:15,606 DEBUG [security.IncludeQueryStringAntPathRequestMatcher] - Matched Path: /web/wicket/bookmarkable/org.geoserver.web.admin.logpage, QueryString: filter=false with /web/
2020-05-05 17:15:15,609 DEBUG [geoserver.ows] - Could not find a layer group named web
2020-05-05 17:15:15,609 DEBUG [geoserver.ows] - Could not find a layer group named web
2020-05-05 17:15:15,609 DEBUG [geoserver.ows] - Could not find a layer group named web
2020-05-05 17:15:15,609 DEBUG [geoserver.ows] - Could not find a layer group named web
2020-05-05 17:15:15,610 DEBUG [org.geoserver] - Thread 142 locking in mode READ
2020-05-05 17:15:15,610 DEBUG [org.geoserver] - Thread 142 got the lock in mode READ


texts: Can NIR minus Red be used instead of NDVI with Drone Imagery

texts:
I have some CIR imagery collected by a vendor using a Cannon Powershot S110 NIR, from an eBee fixed wing platform, about 400 ft AGL. All images were stitched together by the vendor in Pix4D. I want to use it to calculate NDVI for a change detection project in a small riparian corridor. A cursory visual inspection of the NDVI values looks about right, except I'm getting very high values in areas where there are shadows. The red band in these areas is picking up very little visible light, and many more times the amount of NIR.

I've read that NIR minus Red has been used instead of NDVI for drone imagery with this problem, but am not finding much literature on it. It seems like if the red band is compromised, this method is unlikely to produce much more reliable results.

Is this is a viable option?

If not, what other approaches might I try?


texts: Python - gdal.Polygonize produce invalid geometries

texts:
gdal.Polygonize is making invalid geometries and later in processing chain an error is reported stating that there are Ring Self-Intersections. Not one, quite a few so manual solution is not an option in my case. Is there any solution to this problem? I am using Python 2.7 and QGIS. Hot spot (ring self-intersection) is circled in red.

Question is the same like QGIS Raster to Polygons produce invalid geometries, but I am looking for some solution using Python and its libraries.


texts: QGIS and QField show different behaviour regarding a Value Relation with Filter Expression

texts:
I added a Value Relation Widget in the Attributes Form of a shapefile in order to generate a drop-down menu that is supposed to make fieldwork more efficient. The idea is to have an optional filter that limits the options in this drop-down.
More specifically, I have a first drop-down in which I can select plant families. When a plant family is selected, I want to limit the options in the second drop-down (Genus) to the genera that are part of this plant family.
I have very limited experience with SQL, so this might be complete nonsense. But the following filter expression appears to work in QGIS desktop:
"Family" = current_value( 'Family' ) OR
current_value( 'Family' ) IS NULL

When no Family is set, all the options for Genus are available. When I select a Family, the options in Genus are limited to genera that are part of this family.
Now comes the problem:
When I export the project to QField, the filter no longer works. In QField, I can select any Genus at any time (independent of the Family that I have or have not selected previously); moreover it appears that selecting a Genus can even reset the Family value.
Here some screenshots of desktop and mobile version of the project. In both cases I selected Aquifoliaceae as the Family which should limit the Genus options to Ilex in this example. On the mobile device, the options were not limited and selecting Acacia, for example, did reset the Family to NULL. Subsequently, I selected the Family again, with no effect on the Genus.

Is there a way to solve this issue?
Edit: I changed some plant family entries and exported the project again. Now the selection of a Family limits the Genus options, but the Genus doesn't limit the Epithet options. So basically, the problem moved to another level. I'm not sure what causes that issue but I suspect it is a bug in QField?


texts: Displaying field name istead of ID in Attribute Form filter using QGIS

texts:
I have three Postgresql layers:

1- item_line(id,geom,job_id)

2- job(id, name, status_id)

3- status(id, name)

In QGIS I have created two relations: one is between item_line and job, and the other is between job and status.

Then in the Attribute Form of the item_line layer, I configured it to show the job name instead of the job_id, as shown in the following screenshot.

And created a fitler on the job_id widget, based on the status_id, as shown in the following screenshot.

What I want is instead of showing the status_id in the filter combobox, I want to show the status name.

How to Do that?


texts: What's the difference between these two different ways of raster stacking with python QGIS

texts:
I've found two different ways for creating a multiband raster using python and qgis.

Both the methods produce the expected results. The dimension of the multiband raster file created with the first method is double of the dimension (in terms of byte) of the one produced with the first method.

Which method is preferable?

import os
import processing

# first method
processing.runalg("gdalogr:merge", "/path/file1.tif;/path/file2.tif;/path/file3.tif", False, True, 0, 6, "/path/fileout.tif")

# second method
cmd = "gdal_merge.py -separate -o" + " path/outfile.tif " + " /path/file1.tif " + " /path/file2.tif " + " /path/file3.tif"
os.system(cmd)


texts: Get tiles count for shapefile in specified zoom level

texts:
I have polygonal shapefile, which contains several thousands of features. I need count tiles for 18 zoom level which generates for each feature in shapefile.
I find Tile Calculator http://tools.geofabrik.de/calc/#2/28.5307/62.5469&type=Geofabrik_Standard, but there possible get count tiles for one rectangle area (in shapefile this areas many).

How can I calculate the number of tiles for the shapefile at a certain zoom?


texts: Can ArcMap automatically highlight object on map using its label?

texts:
I was wondering if there is a way to highlight a feature (e.g. line) only by selecting its label on the map.

Can ArcMap autmatically highlight this object or do I need to create some "connection" between line and its label?


texts:
Possible Duplicate:
  Can I store an ArcObject inside a BLOB or XML field?

I am saving shape file in database and now i want the symbology selected for the feature layer also get stored in database. I am facing an issue I cannot serialize IStyleGallery, and if its not serialized than how can i store that in database. Anyone got any good solution?


texts:
I am attempting to download a bunch of 4 band basic scene data, the search returns results, however, the asset link is just an empty dictionary.

When I search for the image ID on scenes explorer, the image is available.

    scenes = results.json()['features']
    if len(scenes)!= 0:
        image = sess.get(scenes[0]['_links']['assets'])

there is an asset link:
https://api.planet.com/data/v1/item-types/PSScene4Band/items/20171212_054612_1022/assets/'

However, there is not image at this link nor a link to activate the image


texts: One Polygon Joining to Database with Many matching records

texts:
Using ArcMap 10.1.
Is this possible to have one polygon and join to a database/excel with many matching records?
I would like to have 1 polygon but match to 100 corresponding records in a table.  When I use the identify tool or select this polygon all 100 matching records would display?


texts:
I hate to ask this kind of questions, but i ve been stuck for 2 days now:

I rendered some tiles with mapnik which i want to serve in my offline application which will run from a portable disk, but via browser still. I switched openlayers to leaflet for this particular case because of its nice widgets and seemingly easier setup.

However, i initialize my map and tile layer like this:

var map = L.map('map_canvas');
map.setView([50, 10.0], 6);

var layer = L.tileLayer('./data/map/${z}/${x}/${y}.png', { maxZoom: 16 , tms: true});
map.addLayer(layer);

But what i get is a blank map. If i inspect whats going on with Firebug it seems that tiles are not even being requested (no other errors either). The map is set up right because i can see leaflet map widgets, i suspect its got something to do with the url of the layer. I copied that directly from my other OpenLayers app which works without a problem. Any ideas, or has anyone else stumbled across this?


texts: Why is there Big Red X in QGIS Composer?

texts:
I am working with QGIS version 2.01.  I created a print composer to export the map with legend, scale bar and title to an image or pdf.  When I opened my QGIS project again and went to my print composer there was a big red x which I can't see to get rid of.

Any ideas on went wrong here?


texts:
Is there a way to programmatically access a layer's visibility presets defined in the layer tree view widget (methods, properties, etc.)?

To avoid misunderstandings the part I want to access is the one highlighted in blue in the following screen capture (ex. preset "Test-01" ...)


texts: Generating a raster file using 1D numpy arrays

texts:
Explanation:
I currently have 3 1D numpy arrays, each having a length of over 100000 points. The first contains the X coordinates of an object, the second contains the Y coordinates, and the third one contains the speed of the object at the corresponding coordinates.The 3 arrays are in sequence from the start to end. I am trying to generate a raster file that will match another raster file that I produced using certain images.For simplicity's sake, the output raster size is 400 meters x 400 meters, and each cell has dimensions of 20 X 2. How exactly do I go about generating a raster using my 3 arrays ?

Ideally I would like to do something as simple as:

#xm = x-cooridnate array
#ym = Y-coordinate array
#speed = speed array

xmin,ymin,xmax,ymax = [xm.min(),ym.min(),xm.max(),ym.max()]
#nrows,ncols = np.shape(mining_durations)
xres = (xmax-xmin)/float(2)
yres = (ymax-ymin)/float(20)
geotransform=(xmin,xres,0,ymax,0, -yres)
output_raster = gdal.GetDriverByName('GTiff').Create(outputFileName,400/20, 400/2, 1 ,gdal.GDT_Float32)  # Open the file
output_raster.SetGeoTransform(geotransform)
output_raster.GetRasterBand(1).WriteArray(speed)
output_raster.FlushCache()

But this of course doesn't work.I know that many points will coincide in one cell, but I assume that it will automatically take the average of all the coinciding points.

Do I need to first create some sort of grid for the speed array?

Edit:
I have added sample data, as well as the code I am using.

Sample data: The coordinates are just float values. For example, the first column here is X, the second is Y and the last is some other variable (speed):

52039,57 -54169,66 7,68667678522605
52027,68 -54180,01 6,8999999999869
52039,55 -54169,70 2,12132034354421
52039,59 -54169,71 1,34164078652916
52039,61 -54169,64 0,670820393288984
52039,63 -54169,51 0,600000000013097
52039,69 -54169,49 0,750000000043656
52039,70 -54169,46 0,150000000030559
52039,67 -54169,58 0,75

Each sample sample/row of data is taken every 4 seconds, although this sometimes changes. The X and Y coordinates don't necessarily follow a fixed pattern. The code I am using is shown below.Note that this doesn't produce an output of 400x400 which I originally had in mind, the size is based on the coordinates.

def GenerateRaster(xdata, ydata, zdata):

    binWidth = 2.0
    binLength = 20.0
    xMin = min(xdata)
    xMax = max(xdata)
    yMin = min(ydata)
    yMax = max(ydata)

    statistic,x_edge,y_edge,binnumber = stats.binned_statistic_2d(xdata, ydata, zdata, 'mean',bins = [np.arange(xMin, xMax, binWidth), np.arange(yMin, yMax, binLength)] )
   # statistic = np.ma.masked_invalid(statistic)
    nrows,ncols = np.shape(statistic)

    xres = (xMax-xMin)/y_edge.shape
    yres = (yMax-yMin)/x_edge.shape

    geotransform=(xMin,xres,0,yMax,0, -yres)
    output_raster = gdal.GetDriverByName('GTiff').Create('Test_Raster.tif',ncols, nrows, 1 ,gdal.GDT_Float32)  # Open the file
    output_raster.SetGeoTransform(geotransform)  # Specify its coordinates
    output_raster.GetRasterBand(1).WriteArray(statistic)   # Writes array to the raster
    output_raster.FlushCache()

Although I do produce a raster file, I am not sure if this is correct or accurate. I noticed that the cell sizes are not exactly 20 x 2 which is strange: the dimensions are slightly longer.


texts: How to synchronously execute a queryTask

texts: How can a queryTask be executed synchronously? I'd like a feature set to be returned instead of deferred object.


texts: How to take a set of geohashes and determine which ones are neighbors?

texts:
I'm trying to figure out the optimal way to take a list of geohashes and return a sub-list that contains the geohashes that are direct neighbors or have direct neighbors in common (aka - creating groups of contiguous geohashes).

The only way I can think to do this is by looping through the geohashes, getting neighbors and adding them to a new list. Something like the psudo code below:

startingGeoHashes = ["9qczj5t",
                     "9qczj5s", # directly west of the above
                     "9qczj5e", # directly west of the above
                     "9qczjhh" # not a direct neighbor
]
neighborGroup = [] #This would actually be a list of lists to acccomadate multiple distinct groups
foreach startingGeoHash in startingGeohashes
   if any of startingGeoHashes in getGeoHashNeighbors(startingGeoHash):
      add both to neighborGroups if they aren't already in there

ideally the code above would return two groups:

9qczj5s,9qczj5t,9qczj5e
and
9qczjhh

But I'm wondering if there is a set based approach that could return the subset of direct neighbors without looping through every record and getting the neighbors to compare to everything else in the list.

Is this possible?

Language is not important but I will be implementing this in go if that makes any difference.


texts:
I have two DEMs of an open-pit mine taken a few weeks apart. I would like to estimate the volume that has been excavated during that interval, ie the difference of volume between the two DEMs.

Can I do that with GRASS? If so how?

If not, do you have any suggestion on how to do it? (I would prefer using an open source software...)


texts: Performing classification using composite of Sentinel 2 instead Landsat8 in Google Earth Engine?

texts:
I am looking fo the equivalent of the following code:

    // Load the Landsat 8 scaled radiance image collection.
var landsatCollection = ee.ImageCollection('LANDSAT/LC08/C01/T1')
    .filterDate('2017-01-01', '2017-12-31');

// Make a cloud-free composite.
var composite = ee.Algorithms.Landsat.simpleComposite({
  collection: landsatCollection,
  asFloat: true
});

I know that the equivalent for this i just loading the landsat-2 collection. But at the end I want to use something like this (this is where it fails for sentinel-2, but works for landsat-8):

    // Merge the three geometry layers into a single FeatureCollection.
//classes i uses for the classification
var newfc = urban.merge(vegetation).merge(water);

// Use these bands for classification.
var bands = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7'];
// The name of the property on the points storing the class label.
var classProperty = 'landcover';

// Sample the composite to generate training data.  Note that the
// class label is stored in the 'landcover' property.
var training = composite.select(bands).sampleRegions({
  collection: newfc,
  properties: [classProperty],
  scale: 30
});

// Train a CART classifier.
var classifier = ee.Classifier.cart().train({
  features: training,
  classProperty: classProperty,
});

// Print some info about the classifier (specific to CART).
print('CART, explained', classifier.explain());

// Classify the composite.
var classified = composite.classify(classifier);
Map.centerObject(newfc);
Map.addLayer(classified, {min: 0, max: 2, palette: ['red', 'green', 'blue']});


texts: How to create multi criteria cost surface - Slope reclassification?

texts:
I´m working on thesis  aimed at investigating of  possible routes and paths between Mayans cities in Guatemala. My thesis should help archaeologist to locate and reveal old paths in jungle and speed up all works in field with some scientific background.

I will use 30m DEM and anisotropic LCP approach with Tobler´s hiking equation.
Problem:

I need to create multi criteria cost surface like here:
http://resources.arcgis.com/en/help/main/10.2/index.html#//009z0000001z000000

I don´t know how to create cost raster with correct parameters. I mean I don´t know how to reclassify slope and add costs to it. I tried to use slope

1)how it was created (just derived from DEM using "Slope" function

2)reclassify slope from 1) and, of course, results are different.

I can´t find correct range for intervals and its costs. I need the same, exact range like 0-5° (cost 1), 5-8° (cost 2),... based on scientfic data.

And is this correct way of solving LCP problem? Should I create
multi criteria cost surface and use it as "in_cost_raster" and Tobler equation as vertical factor in Path Distance Tool or use slope (without reclassification) as "in_cost_raster", calculate Path Distance and then add cost of Land use?


texts: How to recover spatial coordinates that correspond to specific value of a raster file?

texts: I want to recover the spatial locations that correspond to specific one value of my raster file. Is there a  way to do that?


texts:
I have a road network with 8.5 million edges but when I run pgr_createTopology it processes edges at a rate of about 1000 edges per second up until 360,000 edges where it slows down and eventually fails with the following error.

server closed the connection unexpectedly
  This probably means the server terminated abnormally
  before or while processing the request.
The connection to the server was lost. Attempting reset: Failed.

I am running postgres 10.1 postgis 2.4 pgrouting 2.5. What would like cause it to fail like this? Or is there any way that I could process the network in 100,000 edge chunks?


texts: I have TIFF raster layers in QGIS project. I need to 'erase' them ( like in MS Paint ), clean them of any smudges from low scan. I do not see a proper 'eraser' tool like that in QGIS. How can I remove unwanted contaminations of layer?


texts: Determining average rainfall for cumulative watershed, where only sub-watersheds defined?

texts:
I have a polygon layer with sub-catchments (sub-watersheds) and a grid of rainfall.  I have used the Zonal Stats as a table to calculate the average rainfall for each sub-catchment, but I now need to get the cumulative average rainfall for each sub-catchment as I move down through the catchment.

I'm not sure how to prepare a layer with the "accumulated" sub-catchments, and then run the Zonal statistics (which would then be on overlapping polygons which I don't think it likes!)


texts: Find network distance for multiple origins and destinations?

texts:
I have a layer of points A (around 30), a layer of points B (around 2000), all recorded as in longitude and latitude. Need to calculate the distance between each pair (Ai, Bj) based on the transportation network. (only distance is needed, not direction, when I say distance I meant shortest distance, but we have a high tolerance to the accuracy, imagine as the result popped out from Google map direction)

I have tried everything, the closest thing I find here is the Network Analyst module in ArcGIS, but the result of Closest Facility (How to find network distances between two sets of origins and destinations?) is not what I wanted. I need a full matrix instead of just nearest facility.
Also I tried Google Distance Matrix API, that did work for me, but the usage limit is 2500 elements per day.

I meant to write a little python code looping the Routing/Direction function in ArcGIS, but i'm really new to python, if you can give me a brief idea how to implement this, that would be much much helpful.


texts: Loading Coordinates without shp file but prj, cpg, hdf, xml

texts:
I want to import the grids for pixel tiles of a MODIS data set and received the data for that from someone else. However, it does not contain a .shp file but only .prj, .cpg., .hdf and .xml.
Is there a way to illustrate those tiles with the data I have?
In the end I need to import them into Google Earth Engine but if there's a way to load it in ArcGIS Pro or QGIS, it already helps a lot!
If it is not possible to load the grid like that, what further information would be needed?


texts: Recommended utility for editing of GEOJSON files

texts:
I have a geojson file with the boundaries and some very basic statistics (population) for all the counties in US. States are enumerated (1,2,...), not named.

I also have a csv file with a comprehensive demographic and health statistics for each county (county is a row, statistics are columns).

I would like to add these values from the csv file as properties, for the county.

Is there a good, open source tool you would recommend for doing this programmatically, and in bulk.
I was looking into Python's geojson.


texts: Setting ArcMap Data Source Refresh Interval?

texts:
In ArcGIS Online there is an option to set a refresh data interval.

Is there a way to accomplish the same within ArcMap?

Specifically, my data source is inside an Enterprise GDB (Oracle).

I can pan the map around and that refreshes the data however I am looking for a way to not interact with the map and have it refresh on its own.

It seems this is possible in ArcGIS Pro 2.0+ but I haven't found a solution yet for ArcMap except for creating a simple python add-in which I would prefer to avoid because it halts other operations inside the ArcMap session.


texts: Delaunay triangulation returning empty geometry in PyQGIS

texts:
I would like to compute the Delaunay Triangulation for polylines. But sometimes it works and sometimes it does not, I tried with hand drawn lines.
lines=QgsProject.instance().mapLayersByName("lines")[0]
for line in lines.getFeatures():
    print(line.geometry().delaunayTriangulation(0.01, False).asMultiPolygon())

I tried a cross product collinearity test but I did not have a product equal to 0 and I still have empty geometry

from PyQt5 import *
from math import *
def modulo(a,b):
    if b!=0:
        return a-(int(a/b)*b)
    else:
        return -1
lines=QgsProject.instance().mapLayersByName("lines")[0]

for line in lines.getFeatures():
     n = len(list(line.geometry().vertices()))
     delaunay = line.geometry().delaunayTriangulation(0.01, False)
     i = str(line.id())
     if n <= 2:
         print("Line ID: " + i + " - not enough vertices for Delaunay triangulation")
     elif delaunay.isEmpty():
         print("Line ID: " + i + " - possible collinear vertices")
     else:
         print(delaunay.asGeometryCollection())

m=len(line.geometry().asPolyline())
pts=line.geometry().asPolyline()
for i in range(0,m):
    i1=modulo((i+1),m)
    i2=modulo((i+2),m)
    x1=pts[i1][0]-pts[i][0]
    y1=pts[i1][1]-pts[i][1]
    x2=pts[i2][0]-pts[i1][0]
    y2=pts[i2][1]-pts[i1][1]
    print(x1*y2-y1*x2)#if collinear then vector product=0
    #print(x1*x2+y1*y2)#scalar product


texts: Changing slopes derived from DEM using ArcGIS Desktop?

texts:
I have a DEM for a region over which I want to simulate an ecological process.  I would like to run several simulations for the same region but I want to alter the topography for the simulations.  Most important, I would like to alter the slope (and probably aspect) in various ways.  I am working in ArcGIS Desktop and I extracted the slope from the DEM but I have not found any way to manipulate the slope or even access the slope values with their coordinates.  Running the process requires several other layers so I would like to remain within the parameters (coordinates, datum, ellipsoid, etc.) of the original DEM.

Is there a way I can alter the slope output, which is currently a raster?

Could I, for example, make the whole region 30% slope with a south aspect?

Is there some other way to generate custom slope files that will be consistent with the other layers?

I am using an ArcGRID DEM from the NED.

Perhaps a different starting point is more appropriate, such as a contour map?


texts: Querying slope data using raster calculator: why are all results true?

texts:
I am using ArcGIS 10.2.2. From a DEM, I have created a new raster showing slope using the slope tool. The color scheme was automatically classified into 9 tiers. What I want to do is query this raster to identify zones above and below 9 degrees to get a Boolean, true/false result. When I use the raster calculator to do this (by inputting "slope" >= 9) this does happen - but the entire area is defined as 1 (true) with no areas defined as false. Other times, the result is a raster with a color ramp representing values from high (1) to low (0) rather than classes.

Is there any way to get a true/false outcome that doesn't count everything as true?

What I intend to do is convert the zones with slope of over 9 degrees to vector polygons and then execute Euclidean distance from these polygons to define areas further away from the slopes as being less under their influence.

I'm somewhat new to GIS. Also, my knowledge of programming languages is non-existent at the minute.

All layers in the data frame are projected using the TM65_Irish_Grid PCS.
No z factor was included when deriving the slope from the DTM.
The slope is measured in degrees.
The slope values range from 0 to 41.29286575 (it has been rounded up a bit in the ToC; under properties the max. is longer - 41.29286575317383).

Screenshot of the DTM extract/
Screenshot of the slope raster/
Screenshot of raster calculator results


texts: Boundary of merged parcel feature does not disappear after merge

texts:
In QGIS, when I merge several parcel features using Merge Selected Features, some of the boundary of parcels still remains.

Steps I used:

Enabled snapping options.
Enabled Toggle Editing.
Select Add Feature.
Created some parcel features.
Clicked on Save Layer Edits.
Goto View>Select>Select Feature(s).
Selected multiple parcel features that I created.
Selected Merge Selected Features.

Now when I save the merge, some boundary still remains.

I am using QGIS version 2.18.4

Some images:

i. Created parcel features with snapping options enabled.

ii. Selected parcel features for merge.

iii. Boundary still remains after merge.


texts: Why does QGIS show too many labels per feature?

texts:
When i check to label polygon feature, one out of a total five polygons shows like 10 labels on top of it. That particular polygon is only one polygon, represented by one row in the table, so i dont see why it shows 10 labels on top of it.....?

I havent selected any advanced label options, that could have maybe created a conflict.

Any idea anyone?

I know that in ArcGIS there is an option; show only one label per feature - BUT there no setting like that in qgis, that i have found.

a screen of the problem


texts: How exactly is the centroid of polygons calculated?

texts: I would like to know how exactly QGIS Geometry tools -> polygon centroid is calculating the point coordinates. For example does it divide the boarder of the polygon into small points, take two pairs of two points (length and width) with maximum distance and take the crossing of these two lines as the centroid?


texts: Selecting rows in a layer using arcpy.SearchCursor

texts:
I'm working on a project to identify points that fall along lengths of road. The roads are stored a single featureclass within a file geodatabase, as are the points. I'm using ArcGIS 10.3.

I'm trying to use an arcpy search cursor, configured with a for statement to iterate through my roads dataset. Within the for loop I'm using the function row.getValue to return a 3 digit short integer stored within a field called 'FacilityID' within the attributes of the roads layer. I wish to then use this value to select the appropriate row within the roads layer. Further down the line I would then use this selected row as the basis of a select by location query. At the moment my code doesn't iterate through selections. I can get it to print all the 3 digit integers within the 'FacilityID' field, but it ceases to work if I use the arcpy.SelectLayerByAttribute_management tool.

As I feel I'm not explaining all of this very well, I will include a cut down version of my code:

import arcpy, os
from arcpy import env
arcpy.env.workspace = "C:\Temp\Test.gdb"

env.overwriteOutput = True

ROADS = "Test_Roads"

FIELD = "FacilityID"

roadcursor = arcpy.SearchCursor(roads_lyr)

for row in roadcursor:
    ROADID = (row.getValue(FIELD))
    print ROADID, type(ROADID)

    expression1 = "\"FacilityID\" = " + str(ROADID)
    arcpy.SelectLayerByAttribute_management(roads_lyr, "NEW_SELECTION", expression1)

    arcpy.SelectLayerByAttribute_management(roads_lyr, "REMOVE_FROM_SELECTION", expression1)

I've changed the selection type to 'clear selection' in the second, and get the same result.


texts: Add external WFS layer to geoserver using REST API

texts:
I am trying to add an external WFS as a datasource to my geoserver via REST API. I already added a new workspace and WFS datasource.

Now I want to publish a layer from that WFS on geoserver. The same works for WMS, but seems to be a bit different and poorly documented for WFS services.

I tried to request

  PUT http://myhost.com/geoserver/rest/workspaces/my_workspace/datastores/my_datastore/layers

the following content as proposed in this post:

<layer>
  <name>my_layer</name>
  <type>VECTOR</type>
  <defaultStyle>
    <name>point</name>
  </defaultStyle>
  <resource class="featureType">
    <name>myFeature</name>
  </resource>
</layer>

This doesn't work out for me, I get a 500: Cannot cast org.geoserver.catalog.impl.LayerInfoImpl to org.geoserver.catalog.DataStoreInfo.

Now I am not sure, if it is required to create a FeatureType for the layer in advance. I added the layer manually via geoserver UI where it automatically detects the feature type from the WFS layer. I would like the same behaviour via the REST API.

Can anyone lead me towards the correct approach?

EDIT:

Maybe that will work, once I have created the featuretype. Just tried to get the featureType description from my manually created layer to adapt it. But there is a lot information I would have to put into the featuretype.xml manually just from the original WFS. I tried to insert the DescribeFeatureType output from the original WFS into this request:

  POST /geoserver/rest/workspaces/my_workspace/datastores/my_datastore/featuretypes

But it complains about the syntax because the xsd:schema tags

<xsd:schema elementFormDefault="qualified" targetNamespace="http://..."><xsd:import namespace="http://www.opengis.net/gml" schemaLocation=".../schemas/gml/3.1.1/base/gml.xsd"/>

Not sure, if this is the expected FeatureType syntax?


texts: Creating Python script to merge broken polylines from Microstation CAD file in ArcGIS?

texts: I have tried using the dissolve tool with unsplit_lines option but it does not merge the lines ... maybe there is a step I am missing?


texts:
I work with infectious disease. Right now i want to create a map showing a prevalence of infected people during an outbreak of the disease.
So I have coordinate data of infected/uninfected people live in the outbreak area.
Because the study area is a small village, there's no big road or boundary that we can use to calculate prevalence. So i created a grid of X/Y : 0.002 (i think it's 200 meter). Then I count GPS point of  infected/uninfected using point in polygon option in QGIS. Then I calculate prevalence in the grid (Infected/total population).
So in the end i get the prevalence in each grid but when i shown in the map. it's not very good because the prevalence showing as square shape.

So i come up with another solution which i'm not sure it's right or not. I created center of each grid and the center contain Prevalence value. I then make interpolation with those center. So now i can escape the square issue but because the interpolation use center point to interpolate, the prevalence shown is not 100% correct to the real location of the infected people.

Is there anyone have any suggestion for me. I'm still in the very beginning level of these mapping stuff.

I want my map to look like this:


texts: Handling QGIS Data Defined Label Position on PostGIS Layer Null Defaulting to 0?

texts:
I just upgraded to qgis 2.2 from 2.0.  In 2.2 there was a minor change to data defined label positioning at least when used with postgis layers.

In Version 2.2 the labeling engine appears to be treating null values in labelX and labelY as 0 instead of null.  It is also doing this in the feature form, but not when I view the Attribute table.  The result of this is that when I have data defined label positioning turned on and I save a feature without moving the label the label is rendered at 0,0 instead of at the feature centroid.

In postgis I have my label coordinates defined as real and qgis is treating them as a double float4.  After looking at another thread I loaded the same label coordinates as first a String, and second as an Integer.  My labels behaved properly when I had the coordinates defined as strings.

Is there any way to fix this problem by changing a setting within QGIS?

My options currently appear to be either changing all of my label data types to strings or writing trigger functions to populate my label position fields with the centroid of each polygon.


texts:
I have a DEM (.asc) file with a resolution of 1000 meters per pixel. Therefore the file has a lot of columns and rows.
I would like to chance the resolution per pixel from 1000 to 3000 meters per pixel. So i would like to resample it.

which gdal command do i use for this?
i know it's gdalwarp but how is the exact command line with the extra options?

if you know how to do it in QGIS I'm also happy with it.


texts: How to correlate NDVI and Land Surface Temperature in QGIS?

texts:
I'm looking for ways to correlate the NDVI values and the land surface temperature values that I have on QGIS. I'm kind of new at using these kinds of programs. I'm working on a thesis and for my results I want to know the correlation between NDVI and LST.
How can I do this in QGIS or how can I do it in any other program?
Keep in mind I don't know much about GIS in general.


texts: Adding geometry to rows with missing X,Y inArcGIS

texts:
I have a table with X, Y data that contains many rows with still missing coordinates. Now I would like to edit it in ArcGIS 10 and simply add the geometry for selected lines after identifying the correct location on a base map.

This seems so simple and basic, but I can't get it to work. Always, adding a point in ArcGIS will just create a new row. With more than 1000 yet unidentified spots, this is really cumbersome.


texts:
I am currently using default simple symbol capabilities or dingbats fonts for point objects.
Is there maybe a SVG symbol library for QGIS that is a bit more extensive than that which is included?
I've found this site http://mapicons.nicolasmollet.com which has some icons but its not what im looking for. These are actually more suitable for google maps.


texts: PgRouting : Shortest path starting from the closest vertex not the closest node

texts:
I am currently working on a routing solution with PgRouting (Dijkstra's algorithm), but the shortest path it creates always start from the nearest nodes generated from the edge table, and sometimes the distance between the node and the point is too long(4km or more) and it makes my calculations wrong.
Is there a way to calculate the shortest path from the closest line vertex to my starting point?
I included a picture to clarify more my problem


texts: How do you show agreement between polygons based on hand-drawn maps?

texts:
I'm somewhat new to GIS, coming from a linguistic background, and now I'm trying to make QGIS do something that I want from that background. I have some hand drawn maps with certain areas of the Netherlands circled. I want to see how much agreement these drawings show. I've dutifully entered these maps in a polygon shapefile, but now I'm stuck.
Now I know (from a brief description in an article) that in ArcGIS this should be relatively easy using Union and Frequency, but I cannot seem to find out how to do this in QGIS. Is it at all possible to make a heat map based on polygons, or do I need to use some other process?


texts: ArcObjects 10.1 IDatabaseCompact lock file

texts:
I'm trying to automate one of task by creating code that will execute the IDatabaseCompact.Compact method to a GDB file. Everything seems to work fine but there is always a .lock file generated after the Compact method is called. I am planning to do something else after the Compact process but it errors out because of the lock file.

Is there any way to close/delete the lock file?

BTW, this is how my process works:

first, I bind and initialize my system using:

ESRI.ArcGIS.RuntimeManager.Bind(ESRI.ArcGIS.ProductCode.EngineOrDesktop)
_ao = New ESRI.ArcGIS.esriSystem.AoInitialize()   _ao.Initialize(ESRI.ArcGIS.esriSystem.esriLicenseProductCode.esriLicenseProductCodeBasic)

I then call my compact method:

Public Sub CompressFileGDB(gdbFilePath As String)
    Dim dbCompact As IDatabaseCompact
    Dim ws As IWorkspace = FileGdbWorkspaceFromPath(gdbFilePath)

    If (TypeOf ws Is IDatabaseCompact) Then
        dbCompact = CType(ws, IDatabaseCompact)
        If (dbCompact.CanCompact) Then
            Dim didFinish As Boolean = True
            Try
                dbCompact.Compact()
            Catch ex As Exception
                didFinish = False
            Finally
                dbCompact = Nothing
                ws = Nothing
                _ao.Shutdown()

                RaiseEvent OnProcessDone(didFinish)
            End Try
        End If
    End If
End Sub

The compact method works, there is a significant change in the GDB file size. I also tried to set my attributes to Nothing after the process but it doesn't seem to help that much.


texts:
I have a folder on my local drive with 1,500+ JSON files (lines) that I need converted to shapefiles. I've tried creating a batch model in Python and ModelBuilder without success.

I'm using ArcGIS 10.2 and PyScripter.

Via KHibma: (this is not yet working for me, maybe someone can chime in)

import arcpy
import os

dirOfJSON = 'c:/temp/jsons'
outputDIrectory = 'c:/temp/jsonSHPs'

arcpy.env.workspace = dirOfJSON

for f in arcpy.ListFiles('.json'):
    j = os.path.join(dirOfJSON, f)
    outputShp = os.path.join(outputDIrectory, os.path.splitext(f)[0] + ".shp")
    print("Processing : {}".format(outputShp))
    arcpy.JSONToFeatures_conversion(j, outputShp)


texts: How to display a lot of images with OpenLayers?

texts:
I would like to display a lot of images at same time. The images are georeferenced (they have extents), with some labels on top of it which, on click, open popup. I cannot say how many exactly because it will depend on user choices.

And I want the user to choose the information he sees. There are some checkboxes he can check or not and it will display more or less images depending on his choices.

My problem is that I have to deal with the "issue" (which actually is not really an issue) described in OpenLayers wiki: What is the maximum amount of layers I can have in my OpenLayers Map?

I was wondering if there is a way to have more than 75 images. Because, the system actually works, except when I want to display more than 75 images: Popups and labels appear below images...

I tried a long time to play with z-index, without success though.


texts: Creating raster from several polygon layers using ArcGIS Desktop?

texts:
I have the following map:

The four layers have the following attribute tables:

SF_UrbanTreeCanopy_84

SF_Building_Footprints_84

GrassBeach

Water

I want to create a raster with the following extent:

Top: 37.822265 dd
Below: 37.708089 dd
Left: -122.522052 dd
Right: -122.351798 dd

That is the extent of the layer Water and with 0.00002 decimal degrees cell size. So the raster should have 5709 rows and 8513 columns. Now:

Values of the cells that are within SF_UrbanTreeCanopy_84 should be
1
Values of the cells that are within SF_Building_Footprints_84 should be
2
Values of the cells that are within GrassBeach should be 3
Values of the cells that are within Water should be 4
Values of other cells should be 5

That is I want to create a ground truth or a classification map.

Is there any way to do such in ArcGIS 10.3?


texts: Limit number of displayed points based on zoom level in GeoServer

texts: My goal is to display large datasets with GeoServer using zoom level limits. My idea is that I only display the parts of the datasets in a certain zoom level, but after each zoom, I want to increase the number of displayed points. Is there any solution using GeoServer? My client will be Leaflet.js.


texts: Layer 'load' event not triggering in CARTO.js

texts:
I amusing CartoDB 3.11.23 version. I have registered layer 'loading' event and 'load' event.
It fires both for the first time. After I zoom in/out the map it fires only 'loading' but not 'load' event. Below is my code -
cartodb.createLayer(this.map, {
            user_name : cartoDBUser,
            type : 'cartodb',
            sublayers : [ {
                sql : <select query>,
                cartocss : '#dummy {polygon-opacity: 0.0;}'
            } ]
        }, {
            https : true,
            no_cdn: this.noCDN || false
        })
        .addTo(this.map)
        .on('done', onDone)
        .on('error', onError);

var onError = function() {
            this.alert('Map Error', 'There was a fatal error loading the map.');
        }, onDone = function(layer) {
            me.cartoLayer = layer;

            layer.bind('loading', function() {
                me.loader.show(); // THIS SHOWS THE LOADING SPINNER
            });

            layer.bind('load', function() {
                me.loader.hide(); // THIS HIDES THE SPINNER. BUT THIS IS NOT TRIGGERING
            });

            me.clearQueryMapArea();

            // Now add the viz
            me.addVisualization(viz);
        };


texts: Assigning values to a grid cell based on number of intersecting features

texts: I need to create a grid that will include a single value for each cell which represents the total number of "stressors" intersecting that cell. So, the final product will be a 1km grid with cells whose values will range from 0 to n (total intersections). Does anyone know of the best way to accomplish this?


texts:
My QGIS is not showing Hindi Fonts. It happened after I formatted my Windows 10. To fix it, I installed Hindi Fonts, Google Input tools, etc. but didn't work. The Hindi fonts work in other software but not in QGIS as shown in the image below.
QGIS version 2.18.3 to 2.18.14

Is there any fix for this?


texts: Why does this HDF4 file converted to .tif with gdal come out at a coarser resolution?

texts:
I'm currently working with large-domain MODIS data. These come in multiple separate HDF4 files which I stitch together into a virtual data set (VRT), reproject into EPSG:4326, sub-set to the domain of interest and convert to .tif. This seems to work as expected.
When I repeat this procedure with a VRT that consists of only a single MODIS tile, the resulting .tif file for some reason comes out at about half the original resolution. Specifying a resolution for the VRT with e.g. -resolution highest does not change this. What am I missing here?
MWE
GDAL 3.0.4 with HDF4 support
Input file (11 MB): https://github.com/wknoben/gis_se/blob/main/MCD12Q1.A2015001.h10v03.006.2018146022714.hdf
gdalbuildvrt modis.vrt MCD12Q1.A2015001.h10v03.006.2018146022714.hdf -sd 1 -resolution highest
gdalwarp -of VRT -t_srs "EPSG:4326" modis.vrt modis_4326.vrt
gdal_translate -projwin -116.55 51.74 -115.52 50.95 modis_4326.vrt modis_4326_cut.vrt
gdal_translate -co "COMPRESS=DEFLATE" modis_4326_cut.vrt modis_4326_cut.tif

Source HDF and resulting tif


texts: How can I carry out a reliable coordinate transformation for point data?

texts:
I have ArcGIS 10.0 with ArcGIS Server (AGS) 10.0. Due to a known bug in the AGS geometry service we cannot carry out coordinate transformations with AGS 10.0. Our Internet client therefore does not have a tool which converts point coordinates between out ETRS89 and the older( but still in use) Gauß-Kruger zone 3 System. A user asked me today to covert some coordinate pairs, and after doing it with arcgis I checked a couple of Online Services which claim to do the same thing. I found that a got 3 different results for the same coordinate pair!

My Test in ArcMap 10 was to simply click with the info Tool in the map projected in ETRS89 on an area in the aerial photo and read off the native coordinates. I then created a new data frame which had the same aerial photos transformed into Gauß-Kruger Zone 3 and clicked on the same point. The following readings were produced:

ArcGIS

Gauß Kruger Zone 3 Dataframe
3570494.458
5877100.795

ETRS89 UTM Dataframe
570395.131
5875189.237

I gave the Gauß-Kruger coordinates into the following 2 sites to see what they produced and got 2 different results.

Berg- und Wanderportal

The World Coordinate Converter

Since I know that AGS Geomtry service transforms incorrectly, how should I test whether the coordinates are correct in ArcMap?


texts: Setting attribute table cell heights in print composer in QGIS

texts: Is it possible (and how) to set the attribute table cell heights in print composer in QGIS?


texts:
I have a hillshade that I would like to use as a background in an ArcPad 10 map.  I get everything set how I want in in ArcMap and when I export the map to use on my Trimble the hillshade doesn't show up right:

The locations of each image are not exact but you get the idea.  I've tried forcing it to TIFF, JPEG2000, and using the native format.  Nothing seems to work.  Initially the raster was stored in a geodatabase, so I tried converting it to a TIFF and exporting it from there - when I do that the image doesn't even load in ArcPad.

I've now also attempted to load the raster directly onto the Trimble and then load it as a layer in an empty APM rather than use the one created when I used the 'Get Data For ArcPad' button on the Arcpad Data Manager Toolbar in ArcMap.  Still nothing seems to work, the raster just will not display.  The raster isn't huge, about 110 MB, 32 bit, unsigned integer, and I've checked the projection about 40 times, everything is in the same projection.

What am I missing about how to get background rasters to display correctly in ArcPad?
Should I be using a different file format?


texts: Selecting by attributes and making new shapefile using Python ogr only?

texts:
I have a problem making my selected features into a new shapefile with the same projection. I am new to python, any ideas?

Below is my sample script.

from osgeo import ogr
import os, sys

driver = ogr.GetDriverByName ("ESRI Shapefile")

fnVec= raw_input("Enter Filename:")
ofn=raw_input("Output Name:")
ctype=raw_input("Crime Type:")
cyear=int(raw_input("Year of incidence:"))

datasource = driver.Open (fnVec,0 )

inLayer = datasource.GetLayer()

#Execute queries(SQL)
CrimeSQL ="SELECT * FROM Crimes WHERE Crime_Type='%s' AND Year='%s'" %(ctype,cyear)
Crime =datasource.ExecuteSQL(CrimeSQL)

Crime_feature = Crime.GetNextFeature()
while Crime_feature:
    #print Crime_feature.GetField('Count')
    Crime_feature = Crime.GetNextFeature()

My problem now is to make the selected points into a new shapefile. I made new scripts, but I cannot view the selected points. Here is my sample code:

from osgeo import ogr, osr
import os, sys

# Get the input Layer
inDriver = ogr.GetDriverByName("ESRI Shapefile")
fnVec= raw_input("Enter Filename:")
ofn=raw_input("Output Name:")
ctype=raw_input("Crime Type:")
cyear=int(raw_input("Year of incidence:"))
inDataSource = inDriver.Open(fnVec, 0)
inLayer = inDataSource.GetLayer()
spatialRef = inLayer.GetSpatialRef()

# input SpatialReference
feature = inLayer.GetNextFeature()
geom = feature.GetGeometryRef()
inSpatialRef = geom.GetSpatialReference()

#Execute SQL
CrimeSQL = "SELECT * FROM Crimes WHERE Crime_Type='%s' AND Year='%s'" %(ctype,cyear)
Filter = inDataSource.ExecuteSQL(CrimeSQL)

# Create the output LayerS
outShapefile = os.path.join( os.path.split( fnVec )[0], ofn )
outDriver = ogr.GetDriverByName("ESRI Shapefile")

# Remove output shapefile if it already exists
if os.path.exists(ofn):
    outDriver.DeleteDataSource(ofn)

# Create the output shapefile
outDataSource = outDriver.CreateDataSource(outShapefile)
out_lyr_name = os.path.splitext( os.path.split( outShapefile )[1] )[0]
outLayer = outDataSource.CreateLayer( out_lyr_name, geom_type=ogr.wkbPoint )

# Add input Layer Fields to the output Layer if it is the one we want
inLayerDefn = inLayer.GetLayerDefn()
for i in range(0, inLayerDefn.GetFieldCount()):
    fieldDefn = inLayerDefn.GetFieldDefn(i)
    fieldName = fieldDefn.GetName()
    outLayer.CreateField(fieldDefn)

# Get the output Layer's Feature Definition
outLayerDefn = outLayer.GetLayerDefn()

# Add features to the ouput Layer
for inFeature in inLayer:
    # Create output Feature
    outFeature = ogr.Feature(outLayerDefn)

    # Add field values from input Layer
    for i in range(0, outLayerDefn.GetFieldCount()):
        fieldDefn = outLayerDefn.GetFieldDefn(i)
        fieldName = fieldDefn.GetName()
        outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(),inFeature.GetField(i))
        outLayer.CreateFeature(outFeature)
print "Process Completed!"
# Close DataSources
inDataSource.Destroy()
outDataSource.Destroy()


texts: How to bulk import gpx files to QGIS and merge into a single shapefile?

texts:
I am looking to create my own 'heatmap' of all the bike rides I have recorded using the app Strava in QGIS.
I have over 450 individual gpx files and would like to merge the track points layers into one single shapefile so I can join them to a polygon layer and symbolise by join_count as seen in this example:

Some of my gpx files contain 'ride' in the file name and others 'run', I'd like to filter the files leaving only those with 'ride'.
How can I achieve this efficiently using QGIS?


texts: How to get URL download in Google Earth Engine?

texts:
In the below code I want to download PM image using getDownloadURL instead of Export.image.toDrive() but returns error. How can I solve this problem?

code link: https://code.earthengine.google.com/2ecd8a1ea1788def648f378eb0f7b6b7
var pm = ee.ImageCollection("ECMWF/CAMS/NRT")
.filterDate('2016-01-01','2021-01-01')
.filterBounds(table)
.select('particulate_matter_d_less_than_25_um_surface')
.mean().clip(table);

print(pm)

var pmDownload = pm.getDownloadURL({
    'image': pm,
    'region': table,
    'filePerBand': 'false',
    'crs': 'EPSG:4326'
});


texts:
I am in the process of designing a story map through ArcGIS Online, and I have recently found a way to embed a 3D Web Scene into the story map as a URL. This allows for me to display my 3D data effectively, but it takes away from my ability to set the extent/zoom to certain views when changing tabs within my story map journal. The map as the main stage content allows for different extents to be set to change the view depending on what you want to direct the user to, however since the web scene can only come in as a URL, it does not have the same level of detailing options for setting views.

I want my story journal to update locations based on different tabs similar to the Mountains of Fire Story Map. Does anyone have any idea how I can implement this?

Right now my theory is that they are using a completely different web scene for each of these different locations. Is there any way to check?


texts: Controlling map details shown using Leaflet and MapBox

texts:
I am building a map based chart using Leaflet (with MapBox). I would like to know if there is a way to control the level of details displayed on the map when zooming in. For example I do not want to see Interstate highway or street names.

This link explains how to do this using google maps.


texts:
I have Geoserver 2.8.2 running on Tomcat7 on an Ubuntu 14.04 virtual machine. I started adding a large number (off the top of my head, about 5000) of smallish (roughly 2 to 10 Mb) Geotiffs as layers using a Python script and the gsconfig module.

After saving a style I was working on, the web interface became unresponsive. I restarted Tomcat but still no response from either GeoServer's web interface or from Tomcat's one.

So I stopped Tomcat and restarted. The whole server became unresponsive to pings or login attempts. The syadmin restarted the virtual machine, with apparently no data loss. Tomcat is up and running but the web interfaces for it and GeoServer are still unresponsive. The logs don't show anything particular, but catalina.out does - I have more than 50Mb of logs looking like this:

11 Mar 11:41:41 INFO [org.geoserver] - Loaded coverage 'es502c_R5_C1_4_T_2795_sec', enabled
11 Mar 11:41:41 INFO [org.geoserver] - Loaded coverage 'es502c_R5_C1_4_T_2795_sec'
11 Mar 11:41:42 INFO [org.geoserver] - Loaded layer 'es502c_R5_C1_4_T_2795_sec'
11 Mar 11:41:42 INFO [org.geoserver] - Loaded store 'cz010c_R9_C1_3_T_2751_sec', enabled
11 Mar 11:41:42 INFO [org.geoserver] - Loaded coverage store 'cz010c_R9_C1_3_T_2751_sec'
11 Mar 11:41:42 INFO [org.geoserver] - Loaded coverage 'cz010c_R9_C1_3_T_2751_sec', enabled
11 Mar 11:41:42 INFO [org.geoserver] - Loaded coverage 'cz010c_R9_C1_3_T_2751_sec'

...and so on. And it is still growing. These are the layers I was adding, but my Python script is definitely not running.

I am reluctant to restart Tomcat again, to avoid freezing the server again.

What is going on? It looks a bit like Geoserver is trying to catch up on its logging, but could it be stuck in a long loop?


texts: Google maps android vector rendering : how do they make it that fast?

texts:
Most of the question is in the title.

If you are on android (I believe its the same on IOS) and use maps you will notice the quality of the maps app and how the vector maps are drawn in a very fluid manner.

From what I understand, they are delivering vector tiles wich are encoded using protocol buffers which is pretty easily achievable on my side.  However I was wondering how they managed to achieve real-time rendering of vector tiles in such a fluid manner.

I am currently using mapsforge to render vector based maps but the rendering times of tiles are pretty far from what I can see in google maps.

Any ideas on how I can approach such awesomness ? What kind of voodoo tricks are they using ?

Thanks for any insight !


texts:
How can I get QGIS project metadata information via PyQGIS?
QGIS version: 3.16.


texts: Leaflet - do I need the rastercoords plugin for an non-geographical map with CRS.Simple and local tiles?

texts:
I got a fictional map set up following this tutorial: https://leafletjs.com/examples/crs-simple/crs-simple.html

But when I try to switch the imageOverlay for map tiles (generated with gdal2tiles-leaflet via docker run --rm -v pwd:/data niene/gdal2tiles-leaflet -l -p raster -z 0-5 -w none /data/<image> /data/<tilesdir>), I only see grey emptiness. What do I need to change to make that work? Or should I better use the rastercoords plugin?

My code currently looks like this:

    var w = 4636;
    var h = 6000;

    //the stuff to be able to work with xy instead of latlon
    var yx = L.latLng;
    var xy = function(x, y) {
        if (L.Util.isArray(x)) {    // When doing xy([x, y]);
            return yx(x[1], x[0]);
        }
        return yx(y, x);  // When doing xy(x, y);
    };

    var bounds = [xy(0, 0), xy(w, h)];

    var background = L.tileLayer('tiles/{z}/{x}/{y}.png', {
        attribution: 'your Map',
        bounds: [[0, 0], [-6000, 4636]],    //the minus is apparently needed to avoid leaflet trying to fetch tiles like -3.png
        tileSize: 256,
        noWrap: true
    });

    var map = L.map('map', {
        layer: [background],
        crs: L.CRS.Simple,
        minZoom: -3
    });

    map.fitBounds(bounds).setMaxBounds(bounds);


texts: simple standalone pyqgis app gives SIGSEGV error on quit

texts:
I am completely new to QGIS and have just recently installed the osgeo4mac version with homebrew. I am still setting up and trying to run a very simple standalone python script.

The following script works fine:

from qgis.core import *

app = QgsApplication([], False)
app.initQgis()

app.exitQgis()

The output I receive suggest that the data providers and layer registry are properly loaded and removed.

I am experiencing problems with the following script:

from qgis.core import *
import time

app = QgsApplication([], False)
app.initQgis()

shp_file = 'test.shp'
layer = QgsVectorLayer(shp_file, 'vlayer', 'ogr')

print layer.isValid()
for field in layer.pendingFields():
    print field.name()

app.exitQgis()

time.sleep(5)
print 'helloworld'

I receive confirmation that the layer is valid, and the script prints out the proper field names in my test shapefile. However, 'python quits unexpectedly' because of a  SIGSEGV error.

The last two lines are for debugging purposes and reveal what I think is a surprising behaviour. This is some of the output:

src/core/qgsproviderregistry.cpp: 234: (clean) [0ms] cleanup:ows
src/core/qgsproviderregistry.cpp: 234: (clean) [0ms] cleanup:postgres
src/core/qgsproviderregistry.cpp: 234: (clean) [0ms] cleanup:spatialite
src/core/qgsproviderregistry.cpp: 234: (clean) [0ms] cleanup:virtual
src/core/qgsproviderregistry.cpp: 234: (clean) [0ms] cleanup:wcs
src/core/qgsproviderregistry.cpp: 234: (clean) [1ms] cleanup:wms
helloworld
src/providers/ogr/qgsogrconnpool.cpp: 36: (QgsOgrConnPool) [5015ms]Entering.
src/providers/ogr/qgsogrconnpool.cpp: 36: (QgsOgrConnPool) [0ms] Leaving.
[Finished in 7.1s with exit code -11]

Any Clue?


texts: How to convert "linear" polygons to lines using NetTopologySuite?

texts:
If I have set of polygons representing country rivers, how to convert them to their linear representation using NetTopologySuite. I did some digging but i was unable to find exact function like Create Centerlines in some other software (Free, but with usability limits).

I could just unpack polygons to linear geometry, but it would be better if it is cleaned up somehow. I am not sure how.


texts:
I'm looking to plot day/night on a Google map, for an arbitrary point in time. I'm familiar with generating map tiles; I'm just looking for an algorithm to tell me whether a particular point on the globe is currently in daylight or darkness, or to otherwise plot the curve of the day/night interface onto the map.

I've done some searching, but it's possible I don't even know enough about the problem domain here to know what terms to search for!

Any ideas? Doesn't have to be perfect -- basically, I'm comparing Flickr geolocation data of sunrise and sunset photos (and their "date taken" timestamps) with reality, and this is to help me visualise it.


texts:
In my data I have to remove some line feature classes, which are those sharing boundaries with two states i.e. New York and Philadelphia.

I have to remove all data of Philadelphia but some feature classes cross the border and fall in New York.

How can I remove those line featureclasses which are in Philadelphia but share a boundary with New York?

attaching screenshot, that might help.

!light pink colour is New Hampshire  and blue line are going toward white colour MA


texts: Is it possible to make "Project Bookmarks" the default when adding a bookmark in QGIS?

texts:
Whenever I add a new bookmark in QGIS (3.14), I have to change the Saved in dropdown from User bookmarks to Project bookmarks. Is it possible to make Project the default?
The User bookmark doesn't even make any sense to me, so I don't know why it is the default. Why would anyone even want to save a bookmark that is not related to some GIS data?


texts:
I've done a little bit of SQL query in MapBasic, but somehow get lost in this problem. I'm hoping that someone can help me out.

I have a table of polygons in 2 groups, F and E. F polygons are larger than E in size. Some E polygons are inside the boundary of some F polygons. Some are not.

How can I select the F polygons where there are E polygons inside its boundary?

It can be very simple, I guess...


texts: Determining if another marker (Point or Polyline) where marker Point dragged in Leaflet?

texts:
I'm trying to create interaction between 2 markers with Leaflet, I would like to be able to drag a marker onto another

Either between one point and one linestring, or between 2 points.

I added the draggable:true option when creating my markers :

var marker = L.marker(latlng, {draggable: true});

So I'm able to move it and catch the dragend event.

But how do I know if there is another marker where I dragend it ?

I guess I could try to find markers in a radius zone from where I release the marker. But this depends on the zoom level I am. If I'm at zoom level 1, two points are "close" if there is only 100 km between them. If I'm at zoom level 18, two points are "close" if there is only 10 meters between them...

Also, I would love to be able to drag a marker onto a LineString between 2 points.
For example :

Point A *-----------------* Point B
                ^
                |
                | Dragging Point C in the middle of the LineString
      Point C   *

But when I'm dragging the Point C marker on the LineString it doesn't fire the onmouseover event of the LineString.

How do I know if there is another marker (Point or Polyline) where I dragend a marker Point ?


texts: cross domain problem in geodjango for openlayers in django development server

texts:
I want to display geojson data on map in geodjango. Below is my code
template/map.html (javascript code). only the map is loading but not the geojson data. Whwn I used debugger the url is gettin geojson data but not displaying on map and I'm not getting any errors.

<script>

        var map,layer;
        var zoom=8;
        var lon=45;
        var lat=5;
        function init(){
             map = new OpenLayers.Map('map');
             layer = new OpenLayers.Layer.WMS( "OpenLayers WMS",
                    "http://vmap0.tiles.osgeo.org/wms/vmap0",
                    {layers: 'basic'} );
                      map.addLayer(layer);
                      map.setCenter(new OpenLayers.LonLat(lon, lat), zoom);

            var featurecollection = {{geojson}};
            <!-- var geojson_format = new OpenLayers.Format.GeoJSON(); -->
            var geojson_format = new OpenLayers.Format.GeoJSON({
                'internalProjection': new OpenLayers.Projection("EPSG:900913"),
                'externalProjection': new OpenLayers.Projection("EPSG:4326")
            });
            var vector_layer = new OpenLayers.Layer.Vector();

            vector_layer.addFeatures(geojson_format.parseFeature(featurecollection));
            vector_layer.addFeatures(featurecollection);

            map.addLayer(vector_layer);

        }

views.py

def geojson_provinces(request):
    conn = psycopg2.connect(dbname="geodjango",host='localhost',user='postgres', password='postgres', port=5433)
    dict_cur = conn.cursor(cursor_factory=psycopg2.extras.DictCursor)
    cur=conn.cursor()
    res = dict_cur.execute("SELECT ST_AsGeoJson(geom) AS json FROM pft_pft LIMIT 10;")

    points=dict_cur.fetchall()
    print (points)
    json_res  = []
    for row in points:
        json_obj  = row['json']
        print (json_obj)
        json_res.append(json_obj)
        print(json_res)
    #json = simplejson.dumps(data)

    dict_cur.close()
    data = {'geojson':  mark_safe(json_res)}
    return render(request, 'map.html', data)

urls.py

 url(r'^geojson/$', geojson_provinces, name='geojson_provinces'),

models.py

class PFT(models.Model):
    species=models.CharField(max_length=50)
    # GeoDjango-specific: a geometry field (MultiPolygonField)
    geom = models.PointField(srid=4326)
    objects=models.GeoManager()
    # Returns the string representation of the model.
    def __str__(self):              # __unicode__ on Python 2
        return '%s %s %s' % (self.species, self.geom.x, self.geom.y)


texts: Geoserver transaction send null date value

texts:
I have a timestamp column in my postgis table, I send data to postgis using geoserver and wmts. I need to send a null value to the timestamp table sometime but I dont know how to do it. Is it possible to send a null value in the transaction?
I have tried to send null or just leaving it blank, but neither seem to work.

Any ideas?

java.lang.IllegalArgumentException: Failed to parse dateTime NULL at:NULL

java.lang.IllegalArgumentException: Failed to parse dateTime undefined at:undefined


texts:
I'm looking for shapefiles that delineate ocean current boundaries: specifically one for the Pacific Eastern Boundary current. This is for masking marine productivity raster data onto the current feature in ESRI.

Anybody have experience with this sort of thing?


texts: I have a sentinel 2 L1C data and I am trying to calculate the Land Surface Water Index (LSWI) using the NIR and SWIR bands. As the Sentinel 2 data is having 2 SWIR bands (B11 and B12), I am confused about which one to use for LSWI calculation


texts: Why does QGIS generate Shapefiles without DBF?

texts:
I need to output shapefiles for archiving, which requires generating .shp, .shx and .dbf files.  When I save a layer as a shapefile it generates the first two, but not the dbf.  All the help I can find tells me how to import a dbf file but not export.

The dbf files will need to contain ID, X and Y fields.

Any ideas?


texts: How to incorporate QGIS custom widgets into plugin?

texts:
Developing a plugin with Qt Designer for QGIS 2.14.1-Essen on Windows 7 and having trouble inserting the QgsCollapsibleGroupBox widget. I receive the following error:

  ImportError: No module named qgscollapsiblegroupbox

I put from qgis.gui import QgsCollapsibleGroupBox into my main.py file but it still doesn't recognise it. Only post I found which mentioned this error is QGIS custom widgets with python plugins  but not exactly sure how to install this qgis_customwidgets.py file.

I would like this plugin to be transferred to other computers but would like to avoid having to install this file each time (if it is indeed required).

Is there another way to use this type of widget or perhaps a method to save the qgis_customwidgets.py into my plugin directory and just call the module?


texts: Returning area of interest image upon lat/long input?

texts:
I would like to know if it's possible to give a latitude and longitude (or bounding box) for a certain location and return an area of interest image, including a target property and a list of nearby schools. Me and my manager went through the Esri docs but couldn't find a way of doing it. If not, any suggestion for a different api to do that?

I would like to get something like this image below where the red point of interest is a property and the black ones are nearby schools:

I am using the esri rest api.


texts: Feature class error in Feature compare (Data Management) analysis

texts:
I have some troubles on performing a feature compare data analysis.
I created a new feature layer with different point variables, see picture. The two different colored points are from two different data sets (1982) and (2014). The attribute tables contain similar variables and my idea is to be able to test the validity of the smaller data set (2014) compared to the (1982) data set on certain parameters to use the more extensive data that it contains.
However, when I perform the feature class analysis it gives a 'feature class extents are different' error.
Within the attribute table, the field names are similar (soil characteristics such as LUTUM, M50 etc). (not in similar order though) However, within the 'Omit_Field' the parameters are different according to the original layer that they are from: 'bhorizont9999.LUTUM' and 'bhorizont5159.LUTUM' for example. I thought that this might be an indicator of that the field names should be changed. Also, it might be due to a different spatial reference?

My main question is really if you think that this is the correct analysis and how I can make sure that I compare the paired profiles that are 'grouped' in the picture.


texts: Why are my polygon features area being calculated as zero?

texts:
I am trying to measure the area of polygons within a larger polygon. However some of them, seemingly without reason are being calculated as having zero area. The highlighted regions are calculated as zero and the non highlighted features are being calculated correctly:

Any ideas on why this is happening? I took the buildings, used buffer to create the borders, and then dissolved them.


texts: Vector tile map packages only in ArcGIS Pro?

texts:
I recently came across the .vtpk packages and I want to ask whether the CreateVectorTilePackage_management is available only in ArcGIS Pro or we can access it from ArcMAP 10.4.1. (or directly from arcpy)?

I tried in my Desktop environment and it looks like it is not there.


texts: Making more than 1 layer tree in opengeo sdk?

texts:
I would like to have more than 2 levels of layers in my geoexplorer.

I use groups, inside of "gxp_layertree", but is it possible to get more groups e.g. a tree of 3 or 4 levels?


texts:
I am trying to get my GRASS plugin working and found that I have to set the path for it, but the path does not exist on my computer.

I tried uninstalling QGIS and reinstalling it and the problem was not fixed, still missing GRASS folder under Applications/QGIS3.10.app/Contents/MacOS.

I have done the homebrew steps and it says everything went as it should but GRASS is not where it should be and I can't find it anywhere else. I am trying to get contour lines from QGIS to Fusion360 as a 3D model, is there another plugin I can use?


texts: Mercator vs. Web Mercator for Large Scale Topographic Maps

texts:
I have 8 terabytes of imagery that is in either lon/lat or UTM projection that I will be using for the basis of large scale topographic maps.

I have been studying the differences between Mercator and Web Mercator, but I am still not sure what projection is more suitable for large scale (scales between 1:24k - 1:63k) topographic maps. I do understand that Mercator is conformal unlike Web Mercator which seems like an advantage.

I would like to project these datasets in UTM as I understand that is the recommended projection by United Nations Cartography Committee for topographic maps because it minimizes angular distortion at large scales (and at detail greater than 1:24k) - and because about half of my data is already in UTM.

At the same time I would like the ability to seamlessly zoom from a large geographic area (~1:250k up to 1:24k. That makes the UTM zones problematic as I will not have the capability of reprojecting adjacent UTM zones into a common UTM zone on-the-fly.


texts: Generate territory polygon from sets of points

texts:
I'm looking for a simple way to generate contiguous "territories" from sets of points. Say, points A, B and C are in territory 1, and D, E and F are in territory 2, then generate polygons that lie roughly halfway between the two sets. No need to take into account routing or anything like that.

In the browser with Leaflet would be great, but failing that, GDAL or QGIS.


texts:
I'm having trouble with an if/else statement in my for loop.

I have a raster Vmax with a constant value, another raster V which has a value of 0 in the beginning and will change during execution of the loop. Furthermore I load two text files with rain and ET-data. The idea is to calculate a recharge depending on the change of V which in turn depends on rain and ET. So something like that in a loop:

AET_t = PET*(V/Vmax)
V_tt = V + prec - AET_t
if V_tt > Vmax:
    recharge = V_tt - Vmax
    V = Vmax
else:
    recharge = 0
    V = V_tt

This is my code:

#import modules and packages
import arcpy
from arcpy.sa import *
import os.path

#set workspace
arcpy.env.workspace = r"C:\Users\kern\Documents\model\test_area"

#load input rasters
Vmax = Raster(r"\input.gdb\Vmax")
V = Raster(r"\input.gdb\V")

#check if Spatial Analyst extension is available
arcpy.CheckExtension("Spatial")
#retrieve license
arcpy.CheckOutExtension("Spatial")

#open rain and PET textfile
rain = open("rain.txt")
pet = open("Epot.txt")
#create list of daily rain and PET data values
rain_list = rain.readlines()
pet_list = pet.readlines()

from itertools import izip, count
for i, p, et in izip(count(),rain_list, pet_list):
    prec = CreateConstantRaster(p, "FLOAT", "5", Vmax)
    PET = CreateConstantRaster(et, "FLOAT", "5", Vmax)
    output_aet = r"C:\Users\kern\Documents\model\test_area\AET.gdb"
    outfile_aet = os.path.join (output_aet, "aet" + str(i))
    output_re = r"C:\Users\kern\Documents\model\test_area\recharge.gdb"
    outfile_re = os.path.join (output_re, "re" + str(i))
    AET_t = PET*(V/Vmax)
    AET_t.save(outfile_aet)
    prj_file = r"C:\Users\kern\Documents\model\model_input_5x5\a_aplis.prj"
    arcpy.DefineProjection_management(outfile_aet, prj_file)
    V_tt = V + prec - AET_t
    recharge = Con(V_tt > Vmax, V_tt - Vmax, 0)
    recharge.save(outfile_re)
    if recharge > 0:
        V = Vmax
    else:
        V = V_tt

rain.close()
pet.close()

And this is the error:

File "C:\Users\kirn\Documents\MSc_Arbeit\model\test_area\WaterBalance.py", line 46, in
    if recharge > 0:
ValueError: The truth value of a raster is ambiguous. Invalid use of raster with Boolean operator or function. Check the use of parentheses where applicable.

It seems the problem is the V = Vmax and V = V_tt. But how can I assign Vmax and V_tt respectively to V in another way?
I found this https://geonet.esri.com/thread/79949 but couldn't really make head of it.


texts: Can I "undo" buffering of a point layer in QGIS?

texts:
I have point layer, where each point is one tree, about 1300 trees in total. I made buffers around each tree, the buffer being the height of each tree in QGIS 2.0.1.
My next step was to use Spatial Query to see if the buffers (the height of each tree) reached a road. The resault was that 60 buffers/trees reached the road.

I now want to show the location of each tree on a map, and I believe that the map would easier to read if it contained points instead of buffers.

Is there a way to go back and make points again from the layer with the selected buffers in QGIS?


texts:
I have .dxf file with contour lines representing elevation. But this contours are too "ugly" and elevation model looks just bad. As you can see on image, contours are not continuous...

Is there any way to clean up my contour lines so generated DEM will looks as it should?

So I tried terrain over TIN method and here is result. Unfortunately that pesky spots are still there. This is created from whole contour file, while first image is just cropped part of it


texts: Add an Excel file with a list of coordinates to QGIS as a layer?

texts:
I've added delimited text layers to QGIS before with columns consisting latitude and longitude separately. But Now I have an Excel file like:

station no parcel no    parcel coordinates
63         215          X1,Y1; X2,Y2; X3,Y3; X4,Y4

So I have all the parcel coordinates in one column for every parcel. I need to add this data with parcels to my previous shapefile with station no data. But I couldn't figure out how to import this coordinates to QGIS.


texts: Is my data locked if deployed elsewhere as a REST service?

texts:
I have a number of Image Services published in ArcGIS for Server 10.3.1.
They are built from mosaic datasets in file geodatabases on a NAS.

The DEM tiles used in these services are required to built another mosaic dataset, on a larger scale.
However, when performing "Add Rasters" to the mosaic dataset I receive a General Function Error.  As if the data I am trying to add is locked in the GDB.

Do I need to shut down the REST services connected to the same rasters I am trying to add to a new mosaic dataset?

Setting up in an SDE is not an option at this time.


texts: how to make a delay to wait for the map to load

texts:
What code could I use to make it wait for the program to load the map first after the zoomIn command?

Here's a part of my code.

for (var z = zoomfrom; z <= zoomto; z++) {
    var scale = MapPanel.map.getScale();
    var zoomlevel = MapPanel.map.getZoom();
    NewZoomLevel = getzoomlevel(scale);
    lat1 = lat2tile(coords.top, NewZoomLevel);
    lat2 = lat2tile(coords.bottom, NewZoomLevel);

    long1 = long2tile(coords.left, NewZoomLevel);
    long2 = long2tile(coords.right, NewZoomLevel);

    for (var x= lat1; x<=lat2; x++) {
         for (var y= long1; y<=long2; y++) {
              list = list + startlink + zoomlevel + '/' + x + '/' + y + '<br>';
        } //y for
    } //x for
    MapPanel.map.zoomIn();

    // this part im hoping to put a code to wait for map to load after zoomIn
}; // z for

I tried using delayed function of extjs but to no avail.


texts:
I have question about QGIS. I have georeferenced a fairly complex geological map (raster file) and would like to vectorize the different geological units. However, the georeferenced map contains text, contour lines, rivers, etc. in addition to the geological units that I would like to convert to shapefiles. I have tried to polygonize the map in bulk, but the resultant Shapefile is too bogged down by the extraneous data that I have no interest in to make any sense of.

Do you have any suggestions of how I can more easily do this without manually tracing all of the boundaries?

I have attached an image of the geological map below to give you all a sense of what I'm dealing with.


texts: Semantic for encoding geometry-event relation in a spatio-temporal dataset

texts:
I'm assembling a database of spatio-temporal events; think deforestation for instance.
I have different types of geometries that link to the events in different ways. I would like to encode that meta-information in some way, but I'm lacking the semantic for that.
Here are the geometry - event relations I have:

A polygon corresponds exactly to the spatial extent of the event
A polygon contains events but not the entire surface of the polygon is concerned by the events (e.g. A municipality that had some cases of deforestation)
A point sampled in the spatial extent of the event
A point loosely linked to the event (e.g. the centroid of a municipality in which there has been deforestation events)

Is there a metadata standard that include this kind of information?
If not, what terms would be suitable to encode these event-geometry relations?
What would be a suitable column name to store that type of meta information?