texts:
I'm trying to calculate relative relief from a raster in QGIS and I have troubles with square kilometer's cirle radius in r.neighbor.
In ArcMap it's done with Focal statistics - circular - map units (564) - range which is almost the same as r.neighbor.
Difference is I can't find a way to make 1 square kilometer radius in r-neighbor.
How do I make this happen?

So, problem with square kilometer is solved. Did anybody work with r.neighbor and knows how to apply equation for square kilometer in r.neighbor?


texts: Making ArcPy exclude frozen / off layers while importing DWG file?

texts:
I'm attempting to update existing GIS feature classes (SDE) using master AutoCAD files but end up importing ALL features present in the AutoCAD drawing regardless of the visibility settings (frozen / off) set in the DWG.

ArcMap recognizes the layer visibility settings if a DWG, how do I make arcpy do the same?

I'm using: arcpy.Append_management(,) to replace previously deleted features in the GIS.

outFC = "<outGDB>\\<outFC>"
inputDrawing = "<path>\\<file>.dwg"
inputLayer = os.path.basename(inputDrawing) + "_layer"
sqlExpression_CAD = "NOT \"Layer\" LIKE 'S_Grid%'" # This is what I want to avoid #

arcpy.MakeFeatureLayer_management(inputDrawing,inputLayer,<layer information>)
arcpy.SelectLayerByAttribute_management(inputLayer, "NEW_SELECTION", sqlExpression_CAD)
arcpy.Append_management(inputLayer, outFC, "NO_TEST")


texts: PBF vector tiles reach the client through Apache HTTP Mod_rewrite, but not shown in Leaflet. How to show these PBF in Leaflet?

texts:
In short:
This is about proxy-ing access to OpenMapTiles (OMT) by Apache Http. The proxy works hence the PBF tiles from OMT reach the client. However Leaflet fails to draw the PBF tiles in the map.

Below is the detailed:

I am responsible in providing OMT service in intranet of my organization. An existing & running internal web app (let's say TheApp) that requires login, currently add a new feature showing map with 2 elements:
- geoJSON queried from DB
- basemap from my OMT

TheApp is PHP running on Apache HTTP on Fedora.

Phase-1:

The javascript for showing the map specifies that client should get the OMT tiles directly from my OMT server. The result is just fine: the map shows geoJSON on top of OMT basemap as expected. Below is the diagram.

Phase-2:

The strategy is changed. The client has no direct access to OMT, and basemap tiles are provided through Apache Http Proxy. The module used in Apache HTTP is mod_proxy. The result is just fine: the map shows geoJSON on top of OMT basemap as expected. Below is the diagram.

Phace-3:

The strategy is changed again (but the diagram is similar to Phase-2). This is because in Phase-2, any browser can get the OMT tiles through the Reverse-proxy in the Apache HTTP and we don't like it. Now in Phase-3 we decide that only browser that logged-in in TheApp can get the OMT tiles.
My colleague (the admin of the Apache HTTP) finds that this requirement can not be achieved by mod_proxy, and hence uses mod_rewrite as suggested by several references. With mod_rewrite he configured that the required resources from OMT (including the PBF tiles) are rewritten to another address.

Now the result is: the mod_rewrite works, as the client re-requests to the specified address, and therefore the client got the PBF tiles.
However Leaflet fails to show the PBF tiles.

How to fix this? (Remember, the client got the PBF tiles). Below is the Dev Console of Chrome.

I suspect that Leaflet can not show the PBF tiles on map because under mod_rewrite the expected identifier of each PBF tile can not be understood by Leaflet. The Z/X/Y naming is still kept in the rewrite result, though.


texts:
The PAMAP dataset comes divided between two different projections: Pennsylvania North State Plane and Pennsylvania South State Plane.

I'd like a great big DEM of the whole state of Pennsylvania. gdal_merge.py can help me make this happen, but I'm going to need to reproject at least one half of the dataset first.

Is there a standard way of choosing how to do such a reprojection? Should I map everything to the North or South Plane or, perhaps, construct an intermediate projection and map everything to that?


texts:
I'm trying to create monthly composites of 2019 sentinel 2 imagery. For most areas this works fine; however, in areas where there are few images or they are cloudy and therefore masked I end up with an empty image for certain months. For those months I want to replace the missing monthly data with the annual greenest pixel image.
I then want to export these as 12 sets of point training data (one for each month) to use in a random forest. The code works; however, rather than give me the greenest pixel image for months with missing data, it returns a missing point. I'm sure in doing something wrong with the ee.Algorithms.If function
https://code.earthengine.google.com/c5c174e5e75fef4ad3269443574e2147


texts: Set of intersection points in trilateration / multilateration

texts:
Like in this post I will get a set of points (3 or more) with latitude, longitude and distance (in KM) information to find out the target location. A distance could be about 500 KM. The circles (unlike the above question) doesn't intersect in the same point which results in a set of points forming a region of a possible location.

Based on this explanation, I managed to write a java function that returns the set of intersection points for a set of circles (x,y,radius in CM) in the Cartesian coordinate system. But I don't know if I can adapt this function to get the intersection points of circles on the earth. Is this possible?


texts: Parse error while inserting values into polygon feature

texts:
I created a polygon as follows

create table boundary (id int, name text);

table created;

select AddGeometryColumn ('',boundary', 'latlon' 4326,'POLYGON', 2);

geometry column added;

then i'm trying to insert values as follows

INSERT INTO boundary (id , name, latlon) VALUES ('1', 'river', GeometryFromTerxt('POLYGON(82.240 16.988,82.243 16.991,82.240 16.988)', 4326));

finally i got an error as follow

  <--parse error at position (position number) within geometry.


texts: Automatic filling of rings with additional polygons in QGIS

texts:
I have a large layer with a couple of thousand polygons mapping inland waterbodies. A lot of these contain rings, mapping islands. I would like to count the number of islands and do some statistics with them. I vaguely remember that there was a tool that could automatically fill such rings with new polygons in QGIS 2.X, but it appears to be missing in QGIS 3.X. I did not find any plugin to do this job either. Is there any tool I have been missing, or can you point me to a different solution?

Edit: Please note that I want to do this automatically - I know, that I can fill rings by hand, but that won't do for thousands of them.


texts:
Where can I access location data for aquaculture sites across Europe?

Ideally there might be, though I haven't been able to see, a dataset held by the EU identifying every aquaculture site (sites such as shellfish farms) across Europe. The only place I have been able to find spatial information on this topic is at Marine Scotland National Marine Plan interactive map (for Scotland only).

The primary countries I need this information for are:
Norway, Sweden, Denmark, Germany, Netherlands, United Kingdom, Faroe Islands.

Currently I have looked through the GeoNames database, queried EEA and been googling away at the problem.

Has anyone seen datasets of this nature for the above mentioned countries?


texts: Calculate vertices' coordinates of a rectangle, given midline and broadth of rectangle

texts:
I have latitude and longitude of two points making the midline of a rectangle. Now I want to determine coordinates of 4 vertices of that rectangle. I know how to calculate on plane geometry by using line equation. But on geographic geometry, I'm new. I already read how to calculate one unit of latitude and longitude decimal degree but I'm still confused how to add them to a given point to calculate wanted points. My problem is like this:

Given:

Broadth of rectangle
Coordinates of two points making the midline, distance between these two points are equal to the width

Result: Coordinates of 4 vertices of rectangle

For example:
The given midline has two points:

A (10.767008,106.665884)
B (10.767715,106.667151)

Broadth is 4-metre long

All I could do now is to calculate the width of rectangle from those two points. How can I use width and height to calculate 4 vertices of this rectangle?

All of my calculated cases are on small scale (part of one street of a city) so I think altitude can be excluded


texts: Download/create address locator styles for Europe?

texts:
Where can I find some address locator styles for European countries? or how can I create my own address locator that contains street, house number and city, and use their aliases?

PS. I have the book: Customizing Locators in ArcGIS 10


texts: ArcGIS Geoprocessing Script - Calling External EXE Works in Desktop - Causes GP Service to Quit Abruptly

texts:
I have an interesting problem on my hands...  I have a generic python script that does a bunch of geoprocessing functions using the standard arcgisscripting module (Yes, I am still at 9.3.1).  At one point in this script, I call out to an external exe.  This external exe does some heavy duty raster processing using gdal17.dll.  After the .exe does its thing, my script continues on doing additional processing using the results from the .exe - still using arcgisscripting functions.

The tool runs beautifully in a desktop environment.  However, as soon as I publish it as an ArcGIS Server Geoprocessing Service - it exits abruptly in the first calling of arcgisscripting function AFTER the gdal-based exe is called.  I can't log or trap an error - it just plain 'ol quits.

I am very aware that mixing gdal in the arcgis geoprocessing environment is looking for trouble.  However, porting the logic from the exe is just not an option at this point - if it is even possible.  The fact the script works fine in the desktop environment gives me a glimmer of hope that it may work in a gp service.

Any suggestions greatly appreciated.


texts: Forcing boundary conditions in interpolated/kriged data?

texts:
I am using pykrige.ok import OrdinaryKriging to generate a spatially continuous field from sparse point data (see image below). The data (glacier ice thickness) are such that there is a hard boundary (red line, glacier margin) at which a specified value must be met (0 ice thickness). I delineated the glacier margin (a portion of these are connected by red lines in the image below) and wrote 0 ice thickness into the input dataset at each of those points.

I then run OK = OrdinaryKriging(eToget, nToget, hToget, variogram_model="linear", verbose=False, enable_plotting=False, ) (where 'Toget' suffixes indicate the data have been combined with the boundary values)  and z, ss = OK.execute("grid", vec_x, vec_y) to produce the output above.
As you can see, the kriging varies too smoothly, and produces non-zero values outside of the glacier margin, which is physically unacceptable.
Is there a way to force the kriging method (or some other spatial interpolation method) to match boundary values, such that everything outside the red lines would have a value of zero?
I know I could do this by masking/clipping the output, but there would be a non-physical discontinuity at the glacier margin if I did that, so I am looking for a way to produce smooth & continuous data up to the boundary, but are forced to match the boundary value.
I feel like this must be a relatively common issue (e.g., interpolating lake bathymetry, where 0 depth must be met along the lake boundaries), but I haven't been able to track down useful information in googling.


texts:
I want to create a light version of global map (only with cities and highways and natural objects). How to use osmosis correctly with multiple tags list? For example I try this command with no success:

osmosis \
--read-pbf file=1.pbf \
--tf accept-nodes place=* --tf accept-nodes highway=* \
--tf accept-ways place=* --tf accept-ways highway=* \
--tf reject-relations \
--write-xml output.osm


texts: Is it possible to geocode street intersections using Nominatim?

texts:
Is it possible to geocode the location of a street intersection using Nominatim?

Example geocode query

Commonwealth Ave at Washington Street, Boston, MA

Successful Geocode by Google

https://maps.google.com/maps?q=Commonwealth+Avenue+at+Washington+street,+Boston,+MA&hl=en&sll=42.351351,-71.118573&sspn=0.054868,0.132093&t=h&hnear=Commonwealth+Ave+%26+Washington+St,+Boston,+Suffolk,+Massachusetts+02135&z=17

Failed Geocode by Nominatim

http://nominatim.openstreetmap.org/search.php?q=Commonwealth+Ave+at+Washington+Street%2C+Boston%2C+MA&viewbox=-337.5%2C85.34%2C337.5%2C-80.5


texts:
I got a list of values and I could print them to the Console. It looks like this:

{
    "chlor_a": 0.12718113166405862
  },
  {
    "chlor_a": 0.37072379355468843
  },

How can I export it to my Drive account in csv format?

This is the code so far:

var bo=ee.FeatureCollection('users/broadsky2008/Border');

var chla = ee.ImageCollection("NASA/OCEANDATA/MODIS-Aqua/L3SMI")
.filterDate('2012-7-04', '2018-1-25')
.select('chlor_a');

var reduceR = function(image) {
  return image.set('number',image.reduceRegion({
  reducer: ee.Reducer.mean(),
  geometry: bo.geometry(),
 }));
};

// Map the function over the collection and display the result.
var mChla=chla.map(reduceR);

// Get a list of the numbers.
var numberList = mChla.aggregate_array('number');

print(numberList);


texts: Adding multiple sources to geometric network using ArcMap

texts:
I am using ArcMap.
I would like to identify all points that are on a river within 10km downstream of a starting point. I have multiple starting points and would like to do this from each of them.
I would like to add multiple sources, from one input file of around 2500 points, to a geometric network.
Do they need to be added as sources? If so, how?
Are there any other ways of doing it?


texts:
I realize there isn't much information here but I have an hdf5 file I want gdal to recognize in python.  I have tried the following:
import gdal

filename = 'I:/ppl/user/GFED/GFED4.1s_1997.hdf5'

os.path.exists(filename)

which returns True
and then I try to open:
dataset = gdal.Open(filename, gdal.GA_ReadOnly)

but this returns None.  I am at a loss on why it is not being recognized, and I am happy to share the file if it will be helpful.


texts:
I found a lot of solutions how to provide nice map tiles from local store without call external APIs (google, osm, and so on...).

Possibility which I found:

data from PostGIS to Geoserver will be styled with SLD and render as
tiles
store map tiles in image format in structured folders for different zoom levels.
store map tiles in mbtiles format and import it into Geoserver

Which one is the best solution? I try 2. and 3. option with Maperitive over command:

generate-tiles minzoom=10 maxzoom=17

into images and also into mbtiles, but label will be bad readable. Generated output will be in lower quality than maps which I can see in Maperitive.


texts: Does anyone know why I'm getting a completely black geotiff as a result of this shapefile cut out of another geotiff in rasterio?

texts:
I am using rasterio to cut out the shape of a lake outlined in my shapefile from a geotiff. When using a shapefile of a singular lake with a geotiff, this works perfectly, but not when cutting out a geotiff of a single lake with a different shapefile that contains a lot of different lake outlines. I'm wondering if there is a problem with the way my code is written or the files I'm using?

Here's my code:

import rasterio
import rasterio.mask
import fiona

with fiona.open("liag_shapes/liag_shapes.shp", "r") as shapefile:
    features = [feature["geometry"]for feature in shapefile]

with rasterio.open("lakesimcoe11.tif") as src:
    out_image, out_transform = rasterio.mask.mask(src, features,
filled=True, crop=True)
    out_meta = src.meta.copy()

out_meta.update({"driver": "GTiff",
             "height": out_image.shape[1],
             "width": out_image.shape[2],
             "transform": out_transform})

with rasterio.open("new_test.tif", "w", **out_meta) as dest:
    dest.write(out_image)


texts: How to clip a .osm.pbf file to a shapefile in R

texts:
I have (1) an  OpenStreetMap file of a large area in .osm.pbf format that I would like to clip with (2) a shapefile .shp of a smaller area using R.

The output output I would like to get is a .osm.pbf file of the intersection between (1) and (2). Any ideas on how to do this in R ?

Reproducible example:

Here is the .osm.pbf file of the larger area
Here is the shapefile file of the larger area


texts: Remove Overlapping Polygons in Google Earth Engine

texts:
I have two feature collections for protected areas in the Amazon Rainforest. One of them is filtered so that only the protected areas within a drawn polygon show, the other is filtered for the entire Amazon Rainforest. Obviously, the one for the entire Amazon rainforset overlaps with the one inside the drawn polygon.
Is there any way to essentially take one group of polygons away from another so that the two are separated?
Here is a link to my code. I want to try and take variable 'AAAProtectedAreas' away from 'AmazonProtectedAreas'


texts: Why does the Extent of NIR Band 8 Landsat image not match other bands, even when cropping to a specified extent?

texts:
I'm trying to load and then stack four bands 2,3,4 and 8 RGB and NIR from Landsat images. I have downloaded and extracted them and they are in .tif format.

# loading the bands
B2 <- raster("./Data/LC08_L1TP_015053_20160413_20170223_01_T1_B2.tif")
B3 <- raster("./Data/LC08_L1TP_015053_20160413_20170223_01_T1_B3.tif")
B4 <- raster("./Data/LC08_L1TP_015053_20160413_20170223_01_T1_B4.tif")
B8 <- raster("./Data/LC08_L1TP_015053_20160413_20170223_01_T1_B8.tif")

# then I want to clip the rasters to the study site
# first I have to load the raster of the Costa Rica provinces
shape <- readOGR("./Data/gadm36_CRI_1.shp")

# Then add that projection to the shapefile
shape.proj <- spTransform(shape, CRS(proj4string(B2)))
proj4string(shape.proj)

# Combining the two provinces that I want : Heredia and Alajuela
p1 <- shape.proj[shape.proj$NAME_1=='Heredia', ]
p2 <- shape.proj[shape.proj$NAME_1=='Alajuela', ]
p2$NAME_1 <- NULL
x <- bind(p1, p2)

# Crop raster by extent of shapefile
studysiteB2 <- mask(crop(B2, x), x)
studysiteB3 <- mask(crop(B3, x), x)
studysiteB4 <- mask(crop(B4, x), x)
studysiteB8 <- mask(crop(B8, x), x)

# now stack the Bands together
image <- stack(studysiteB2,studysiteB3,studysiteB4,studysiteB8)

But then I am getting this error message

Error in compareRaster(x) : different extent

When I look at the extent of the cropped bands 2,3 and 4 all have the same extent: e.g.

extent(studysiteB2)

class       : Extent
xmin        : 710385
xmax        : 860325
ymin        : 1087635
ymax        : 1225725

But the extent of B8 cropped image is .5 more on the xmin / xmax

extent(studysiteB8)

class       : Extent
xmin        : 710392.5
xmax        : 860332.5
ymin        : 1087628
ymax        : 1225733

This is the same on the original images as well. Any ideas on why this is happening and how to solve it?

Here is a DropBox link to the images: https://www.dropbox.com/sh/fk2acf07lc9437l/AAD-DCopQxFFIwHRWAfew2PMa?dl=0


texts: Error parameter invalid in downloading using google earth engine

texts:
I am trying to download a soil moisture layer for a single day from google earth engine using the following code:
var dataset = ee.ImageCollection('NASA_USDA/HSL/SMAP_soil_moisture')
                  .filter(ee.Filter.date('2010-02-01'));
var soilMoisture = dataset.select('ssm');

// Create a geometry representing an export region.
var geometry = ee.Geometry.Rectangle([-76.5060023341, 3.67837906194, -76.0152561351, 3.16455715204]);

// Export the image, specifying scale and region.
Export.image.toDrive({
  image: dataset.select('ssm').first(),
  description: 'soilMoisture',
  scale: 250,
  folder: 'google_earth',
  region: geometry,
  fileFormat: 'GeoTIFF'
})

but I got the following error: Error:  Image.clipToBoundsAndScale: Parameter 'input' is required.
I have checked the records of the NASA_USDA/HSL/SMAP_soil_moisture and it goes from Jan 1, 2010 - Nov 16, 2020, therefore it should be a tiff for that day.
Any idea why it is not working?


texts: Error when using variable in SQL clause with arcpy.MakeQueryTable_management

texts:
When I hard code a date my SQL query works fine

where_week = "dtrepor > ‘2017-04-21 07:00:00’"

but when I try to use a variable it booms out.
This variable prints to the exact same date and time but doesn’t work when substituted

my_OneWeekAgo = my_now - datetime.timedelta(days=7)
SQLdateFrom = my_OneWeekAgo.strftime("%Y-%m-%d 07:00:00")
print SQLdateFrom

2017-04-21 07:00:00

This errors out

where_week = "dtrepor > 'SQLdateFrom'"

Here is the error

An underlying database error occurred. [C:\Users\mandrews\Documents\WeeklyIncidentMaps\WeeklyVehicleTheft.gdb\MyRecords]
Failed to execute (MakeQueryTable).


texts: Extending polylines with points from Excel/Access file

texts:
I have a .shp-file filled with polylines, where some of them are currently being drawn as a straight line from A-B, where in reality they curve a lot along the way.

In the shapefile, I have an "ID" that matches with an Access/Excel table where I have the coordinates for all the vertices that need to be added to the lines in order for them to have their correct position.

I need to change the lines from straight lines to lines with a lot of vertices using coordinates from my Acces/Excel-file. CRS is the same for both files, and there are about 15.000 vertices in total.

I am using ArcMap 10.6.1, Advanced license, with Spatial/3D Analyst.


texts: Installing ArcGIS for Desktop with existing ActivePython install?

texts:
I just did a clean OS install and installed ActivePython 2.7.2

I want to use this Python install for use with ArcGIS 10.1

I just installed ArcGIS 10.1 and it did a new Python install.

What do I need to do to get Python for ArcGIS to install on my current Python?


texts: I cannot find this information on the internet. Is there a way to migrate project from QGIS to ArcMap? At least layer symbology. Shapefiles not the issue. On the internet I find the contrariwise info, where usually people are looking a way from ArcGIS to QGIS.


texts: Create a flow network from a flow direction ARC GIS

texts:
First, I created a flow_direction from a DEM, then I applied corrections to this flow direction.
Now my question, how to create a flow network (order flow) from a flow direction, that I already corrected?
What tool to use on ARC GIS: ArcToolbox > Spatial Analyst Tools> Hydrology> ....


texts: How to write PostreSQL source in virtual layer in QGIS

texts: In QGIS 2.14, I am trying to create a connection with my PostreSQL database using "Virtual Layer". But in the "Virtual Layer" window, how should I write the parameters in source field ?


texts: QGIS install on OSX with homebrew: "You have other linked versions!"

texts:
I have an OSX Mojave setup where I've already installed PostGIS and GDAL using Homebrew.

Now I'd like to install QGIS. I'm trying:

brew install qgis

But... I'm getting lots of warnings of the form:

osgeo-postgis: You have other linked versions!

I'm using PostGIS in various projects, so I would prefer to avoid breaking / version-changing anything that's currently running if possible, though I can do so if absolutely necessary.

What do I need to do to install QGIS?

Here's the full output:

==> Installing osgeo-qgis from osgeo/osgeo4mac
osgeo-proj: You have other linked versions!

Unlink with brew unlink proj or remove with brew uninstall --ignore-dependencies proj

osgeo-libgeotiff: You have other linked versions!

Unlink with brew unlink libgeotiff or remove with brew uninstall --ignore-dependencies libgeotiff

osgeo-libspatialite: You have other linked versions!

Unlink with brew unlink libspatialite or remove with brew uninstall --ignore-dependencies libspatialite

ant: Java 1.8+ is required to install this formula.
JavaRequirement unsatisfied!
You can install with Homebrew Cask:
  brew cask install java
You can download from:
  https://www.oracle.com/technetwork/java/javase/downloads/index.html
osgeo-netcdf: You have other linked versions!

Unlink with brew unlink netcdf or remove with brew uninstall --ignore-dependencies netcdf

osgeo-postgresql: You have other linked versions!

Unlink with brew unlink postgresql or remove with brew uninstall --ignore-dependencies postgresql

osgeo-gdal: You have other linked versions!

Unlink with brew unlink gdal or remove with brew uninstall --ignore-dependencies gdal

osgeo-postgis: You have other linked versions!

Unlink with brew unlink postgis or remove with brew uninstall --ignore-dependencies postgis

osgeo-grass: XQuartz 2.7.11 (or newer) is required to install this formula. X11Requirement unsatisfied!
You can install with Homebrew Cask:
  brew cask install xquartz
You can download from:
  https://xquartz.macosforge.org
mesa: XQuartz 2.7.11 (or newer) is required to install this formula. X11Requirement unsatisfied!
You can install with Homebrew Cask:
  brew cask install xquartz
You can download from:
  https://xquartz.macosforge.org
osgeo-openscenegraph: XQuartz 2.7.11 (or newer) is required to install this formula. X11Requirement unsatisfied!
You can install with Homebrew Cask:
  brew cask install xquartz
You can download from:
  https://xquartz.macosforge.org
freeglut: XQuartz 2.7.11 (or newer) is required to install this formula. X11Requirement unsatisfied!
You can install with Homebrew Cask:
  brew cask install xquartz
You can download from:
  https://xquartz.macosforge.org
osgeo-proj@5: You have other linked versions!

Unlink with brew unlink proj or remove with brew uninstall --ignore-dependencies proj


texts: OpenLayers. Slow map after render heatmap even when I clear the layer source

texts:
I'm working in a map with about 10 layers and it works fine, when I create a heatmap it becomes slower, that's something I was expecting, because the source has about 2000 points. When I want to make the heatmap disapear I just clear the layer source, but the map continues very slow.

Check out my code:

let map = new ol.Map({
    controls: ol.control.defaults({attribution:false}).extend([
        expandedAttribution
    ]),
    target: document.getElementById('map'),
    renderer: 'canvas',
    overlays: [overlayPopup],
    layers: layersList,
    view: new ol.View({
     extent: extent,
     maxZoom: maxZoom,
     minZoom: minZoom,
    }),
    interactions: ol.interaction.defaults({
        doubleClickZoom: false
    })
});

let data = new ol.source.Vector();
//I create the heatmap with a empty vector
let heatMapLayer = new ol.layer.Heatmap({
    source: data,
    radius: 10
});
map.addLayer(heatMapLayer);

//*******************************************************************
//To add the source
let coord = ol.proj.fromLonLat([ocor.longitude, ocor.latitude]);
let lonLat = new ol.geom.Point(coord);
let pointFeature = new ol.Feature({
    geometry: lonLat,
    weight: 20 // e.g. temperature
});
dataHeatmap.addFeature(pointFeature);
heatMapLayer.setSource(dataHeatmap);

//******************************************************************
//and to remove the heatmap:

heatMapLayer.getSource().clear();


texts: Using multiValue in Python Toolbox with parameters interaction

texts:
I want to select out the records from A.shp if the where clause is satisfied. First, get the unique values (e.g., [2011,2012,2013,2014]) from a field (e.g., Year) and the field is from the A.shp, and then give the user to select the values as he wants. Suppose he is interested in the dataset from the year 2011 and 2012. Based on the values selected, delete other records (Year == 2013 or 2014 in this case) in the A.shp.

The code below has several problems.
1. Select All and Unselect All option in the tool window don't work.
2. The select out values don't change if the user unchecked several of them.
3. The UpdateCursor part removes all records rather than reflecting the user interests.

import arcpy

class Toolbox(object):
    def __init__(self):
        """Define the toolbox (the name of the toolbox is the name of the
        .pyt file)."""
        self.label = "Toolbox"
        self.alias = ""

        # List of tool classes associated with this toolbox
        self.tools = [Tool]

class Tool(object):
    def __init__(self):
        """Define the tool (tool name is the name of the class)."""
        self.label = "Tool"
        self.description = ""
        self.canRunInBackground = False

    def getParameterInfo(self):
        """Define parameter definitions"""
        #fc
        fc = arcpy.Parameter(
            displayName="FC:",
            name="fc",
            datatype="Feature Class",
            parameterType="Required",
            direction="Input")

        # Plat category field parameter
        fld = arcpy.Parameter(
            displayName="Field from the Plat Features",
            name="fld",
            datatype="Field",
            parameterType="Required",
            direction="Input")
        fld.parameterDependencies = [fc.name]

        # MultiValue
        values = arcpy.Parameter(
            displayName="Values",
            name="values",
            datatype="String",
            parameterType="Required",
            direction="Input",
            multiValue=1)

        params = [fc, fld, values]
        return params

    def isLicensed(self):
        """Set whether tool is licensed to execute."""
        return True

    def updateParameters(self, parameters):
        """Modify the values and properties of parameters before internal
        validation is performed.  This method is called whenever a parameter
        has been changed."""

        if parameters[0].value and parameters[1].value:
            fc = str(parameters[0].value)
            col = str(parameters[1].value)
            vals = sorted(set(row[0] for row in arcpy.da.SearchCursor(fc,[col]) if row))
            parameters[2].filter.list = vals
            parameters[2].value = vals
        return

    def updateMessages(self, parameters):
        """Modify the messages created by internal validation for each tool
        parameter.  This method is called after internal validation."""
        return

    def execute(self, parameters, messages):
        """The source code of the tool."""
        values = parameters[2].valueAsText
        list = values.split(';')

        for x in list:
            arcpy.AddMessage(x)

        fc = parameters[0].valueAsText
        field = parameters[1].valueAsText
        with arcpy.da.UpdateCursor(fc,[field]) as cursor:
            for row in cursor:
                if not row[0] in list:
                    cursor.deleteRow()

        return

Any suggestion?


texts:
I am working on QGIS - 1.8 Version and frequently getting an error while using Road Graph and Shortest Path Opetions.

Kindly Suggest or provide solution..

With Regards

Prashant Persai


texts: Access secure Geoserver WMS from ArcGIS using LDAP Domain user

texts:
When I use public workspace:

http://(myserver):8080/geoserver/public/wms?

The layers are displayed in ArcMap, but when I switch to a private workspace (imagery_mp) using:

http://(myserver):8080/geoserver/imagery_mp/wms?
I put my credentials and I get: No resource could be found at that address. (status code 404).

http://(domain/user:password@)(myserver):8080/geoserver/imagery_mp/wms?
I get: Internal server error. (status code 500).

Then I try:
http://(domain\user:password@)(myserver):8080/geoserver/imagery_mp/wms?
I get: No resource could be found at that address. (status code 404).

Then I try:
http://(myuser@mydomain.com:password@)(myserver):8080/geoserver/imagery_mp/wms?
I get: Unidentified error

QGIS displays correctly using boxes for credentials, my setup is ArcMap 10.4 and GEOSERVER 2.11.1.


texts: Splitting polylines according to coincident (or touching) points with other polylines?

texts:
The question is unfortunately not answered yet.

I'm looking for a method or python script to split polylines according to coincident points with other polylines, but avoid splitting the ones which just intersect (without a sharing point like the planarize or intersection tool in ArcGIS would produce). The "touches (second_geometry)" function might be a way but I have at the moment no clue how to solve it.

The image illustrates my aim. The yellow line hasn't changed, because it has no coincicent points with the red line. While the red line as well as the dark blue line is split into 4 different polylines, because both share a coincident point.

Thanks for any answer in advance!


texts: Modifying a WMS feature in POSTGIS with Openlayers

texts:
I have a WMS layer that i'm using to display complex features to the users. It was really the main reason i chose GeoServer, because rendering all the vectors, made it painfully slow. Like GeoServer converts all features to raster tiles things got much faster.

This works fine to visualize the information and i even draw a vector over the selected shape, using the features geometry attributes to highlight it like so:

But I want also to have an edit mode to change the attributes. Like this:

I searched other discussions and everybody is talking about WFS-T to change vectors, and it seems to me a bit strange to limit the editing of features only to vectors.

If one can identify the feature even if it's raster using WMSGetFeatureInfo to create the popup is it so difficult to change its attributes?

I have the id of the selected feature so can one connect and update directly the POSTGIS table? Is this possible?

Update table where id = etc etc ?

What kind of approach do you recommend?


texts: QField: Draw line/polygon freehand using a digital pen on Android Tablet

texts: I have a Samsung Tab Active (Android) with pen and would like to draw lines or polygons freehand. Not just tabbing and adding a new node. Is that possible or planned?


texts: Locks not being released in Linux held File Geodatabases

texts:
We have a system with 2 SOC servers, and one SOM server, acting on a shared storage area.
We have had performance issues for years, and believe we have worked out why.
There are several file Geodatabases on the Linux box, and on these fgdb, were massive amounts of unreleased lock files (up to 250,000 lock files per gdb).
We're currently running diagnostics on the servers.
I was sure ArcGIS would be releasing the lock files when they were no longer needed.
Even hard booting the system led to these lock files still being kept in Linux.
I am no linux expert, and have put this to the linux team, but does anyone know of this, have experienced this or seen this before?


texts: Running Python scripts (with parameters) within another Python script with ArcPy?

texts:
A common coding pattern used in AML was to run an AML (with parameters) inside another AML.

An application that I'm currently developing would benefit from being able to run a Python script (with parameters) within another Python script.

However, this does not seem to be at all straightforward.

Using ArcGIS 10, I'm experimenting with wrapping the "inner" Python script into an ArcGIS tool that has the parameters.  I thought it would be a simple matter to have the "outer" Python script use arcpy.ImportToolbox to import the toolbox and then run the tool(s) within it.  However, in testing so far all my attempts to run the "inner" tool from the "outer" script appear to simply skip the "inner" tool (no error is thrown).

Here is some test code to try and illustrate better what I am trying to describe.

My testinner.py script is:

inputString = arcpy.GetParameterAsText(0)

newFC = "C:\\Temp\\test.gdb\\" + inputString
arcpy.Copy_management("C:\\Temp\\test.gdb\\test",newFC)

My testouter.py script is:

import arcpy

inputString1 = arcpy.GetParameterAsText(0)
inputString2 = arcpy.GetParameterAsText(1)

arcpy.ImportToolbox("C:\\Temp\\test.tbx")

arcpy.testinner_test(inputString1)

arcpy.testinner_test(inputString2)

For testinner.py its tool needs a single String parameter.

For testouter.py its tool needs two String parameters

The two tools are placed in a test.tbx.

The test.gdb just needs a single empty feature class called test.

Once you have the above assembled, running the testinner tool with a string like 'abc' passed in as its parameter should result in feature class 'test' being copied to one called 'abc' OK.

But when you try running the testouter tool with two strings like 'uvw' and 'xyz' as its parameters, the testinner tool within testouter.py seems to run OK once, but sends ArcMap 10 SP2 on Vista SP2 to a Serious Application Error when trying to use it the second time.

The same test using Windows XP SP3 and ArcGIS Desktop 10 SP2 also produces a Serious Application Error at the same point.


texts: Received a python error when I clicked on new catalog in QGIS Geoserver

texts:
What shall I do if I Python error message is there when I clicked on new catalog in QGIS Geoserver.

I have the same message: Python error: An error has occurred while executing Python code: See message log (Python Error) for more details.


texts: Showing polyline direction in ArcGIS for Desktop (ArcMap)?

texts:
How can I show the direction of a polyline as it was digitised?

I need this to check if my network analysis is obeying one way streets.

I know of a tool in Xtools that can do this but my trial period just finshed.


texts: Setting color assignment order with ArcPy?

texts:
I have a shapefile feature class with the following attribute table:

I set the shapefile symbology using these lines:

lyr.symbology.valueField = "GRIDCODE"
lyr.symbology.addAllValues()

What i get is this:

The color range is not in ascending order. It is starting at "10".

When I click on the color range ("Farbverlauf"), it changes to an ascending order, starting at "2". However, I have too many shapefiles to do this for each single one.

How can I adjust this color assignment in ArcPy?

EDIT

I'm using this workaround now:
How To: Add leading zeros to values in a field in an attribute table

Example code: fill each row of field "GRIDCODE" of layer "lyr" with "002"

with arcpy.da.UpdateCursor(lyr, "GRIDCODE") as cursor:
    for row in cursor:
        row[0] = (str(2).zfill(3))
         cursor.updateRow(row)


texts: SIP failed to load Ubuntu 20.04.1 QGIS 3.16.1

texts:
Recently I updated my Ubuntu to 20.04.1 and QGIS to 3.16.1. If I start QGIS now I get this message:
Das SIP-Modul konnte nicht geladen werden.
Die Python-Unterstützung wird abgeschaltet.

Traceback (most recent call last):
  File "", line 1, in
  File "/usr/lib/python3/dist-packages/qgis/__init__.py", line 70, in
    from qgis.PyQt import QtCore
  File "/usr/lib/python3/dist-packages/qgis/PyQt/QtCore.py", line 24, in
    from PyQt5.QtCore import *
ImportError: /home/kai/.local/lib/python3.6/site-packages/PyQt5/QtCore.so: undefined symbol: _ZdlPvm, version Qt_5

Python-Version:
3.6.9 (default, Oct  8 2020, 12:12:24)
[GCC 8.4.0]

QGIS-Version:
3.16.1-Hannover 'Hannover', b381a90dca

Python-Pfad:
['/usr/share/qgis/python', '/home/kai/.local/share/QGIS/QGIS3/profiles/default/python', '/home/kai/.local/share/QGIS/QGIS3/profiles/default/python/plugins', '/usr/share/qgis/python/plugins', '/usr/lib/python36.zip', '/usr/lib/python3.6', '/usr/lib/python3.6/lib-dynload', '/home/kai/.local/lib/python3.6/site-packages', '/usr/local/lib/python3.6/dist-packages', '/usr/lib/python3/dist-packages']

The question I found is over 8 years old, is there now a day an easy solution?

If i execute this command:
ls /usr/bin | grep python

i got these lines back
dh_python2
python
python2
python2.7
python2.7-config
python2-config
python3
python3.8
python3.8-config
python3-config
python3-futurize
python3-pasteurize
python-config
x86_64-linux-gnu-python2.7-config
x86_64-linux-gnu-python2-config
x86_64-linux-gnu-python3.8-config
x86_64-linux-gnu-python3-config

At the error message python version 3.6.9 is searched.
Is my python to new, for qgis ?


texts: Drawing Queried Geometry using ArcGIS API for JavaScript?

texts:
I am new to all of coding and I have a rest point with all of my park polygons, what I want to do is be able to  select a park  name from a drop down menu and query out and display the single park that is chosen. I believe I have everything set up to query my geometry, now I don't know how to get it to draw the queried feature. right now I have 2 different drop downs. According to my professor i have the selection stored in "selectedParkFeature" variable, and all I have to do now is to draw it on the map. I have researched and I am stuck. I am doing all of this on ArcGIS API for JavaScript 3.17.

    query("#selectParkName").on('change', function (evt) {
    var tempParkName = evt.target.value;
    getParkFeature(tempParkName);
});

//function to select the park
function getParkFeature(in_park_name) {
    var query = new Query();
    query.where = "NAME='" + in_park_name + "'";
    query.outFields = ["*"];
    park.queryFeatures(query, function (featureSet) {
        console.log(featureSet);
        selectedParkFeature = featureSet.features[0];
    });
}


texts:
I have a netcdf4 data set that contains array and information about latitude and longitude which can be accessed with the following R command:

ncvar_get(dataname, 'lon');
ncvar_get(dataname, 'lat');

where

dataname <- nc_open(paste(filename,".nc",sep=""))

I have another file containing the location (with latitude and longitude also) of the grid contained in dataname.
I'd like to implement interpolation (spline, kriging, etc.) using these two files but can't find any tutorial which is clear about the subject.

Could anyone guide me?

(package that should be used, I have seen gstat and gam but examples provided are far away similar to what I want to do, practical implementation)

Thanks in advance.


texts: Image not loading properly with OpenLayers Map in Chrome

texts:
I have an application using OpenLayers 2.13 that load serveral maps. The first image loads and displays fine. The second dosen't load properly. See image below:

I am new to OpenLayers. The application uses the same code to load and build all the maps but only the first one seem to load properly. It seems to work just fine in IE and FireFox but in Chrome it loads and displays like the above image. Any help would be greatly appreciated.


texts: Difference between quotation marks ('single' vs "double") in QGIS

texts:
I am not a new QGIS user, however, I am still confused concerning the quotation marks syntax ('single' vs "double").

I know something about the quotation marks and the way they are applied in Python scripting (i.e. comments, strings etc.).

There are even some related questions:

Single quotes vs. double quotes in Python [closed]
Is there any difference between “string” and 'string' in Python? [duplicate]

Nevertheless, I would like to clarify the details regarding the quotation marks ('single' vs "double") in the scope of QGIS.

For instance, if I type something in the Expression dialogue I will have a different output.

"something", it is written as

The Output is NULL.

'something', it is shown as

The output is 'something' as a string.

So, my questions are:

What is the syntax/meaningful difference between quotation marks ('single' vs "double") in QGIS?
Are there any other circumstances when other types of quotation marks are used in QGIS, for instance, triple quotation marks as in Python?


texts: QGIS + GRASS: r.in.gdal "View output" produces "cellhd file (...) does not exist"

texts:
I'm working with a recently installed QGIS (2.10.1) + GRASS (6 v2.0) package, and have been working with some rasters and vectors in QGIS that I'd now like to play with using the GRASS tools.

I've set up my GRASS location, mapset etc, and made sure the region and projection are all correct - no issues there. I've then used r.in.gdal.qgis to load a raster from my QGIS environment, and that's proceeded without giving me any errors. However, nothing appears on the screen after that command is completed, and when I click "View output" I get an "Invalid Layer" error. The message I'm getting (with some URI parts removed) is as follows:

GRASS provider: cellhd file D:...\GRASS/Hopefield_GRASS/Hopefield_GRASS/cellhd/aster_GRASS does not exist

Raster layer: Provider is not valid (provider: grassraster, URI: D:...\GRASS/Hopefield_GRASS/Hopefield_GRASS/cellhd/aster_GRASS

I'm thinking maybe there's an issue with the path having both forward and backward slashes in it - perhaps that's getting GRASS confused?


texts: setFilterExpression not giving results because of space?

texts:
I am trying to select data based on QgsVectorLayer() getting it's features from layer.getFeatures().next() method.

When I pass expression in that using QgsFeatureRequest().setFilterExpression(expression) where NAME column contains a space, like

NAME = 'ABC XYZ'

the filter doesn't return any rows, even when such value is present.
When I pass

NAME = 'DEF'

It works fine.

Can anyone help?


texts: How to highlight marker dynamically with leaflet and leaflet draw

texts:
I am trying to make a selection via drawing a rectangle (or any shape) with leaflet draw. The selected features should be highlighted. I saw that there are various functions within the leaflet source that could help (contains from Bounds class and _containsPoint from Poly) however I keep getting an error saying that layer does not have this function. I have tried putting the draw layer into a new variable by creating a new L.Polygon object and calling _containsPoint(), but this gives a similar error - saying that the function does not exist for that polygon.

var marker1 = L.marker([51.441767, 5.470247],{
    draggable: true
}).addTo(map);

var someLayerGroup = L.featureGroup();
marker1.addTo(someLayerGroup);

map.addControl(new L.Control.Draw({
    edit: {
    featureGroup: drawnItems,
        edit: false
    },
    draw: {
        marker: false,
        circlemarker: false,
        polyline: false,
        circle:  false,
        polygon: {
          allowIntersection: false,
          showArea: true
        }
    }
}));

map.on(L.Draw.Event.CREATED, function (event) {
    var layer2 = event.layer;
    var type = event.layerType;
    if(type==='rectangle'){
        //do rectangle specific actions
        someLayerGroup.eachLayer(function(layer){
            if(layer2.contains(layer)){
            //set style of the point
            layer.setStyle({fillColor: '#CCCC00'});
            }
        });
    }
    drawnItems.addLayer(layer2);
});

JsFiddle:
http://jsfiddle.net/ve2huzxw/2109/


texts: static map api "Invalid or nonexistent map configuration token "

texts:
Im trying to create static map and im having problem:

  Invalid or nonexistent map configuration token

When I publish map I get two links and they work fine with token in link https://rokiszb.carto.com/builder/ec62c84b-c6ec-41b3-9468-353b608c88fe/embed . But this token doesnt work for my static map link https://rokiszb.carto.com/api/v1/map/static/center/ec62c84b-c6ec-41b3-9468-353b608c88fe/10/54.1/24.1/1500/1500.jpg


texts: Calculating monthly mean temperature from TerraClimate

texts: I'm working on a hydrology project and need to obtain monthly average temperature data for a watershed. I'm using the TerraClimate dataset, but it only contains monthly average high and low temperatures. Where can I download monthly mean temperature data, or how can I calculate it?


texts:
I am trying to rotate the Raster Layer loaded by GDAL plugin.
var layer = MainMap.Layers.SelectedLayer as IMapRasterLayer;
layer.DataSet.RotateAt(angle, CenterCoordinate);

This did not rotate keeping the center as it is given in the parameter of RotateAt.
I need something like below


texts: Postgis query not using the spatial index

texts:
I've been using postgres8.3/postgis1.5 for some time, and it works well. I've moved my entire database to a new, faster machine and upgraded to postgres9.3/postgis2.1 using pg_dump/postgis_restore.pl, which also worked well. Everything else seems to be running faster, but a particular query is now taking extremely long to complete - so naturally I checked that it was correctly indexed and the tables were vacuum+analyzed. VACUUM ANALYZE public.places (the_geom) doesn't change the result. Using EXPLAIN, I see that the 9.3/2.1 instance isn't even using a spatial index.

The public.places table is identical in both databases; same structure and exactly the same records (236,884 rows).

I don't doubt that I can write a more efficient query, but before I launch into that I'd like to at least get it running as fast as it was in 8.3/1.5 (700ms vs 580000ms). Suggestions? Does anyone see my error? What has changed in postgis2.1 that keeps the spatial index from being utilized in this query?

The sql:

SELECT max(foo.pname) as pname, split_part(attribute, '|', 2) as zipcode FROM places,
(SELECT features.pname as pname, ST_Buffer(the_geom, -.001) as the_geom FROM places p,
    (SELECT pname, pid FROM search WHERE gid = '50873' LIMIT 1) as features
WHERE p.pid = ANY(features.pid)) as foo
WHERE places.ptype = 'Zip Code Areas'
AND ST_Intersects(places.the_geom, foo.the_geom)
GROUP BY zipcode
ORDER BY zipcode;

The table schema:

CREATE TABLE public.places
(
  pid integer NOT NULL DEFAULT nextval('places_pid_seq'::regclass),
  ptype text,
  pname text,
  the_geom geometry,
  sid integer NOT NULL,
  pcountry text,
  woeid integer,
  attribute text,
  CONSTRAINT places_pkey PRIMARY KEY (pid),
  CONSTRAINT enforce_dims_the_geom CHECK (st_ndims(the_geom) = 2),
  CONSTRAINT enforce_srid_the_geom CHECK (st_srid(the_geom) = 4326)
)
WITH (
  OIDS=FALSE
);
ALTER TABLE public.places OWNER TO xxxx;
GRANT ALL ON TABLE public.places TO xxxx;
GRANT SELECT ON TABLE public.places TO xxx;

CREATE INDEX places_pcountry_idx ON public.places
  USING btree (pcountry COLLATE pg_catalog."default");

CREATE INDEX places_ptype_idx ON public.places
  USING btree (pcountry COLLATE pg_catalog."default");

CREATE INDEX places_ptype_key ON public.places
  USING btree (ptype COLLATE pg_catalog."default");

CREATE INDEX places_spatial_index ON public.places
  USING gist (the_geom);

CREATE INDEX places_woeid_idx ON public.places
  USING btree (woeid);

EXPLAIN 8.3/1.5:

EXPLAIN 9.3/2.1:


texts: Where did the polygons go after projecting a map in QGIS?

texts:
I'm fairly new at GIS know-how, but I've run into a problem with a projection of a world map.  I have a shapefile of all countries I downloaded from Natural Earth.  After projecting it into a orthographic as suggested in this other question, I changed the lat/long so that the Atlantic would be somewhat in the center.

However, after changing it, the United States, Mexico and Canada disappear from the map, as you can see here:

When I reproject the map to a long over the US, the polygons appear again.  Any suggestions?


texts: Random sample of points to calculate K function on a large point set

texts:
For a point pattern analysis in R, I want to study a relatively large point dataset (about 1 million points) with a cross K function to compare different categories of points. Including a Montecarlo simulation for envelope, the calculation of K on this dataset is very intensive, and I would like to sample the data to be able to calculate it faster.
How can I estimate the error in the result of K if I take a random sample of points (e.g., 1% or 5% of the data)? Can it be derived mathematically or should I run a Montecarlo simulation to see how much the results differ?


texts: Using Google Maps to geocode for ArcGIS?

texts:
I've been trying to use Google Maps for geocoding using the steps described in this tutorial:
http://sites.tufts.edu/gis/files/2013/11/Geocoding-addresses-with-Google-Maps1.pdf

I'm a little confused because what I'm trying isn't working, and all the additional documentation I have found online about using Google Maps for geocoding seems to require scripting.  Is the above tutorial out of date?  If not, what am I missing.

Here's what I've tried.

First I create a .csv file with addresses like this:

Name,Location,Address,City,State,Zip_Code
Boston City Hall,Massachusetts,1 City Hall Square,Boston,MA,02201
Massachusetts State House,Massachusetts,24 Beacon St.,Boston,MA,02133

Then I open up MyPlaces in my Google account and import the file.  This works as expected.

Then I export the layer as a .kml.

In ArcMap, I try to use KML to Layer to convert the .kml, but the conversion fails with the error:

ERROR 000401: No features were found in processing

I've opened up the .kml files with a text editor, and I notice that there don't seem to be any coordinates present, which I assume is related, but I don't actually know.  What is going on here, and what is the solution?


texts: Determining orientation angle of each street on dataset in ArcGIS Desktop?

texts:
I am trying to join centerline attributes to a blockface shapefile using the "near" tool. The problem is that on many streets, the tool gets confused because of intersections and varying width. My idea is to select streets that are oriented similarly and run the near tool, omitting any street segment that might confuse the tool. I need to affix the angle of orientation to each individual segment in both the centerline and the blockface layers.

I'm not very good at python and am hoping to find an answer within ArcGIS 10.5


texts:
I am using R for doing a map of groundwater levels which means I have some measurements at some coordinates in the Gauss-Krueger system and I interpolated between these measurements.
I want to plot this interpolated distribution map on top of a satellite image of that area. Using ggmap, I can only retrieve the map for the corresponding Lon-Lat coordinates, but I actually would like to plot the whole thing in the GK system. Is it possible to do this, or alternatively, convert the axes of the satellite map somehow? I thought about converting the GK coordinates to Lon-Lat, but I am not sure about the accuracy of that, since the area of interest is pretty small, maybe 300x300 meters.


texts: Colouring points in the same colour as their nearest point/polygon from another layer in QGIS

texts:
I got two files, a multi-polygon file with polygons and points (File 1) and a file with points only (File 2).
Our task is to fill the color of the points from File 2 with the respectively nearest polygon or point from File 1.
I use QGIS. I'm almost sure there is an expression or method to do just that!
We also got another line layer that shows the distance between the points and the respectively nearest point/polygon (File 1), so I could technically color it all manually, but I want to learn how to automate my workflow.

EDIT 1: @Babel’s answer does not work for me so far.
Here are the wrong symbols:
https://ibb.co/ZXSkk2y
My two layers do not have the column "Id" but @osmId in the attribute table instead, but changing the "Id" with "@osmID" doesnt work either :/ https://ibb.co/020zD61
This is my command, which should be correct...: https://ibb.co/Fwy10Xk
I also tried to activate the current expression here, didnt work either: https://ibb.co/Y3bz13p
And lastly, I should add that that file 1 is a multi-polygon file with polygons as well as points...maybe theres the error?

EDIT 2: Those are my values from Layer/File 1 (points): https://ibb.co/r2wh0MX
and those are my vaules from Layer/File 1 (polygons): https://ibb.co/Df3dfgJ
Both definited by @osmId because there is no column "Id"
Now I want that the points from Layer / File 2 (Kindergarten_Karlsruhe) have the nearest colour either from the point from Layer/File 1 (Krankenhaeuser_Karlsruhe) (points) or from Layer/File 1 (Krankenhaueser_Karlsruhe) (polygons).
Layer 2 (Kindergarten_Karlsruhe) should receive the colour from either of those two geometries (polygons/points) from Layer 1 (Krankenhaeuser_Karlsruhe) which the respecitve kindergarten point is closest to.
Here you can see both Layers again: https://ibb.co/n3T03m2
Now, if type in the formula instead of "Id" "@osmId", there is an error messing appearing... see also: https://ibb.co/12jMZV2

EDIT 3:
This is correct, I have two (grouped) layers, both named Krankenhaeuser_Karlsruhe, once as polygon and once as points. I renamed both of them and set the colour settings from both again and tried to type in the command - still with an error!
See: https://ibb.co/4ZS08qQ
Error Message: https://ibb.co/vYSGNxn


texts: Are there legitimate companies that sell old versions of map info.  I use v9 and it will not load onto a new pc.  I don't want to buy the latest version as I fear it will be unrecognisable


texts: Clip linestring layer by polygons (partly overlapping)

texts:
I have a linestring layer that I want to clip with a polygon layer. So far so good. However, as these polygons are river catchments starting from a sampling site, they are sometimes overlapping (due to multiple sampling sites on a river). When this is the case, I want the clipped linestring to be duplicated according to the number of overlapping polygons (i.e. I want to have a linestring for each polygon within each polygon).

To illustrate you my example I've uploaded two sample shapes here. It contains 3 overlapping polygons and a linestring. My actual data consists of ~ 2500 polygons and a stream network for whole Germany.

One possibility that comes to my mind would be to loop through all the polygons' unique ID and clip them individually with the line network. However I ask myself, if there are already functions in QGIS, PostGIS, or R that can do this?

using the loop strategy in R I came with the code below.
It's basically working, but I still have problems when a line has multiple segments within a polygon and it's extremely slow compared to the clip function in QGIS for example.

require(sf)
ezg_fin = st_read('polygons.shp')
gew_typ = st_read('line.shp')

todo = ezg_fin$site_id

line_list = list()
for (i in 1:length(todo)) {
  catchm = ezg_fin[ezg_fin$site_id == todo[i], ]
  line_sub = gew_typ[catchm, ] # 1st spatial subset the stream segment (fast!) - only possible 'cause stream network is segmented
  message(i, '. ', 'Clipping river at: ', todo[i])
  line_clip = st_intersection(line_sub, catchm) # 2nd clip only the subsetted segement (whole stream network would take way too long!)

  line_list[[i]] = line_clip
}

fin = do.call(rbind, line_list)


texts: Split and count overlapping features within a layer

texts:
I would like to split a shapefile in qgis according to how the polygons overlaps.

See pic for explanation:


texts: Generating Queried CSV from shapefile using ArcGIS API for Python?

texts:
I'm trying to create a csv from a query of a feature layer hosted in ArcGIS Online. The csv only needs to contain one field, where the query is met, geometry is not needed. I'm having issues getting the queried data moved into a dataframe and exported into a csv. From what I've read, you can use the SpatiallyEnabledDataframe to query and export data from a feature, but I'm getting stuck generating the dataframe.
from arcgis.gis import GIS
from arcgis.features import FeatureLayer
from arcgis.features import manage_data
from datetime import date, time
import pandas as pd
from arcgis.features import SpatialDataFrame

gis = GIS("https://------", username="---------", password="---------")
item_id='---------'
get = gis.content.get(item_id)
feeder = get.layers[0]

sdf = pd.DataFrame.spatial.from_layer(feeder)
df = feeder.query(where="feeder_status = 'Closed'").sdf
df.head()

This throws the following error: AttributeError: type object 'DataFrame' has no attribute 'spatial'
Is this error thrown because pandas is looking for the features geometry?


texts: Crop two rasters to common extent where pixelvalues are not nodata

texts:
I hope this doesn't get too confusing. I have two raster that overlap in a certain region. Like the following:

What I would like to have, are two raster that only have values in the region where both rasters intersect, but also have values and no NA. Like the extent in the (really poorly drawn) image below:

What I did so far was using the gdaltindex command and compute a bounding box for each raster and intersect the, but than I get values for the red raster in the blue parts (and there I would like to have nodata in both resulting rasters)

I hope the explanation is sufficiently clear. Sorry about the terrible images...
The solution could be on the command line with gdal, qgis, with rasterio in python, something in R....

Thanks a lot:)


texts: Is the data supplied to Google Maps or leaflet.js javascript rest API is private

texts:
I have a dataset which has to be used for internal/intranet web application. Is the data supplied to google map javascript api to draw map as well as to show meta information through markers is private (data used intenally only ) and does not go to google server. Example : on mouse over newyork city marker/circles/point it should show private meta information .This info i dont want to be sent to public webservice/javascripts.I understand i need to send lat long geo information which is ok but dont want my custom data to travel to public and comeback. Does google map api safe for my usecase.

else i would use something like d3.js maps which assures me that all the data used are within the internal scope.


texts: Calculating highest point, lowest point, average altitude in specified area around click position in Cesium?

texts: I have successfully computed the height of a location where user clicks in cesium but my next step is a bit complex. What I want to do is that when user clicks at some point on globe, at any zoom level, I want to calculate highest point, lowest point and, average altitude around click position (say a circular area of 1000 meters).


texts: How to fix incorrect coordinate data for the mouse's position

texts:
I am running Quantum GIS 1.8.0 on a 64-bit Windows platform. I've imported a data set that plots correctly when overlaid with a Google Satellite map, so I know that the layer/project CRS is okay (since my points are being projected correctly). The issue is that when I hover my cursor over any point on the map, the reported coordinates are entirely different from the coordinates used for plotting the points. How do I fix this or where can I get information on how to resolve this? I've checked the QGIS bug repository, searched this forum, and reviewed the User Guide.

example

plotted point:
-76.83617103
39.21012998

map coordinates at mouse cursor position (same point as above):
-8553363.45
4751815.77

Thanks!


texts:
I need to change the measuring interval ("resolution") of tracks in gpx files. The reason is that I want to try a QGIS plugin that requires a specific interval.

I took a look at gpsbabel but from what I can see it only allows simplification by setting max number of points. I need to set points per second.

Any suggestions on how to achieve this?


texts: Find all of the ways that intersect a linestring

texts:
I'm new to Open Street Map and PostGIS, and am trying to figure out how to
find all of the ways that intersect a linestring (i.e. all of the streets that
intersect a route).  Is there a way to do this that doesn't involve calling
ST_Intersects on every row of the table?


texts: Calculating percentages in ArcGIS field calculator gives TypeError: 'float' object is not iterable?

texts:
I am trying to calculate the percent AREA in each class.
I have created a new percentage field and have tried to populate it using the following equation in python in the calculate field tool.

percent = float(!Area_ha!)/ sum(!Count!)

in other words I want to take the area of each record and divide it by the total area.

I am getting the following error

  File "expression", line 1, in module

TypeError: 'float' object is not iterable


texts:
How can I overlay the country shapes from this Leaflet example onto a WebGL Earth globe?
I've tried replacing geojson = L.geoJson(euCountries).addTo(map); with geojson = WE.geoJson(euCountries).addTo(earth);.
This produces: TypeError: WE.geoJson is not a function
I am adapting this code, importing the data and making the change listed above, to produce this:
<!DOCTYPE HTML>

    <title>WebGL Earth API: Markers</title>


texts: Removing black lines between merged tiles?

texts:
I merged several rasters, but there is a problem with lines between sheets.

Raster is created from LiDAR datasets. Rasters were produced from different datasets via this application: https://iaps.zrc-sazu.si/en/rvt#v and then merged.

For merging I used Raster > Miscellaneous > Merge in QGis 2.18.15.

Please see snapshot of a problem:

And the Merge function pop-up:

I found that even between unmerged sheets there is some minimal space:

How do I get rid of those lines?

Probably the original sin is already in the input.


texts: QGIS/PostGIS finding shortest route going through all the polylines at least once

texts:
I have a QGIS 3.10 map with two main PostGIS 3.0 layers:

Zones (Polygons): the town is divided in 26 zones. Each municipal street sweeper only works inside his zone.
Street lines (Linestrings): These are the paths municipal street sweepers follow when moving around the town.

I would like to get for each zone  (each street sweeper), which is the most optimized path in distance (and thus in time). There are several routing alternatives and I need to get the most efficient one, the one that following it would walk the lowest distance. The only condition is that it has to walk through all the lines at least once.
I've only found algorithms for finding shortest path from point to point but nothing similar to what I want...
Is there any QGIS functionality, plugin or PostGIS query that would help me with this task?


texts: select living-end-streets from a local osm-road-extract with QGIS

texts:
I am looking for a method or query to select living-end-streets from a local osm-road-extract.

Living-end-streets are dead-end-streets for cars but cyclists and pedestrians can continue on their way.

So I want to select roads that share their endpoint with one (ore more) "highway"  IN  ( 'footway' , 'footpath' ,'steps', 'cycleway','pedestrian','track','bridleway') but NO other highway types.

Preconditions are:

Preferrable Using built-in functionalities or plugins from QGIS 2.14 or later
No PostGIS
Only python if no other solution is possible - like it seems to be

Eg:


texts: Extracting OpenStreetMap data for area from 2 neighboring countries?

texts:
I need to prepare offline background data for my area of interest which is located around border between Czech Republic and Germany. I need to convert the OSM data to spatialite and then to use in QGIS as Google Maps like background.

I already have working method (from @Underdark) for how to convert OSM to spatialite DB and also the visual styles.

However I would like to avoid downloading OSM data for whole Europe - from Geofabrik.de it is about 21 GB.

Is it better to first somehow merge OSM files of Germany and Czech Republic or is it easier to perform Spatialite extraction (so would have smaller file and less data than for whole country) and then merge the Spatialite datasets?


texts:
On Windows 7, QGIS 2.12 Server returns no features on a GetFeature WFS request while WFS DescribeFeatureType request and WMS requests are OK.

Does anybody have similar problems with QGIS server 2.12 ?


texts: In OpenLayers 3, transformExtent to EPSG:8155 produces error

texts:
I'm trying to get the extent of a map using OpenLayers 3, so that I can use it in a CQL filter to limit the area of the WMS request to the area currently shown in the map extent. To do this, I have to transform the extent into the source layer's SRID coordinates, but I am getting an error when I try this. The code I'm using is:

var currentExtent = map.getView().calculateExtent(map.getSize());
transformedExtent = ol.proj.transformExtent(currentExtent,'EPSG:3857','EPSG:8155');

I get the following error bubbling up in Firebug, which points to line 68 in http://openlayers.org/en/v3.11.1/build/ol.js:

  TypeError: c is undefined function Je(b,c){var d=b.a,e=c.a,f;d in
  Ke&&e in Ke[d]&&(f=Ke[d][e]);void 0===f&...

If I transform from 3857 to 4326 using the same syntax, otherwise, it works as expected, giving me coordinates that make sense for 4326:

transformedExtent = ol.proj.transformExtent(currentExtent,'EPSG:3857','EPSG:4326');

8155 is the WYLAM SRID, and was added somewhat recently.  It looks like some of the formatting information is missing on spatialreference.org.

Any ideas as to whether this is a problem that needs to be addressed in OpenLayers, Spatial Reference, or elsewhere?


texts: Clipping Feature by using irregular polygon in ArcObjects?

texts:
I want to clip the features based on selected polygon.I have the two layers "suburban" and "lot". I want to clip features of lot layer based on selected polygon from suburban layer.
I have tried manually in arc map 10.3 using Analaysis tools->Extract->Clip. It clipped lot features only which are underlying on selected polygons of suburban layer. But it's not working correctly when trying the same from ArcObjects .Net using geoprocessor interface.

It's not clipping the lot feature based on selected polygon of suburban layer.Instead of clipping the lot features based on all the polygon of the suburban layer.

Expected output:

I am also tried ITopologicalOperator QueryClipped functionality.It clip the features only in rectangular (envelop) size.It's not exactly clip the features based on irregular polygon boundary.

I am using the following code.

Private Sub lotClip(ByVal pInFClass As ESRI.ArcGIS.Geodatabase.IFeatureClass, ByVal clipGeometry As IGeometry, ByVal pClipFClass As ESRI.ArcGIS.Geodatabase.IFeatureClass, ByVal outputName As String)
        Dim clipperGeo As New ESRI.ArcGIS.AnalysisTools.Clip
        Dim Result As ESRI.ArcGIS.Geoprocessing.IGeoProcessorResult2
        Try
            Using releaser As New ESRI.ArcGIS.ADF.ComReleaser
                releaser.ManageLifetime(clipperGeo)
                clipperGeo.in_features = pInFClass
                clipperGeo.clip_features = clipGeometry
                clipperGeo.out_feature_class = outputName
                Result = runClipTool(clipperGeo, Nothing)
                If Result Is Nothing Then
                    System.Windows.Forms.MessageBox.Show("Could not clip dataset")
                    Exit Sub
                End If
                ReturnObjectfromResult(Result, "Feature Class")
            End Using
        Catch ex As Exception
            System.Windows.Forms.MessageBox.Show(ex.ToString, "Clip error")
        End Try
    End Sub
    Friend Function runClipTool(ByVal Process As ESRI.ArcGIS.Geoprocessor.IGPProcess, ByVal TC As ESRI.ArcGIS.esriSystem.ITrackCancel2) As ESRI.ArcGIS.Geoprocessing.IGeoProcessorResult2
        Dim Result As ESRI.ArcGIS.Geoprocessing.IGeoProcessorResult2
        Dim GP As New ESRI.ArcGIS.Geoprocessor.Geoprocessor
        Try
            Result = CType(GP.Execute(Process, Nothing), ESRI.ArcGIS.Geoprocessing.IGeoProcessorResult2)
            If Result.Status <> ESRI.ArcGIS.esriSystem.esriJobStatus.esriJobSucceeded Then returnMessages(Result, "Geoprocessing Error")
            GP.ClearMessages()
        Catch ex As Exception
            returnMessages(Result, "Fail")
            System.Windows.Forms.MessageBox.Show(ex.ToString, "Run Geoprocessor")
        End Try
        Return Result
    End Function
    Friend Sub ReturnObjectfromResult(ByVal result As ESRI.ArcGIS.Geoprocessing.IGeoProcessorResult2, ByVal Type As String)
        Dim GPVal As ESRI.ArcGIS.Geodatabase.IGPValue
        Dim InMemFC As String
        Dim GPUtil As ESRI.ArcGIS.Geoprocessing.IGPUtilities3 = New ESRI.ArcGIS.Geoprocessing.GPUtilities
        Dim pOutputFc As ESRI.ArcGIS.Geodatabase.IFeatureClass
        Try
            GPVal = result.GetOutput(0)
            InMemFC = GPVal.GetAsText()

            Select Case Type
                Case "Feature Class"
                    pOutputFc = GPUtil.OpenFeatureClassFromString(InMemFC)
                Case "Table"
                    GPUtil.OpenTableFromString(InMemFC)
            End Select
        Catch ex As Exception
            System.Windows.Forms.MessageBox.Show(ex.ToString, "Return FeatureClass error")
        End Try

    End Sub

How do I clip the features based on an irregular polygon?


texts: Installing ArcGIS Desktop on web server?

texts:
I am working on a spatial tool that we do not want to publish as a geoprocessing service due to certain organizational needs; instead we are publishing up a .NET MVC application that calls a Python script. This Python script utilizes ArcPy to do some spatial analysis and thus it needs ArcGIS Desktop or ArcGIS Server to run. Since the .NET application sits on the web server once it is published, that is also where ArcGIS Desktop would need to be installed so the script can run.

Aside from probably not being best practice, can anyone forsee any problems with this?


texts: Raster image outside Grass screen, How to get it inside from start?

texts:
I'm new in Grass and need some help!

I have created a LOCATION based on the coordinate system in which I want my data to be projected. So far I'm being able to Import & display the raster image (Tiff Image). However, to see this image have to click several times in zoom-in/out, I understand this is because the data is somewhere out of the "screen" but still inside the LOCATION. Is there to get this image inside de screen from the start?

Location created:
proj=UTM
datum=wgs84
zone=30
towgs84=0,0,0,0,0
no_defs=defined

The Geo-TIFF image (MOD13Q1-NDVI and MOD13Q1-EVI ): was reprojected using MRT-Modis:
-UTM
-Resampling: nearest Neighbor
-utm zone: 30
-datum=wgs84
-UL 39.9999 / -13.05407
-LR 29.9999/ 0.0

Despite creating few locations to get this data "inside the screen", couldn't make it. Don't know where my mistake is.


texts:
I am new to GIS and am trying to create my first map.

Could anyone recomend any good tutorials which are clear and simple? I am having difficulty importing my open source map information into my QGIS propgram. (I have had a look at the tutorial and online)


texts: How to toggle ArcGIS webmap layers using ArcGIS Javascript API

texts:
I've created a webmap using ArcGIS Online that I want to render in a web app using the ArcGIS for Javascript API.  In ArcGIS Online, I've saved the webmap with all of the layers off (visibility: false).  I want set the visibility of one layer ("Layer X" below) to true using the API.

    arcgisUtils.createMap(mapid, "map").then(function(response){

         layers = response.itemInfo.itemData.operationalLayers;

         dojo.forEach(layers, function(layer){

             if (layer.title == 'Layer X') {
                 layer.setVisibility(true)**;
             }
             console.log(layer.title + ", visibility: " + layer.visibility);

         });

     });

**I've tried using .show() as well.

Using the code snippet above, the visibility of "Layer X" doesn't change.  The console log reports the layer title & visibility of all the layers except the one I'm trying to show.

The response object returned from the .createMap method is a Deferred object.

Is it possible to toggle the layers of a webmap using the ArcGIS for Javascript API, or is the response object immutable?


texts: Options for precise relative position calculations

texts:
I'm investigating an industrial application that has a requirement to calculate with precision the relative distance on the ground between two outdoor points in the following scenario:

The first point will vary in height from 10m to 150m off the ground
The second point will always be around 3m - 5m off the ground
The points will be separated by 0m to 3m on the ground plane
I'm only interested in the separation on the ground, not the vertical separation
The precision should be within 10cm and I don't care about the absolute position accuracy.
Both points could be slowly moving around any axis before stabilising.

Looking into options available, I'm thinking standard GPS should be able to give me this type of precision if I compute the relative differences in position, which should hopefully compensate for the absolute positioning errors that GPS normally has, so I can avoid the complexity / cost of DGPS system.

I'm unsure if the differences in vertical height between the 2 points will impact the precision, and I also assume I have to use the a quality GPS receiver of the same type for both positions (to remove any differences caused by the receiver).

If this won't work, do I have any other options for this use-case (eg. combining an IMU with the GPS data, DGPS)?


texts: How to a sync file geodatabase across multiple office locations?

texts:
I am planning to have a work project with multiple file geodatabases, fgdb (consisting of images, grids, feature class, etc) and these will be passed along to three regional offices and will sit on their individual file servers (all servers are linked via distributor file system). And I was told that I can set them up to sync up during the after office hours.

There will be a "parent" fgdb which is sitting in my office and the rest are "child" databases. Of course, this will bring more complex questions, but allow me to ask below questions first.

Do I need any add-on to do it?
If no, where do I set this up?
Will Symantec back up which run at night be affected or affecting this sync?
Is there a way to bring non-ESRI files into fgdb ? (eg. .xls, .dat etc)

FYI. All users, are on 9.3.1 , with ArcMap only, without other ESRI products.


texts: How can we get the soil type, temperature and rain from API 2.0?

texts:
We were using the v1 of the SoilGrids REST API and previously we were able to get the soil type, rain and temperature. We used the following values:
For the API v1, the query required only the latitude and longitude of the location we wanted to know about.
GET https://rest.soilgrids.org/query?\
     lat=xxxxx&lon=yyyyy\
     &attributes=TAXOUSDA,TMDMOD,PREMRG\

And then I used the response to pull the values of soil type, temperature, and precipitation.
$soil_properties[ 'soil_type' ] = $result->properties->TAXGOUSDAMajor;
$soil_properties[ 'temp' ] = $result->properties->TMDMOD_2011->M
$soil_properties[ 'rain' ] = $result->properties->PREMRG->M

Now the API v1 is returning a 404. The new version (API 2.0) doesn't appear to provide these values.
Where can I get them?


texts: Annotation stored in a database is not showing my leader lines

texts:
I am using 10.2.2, I initially tried to use just plain annotation to label my features. This however was not useful as it was very laggy and took at least 30 sec or more to choose each individual annotation. ESRI help did not have a clue as to why this was happening. We tested and showed that storing the annotation in a database was much quicker. This being said, I really need it to have leader
lines, when I initially labeled the features I did make the label the way I liked with the leader lines, once it got into the annotation
format in the database the leader lines mostly disappear and the ones that show up disappear when I adjust the label.

I am using Microsoft Windows 7 Enterprise

Any clue to why this is happening would be wonderful! This map was needed last week!

Thanks!


texts: Get more than 2 aggregate results in QGIS

texts:
I have an indexed grid and am using the following in field calculator to get a grid reference for each grid a road intersects.

It works fine except when the line intersects more than 2 grid cells. Is there a way to write this so that it finds all grid intersects?

aggregate(
layer:= 'Grid',
aggregate:='concatenate',
expression:=Name,
concatenator:=', ',
filter:= intersects(
geometry(@parent),$geometry)
)

So the road in the image below should have Grid_Ref = AC13,AD13,AE13,AE14,AE15

Also is it possible virtually merge the dataset this is being run on so that it treats all segments that have the same road name as one?


texts: Select layer by location in ArcPy ERROR 000368

texts:
I am trying to create a scrip that create a list of all features in my GDB, then select objects that intersect with a user given shapefile, and finally copy selected features to a new location. After I have successfully created a list with features within my GDB, I am having troubles selecting objects that intersect with the shapefile provided by the user.
My code looks like this:
import arcpy
import os
from arcpy import env

nuevaAOI = arcpy.GetParameterAsText(0)

pathBase = r"C:\Users\SCGEO\Documents\SIANCERO\InfoTematicaEIA.gdb"
arcpy.env.workspace = pathBase

featureDataSets = arcpy.ListDatasets(feature_type='feature')

featureClasses = []

for fds in featureDataSets:
    for fc in arcpy.ListFeatureClasses(feature_dataset=fds):
        featureClasses.append(os.path.join(pathBase, fds, fc))

for f in featureClasses:
    arcpy.SelectLayerByLocation_management(f,"INTERSECT", nuevaAOI, 2000, "NEW_SELECTION")

After I run the code an error message is prompt stating:
"Parameters are not valid.
ERROR 000368: Invalid input data.
Failed to execute (SelectLayerByLocation)."

Can anybody give me a hint on this?


texts: Writing PDF anchor point into text element of layout using ArcPy for ArcGIS Pro?

texts:
I am currently working on an application that uses ReportLab to assemble PDF reports that include maps and tables (both created using ArcMap and ArcPy) and that application may soon be converted to use ArcGIS Pro and ArcPy instead.

As powerful as ReportLab is, much of the functionality that we use from it can be achieved using ArcGIS Pro and ArcPy, and it would simplify our deployment if ReportLab could be eliminated.

The requirement that I cannot meet without it is to embed a hyperlink in one single page PDF file that references an anchor point on another single page PDF file, and to have the hyperlink work when those two files are appended using appendPages on the PDFDocument object of ArcPy for ArcGIS Pro.

For the first step of this workflow, is it possible to write text into a Text element of one layout that can be used as an anchor point (e.g. Fig.1 is here) and have that anchor point remain usable when the layout is exported as a PDF file?

Changing the text in a text element is straightforward from ArcPy but it is the formatting to use within that text element which I seek, if it exists.


texts: Projecting EPSG:4326 with OpenLayers 2 and GeoJSON?

texts:
I'm trying to project some GeoJSON data but it ends up at 0,0. I've looked at other posts and tried those solutions to no avail.

Can anyone help?

Here is the code:

  <title>OpenLayers Example</title>


texts: How do you add a long list of rasters together in Python?

texts:
If you have produced a long list of rasters (lets say 50) and want to add them together, is there any shorter way to write in python rather than typing all the inputs, like so:

outRas = inRas1 + inRas2.........+ inRas50


texts: Finding problematic rasters in large set of rasters using ArcPy?

texts:
Just wondering if anyone knows how to quickly find rasters that cause errors when I attempt to process them in cell statistics (as I am trying to sum them). I have many hundreds of rasters and don't want to trial putting them 1 by 1 into ArcMap to see which one makes it crash.

I am unfamiliar with python, so am unsure of how to approach this.
I have tried this python code within arcmap, however nothing appears after I hit enter

import arcpy
from arcpy import env
# Set the current workspace
arcpy.env.workspace = "F:\Uni_Work\WOFS data all folders"
# Get and print a list of GRIDs from the workspace
rasterlist = arcpy.ListRasters("*", "All")
for raster in rasterlist:
    print rasterlist


texts: Display colour coded road traffic data in QGIS

texts:
I would like to colour code a shapefile containing road polylines to represent the heaviest flows of traffic. The shapefile contains attribute data linked to each section of road where road traffic counts were taken (count points). There are two fields that contain Road names and a unique count point ID.

The second shapefile contains point geometry representing the location of each count point. The attribute data contains various fields including traffic numbers along each section of road, road name and a corresponding count point ID.

How do I A) best combine the two sets of attributes to then be able to B) colour code the busiest roads in the darkest red?


texts:
I have a SQL-Layer that is defined by the following statement:

select b.id,
       b.description,
       b.is_event,
       geodesic_buffer(b.location_point, 3000)
from bus_stop b
where b.event_id = 7

Is there any way to add a parameter to the SQL, so that it is easier to change, like where b.event_id = :event_id


texts: GDALwarp converting data from 'Planar' to 'Geographic (degree)'

texts:
I have a dataset in nc format (Daymet weather data) using gdal translate. I convert them to GeoTIFF, then using GDALwarp I try to change the coordinate system type to geographic because I want to have latitude and longitude. It does not give me the right lat and lon. I am not sure what I am doing wrong

gdal_translate -of GTiff -a_srs "+proj=aea +datum=WGS72 +lat_1=25 n +lat_2=60n +lat_0=42.5n +lon_0=100w" NETCDF:"prcp_11742_2015.nc":prcp prcp_11742_2015_aea_lat.tif

then

gdalwarp -t_srs EPSG:4326 prcp_11742_2015_aea_lat.tif prcp_11742_2015_aea_lat001.tif

the source file is in this format

The North American Daymet projection system and parameters:
Projection System:  Lambert Conformal Conic
Parameters:
projection units:   meters
datum (spheroid):   WGS_84
1st standard parallel:  25 deg N
2nd standard parallel:  60 deg N
Central meridian:   100 deg W
Latitude of origin: 42.5 deg N
false easting:      0
false northing:     0

The lcc is the projection of source file and i I wanted to change the projection to the equal distance area too that is why I used the proj =aea, also one problem i have is i want to have the lat/lon in the geotiff file but i do not see the right lat/lon in the geotiff file. the corners are as follow but after the conversion

Corner Coordinates:
Upper Left  (  158000.000,  -45000.000) ( 97d59'54.18"W, 42d 3'26.22"N)
Lower Left  (  158000.000, -263000.000) ( 98d 3'45.07"W, 40d 0' 2.98"N)
Upper Right (  326000.000,  -45000.000) ( 95d52'21.59"W, 41d59'56.50"N)
Lower Right (  326000.000, -263000.000) ( 96d 0'17.11"W, 39d56'40.19"N)
Center      (  242000.000, -154000.000) ( 96d59' 5.65"W, 41d 0'18.37"N)

i have this

Corner Coordinates:
Upper Left  (-126.0746130,  -0.4065400) (126d 4'28.61"W,  0d24'23.54"S)
Lower Left  (-126.0746130,  -2.3815368) (126d 4'28.61"W,  2d22'53.53"S)
Upper Right (-124.5595470,  -0.4065400) (124d33'34.37"W,  0d24'23.54"S)
Lower Right (-124.5595470,  -2.3815368) (124d33'34.37"W,  2d22'53.53"S)
Center      (-125.3170800,  -1.3940384) (125d19' 1.49"W,  1d23'38.54"S)


texts:
I'm beginner in python and learning it with examples. I have 100s of dbf tables like

    Type value1 value2           Type value1 value2
       A      3      1    and       A      5      7
       B      4      2              A      6      8

I want them to simply merged in single file like

    Type value1 value2
       A      3      1
       B      4      2
       A      5      7
       B      6      8

to get this from arcpy I used the code as below but it doesn't return any results or any errors. I have made "alldata.dbf" at that place with exact column to work with but the updateCursor doesn't update anything on that file. Can anyone help? or is it blunder to do thing like this in python instead of other easier method?

import arcpy

arcpy.env.workspace = r'D:/dbflst/'
dbfList = arcpy.ListTables()

finalTable = r'D:/alldata.dbf'

for dbf in dbfList:
dbfFile = arcpy.env.workspace + '/' + dbf
dbfData = {}

# read in data from DBF and create dictionary
with arcpy.da.SearchCursor(dbfFile, ['Type','Value1','Value2']) as cursor:
    for row in cursor:
        dbfData[row[0]] = row[1]
del cursor

# write data from dictionary to final table
with arcpy.da.UpdateCursor(finalTable, ['Type','Value1','Value2']) as cursor:
    for row in cursor:
        row[1] = dbfData[row[0]]
        cursor.updateRow(row)
del cursor


texts: gdal2tiles.py gives ERROR 6 about EPSG:900913 on fresh Ubuntu 14.04 install

texts:
I have a fresh install of Ubuntu 14.04. I installed gdal with sudo apt-get install gdal-bin python-gdal. After this I have: GDAL 1.10.1 released 2013/08/26.

I am trying to tile the weather radar images from NOAA. The following are my steps and the error message that is displayed.

$wget http://radar.weather.gov/ridge/Conus/RadarImg/latest_radaronly.gif
$wget http://radar.weather.gov/ridge/Conus/RadarImg/latest_radaronly.gfw
$gdal_translate -of vrt -expand rgba ./latest_radaronly.gif temp.vrt
$cp ./latest_radaronly.gfw ./temp.wld
$gdal2tiles.py -s EPSG:3857 -w none -r near -z 3-5 -n ./temp.vrt ./temp/

The errors shown are then:

ERROR 6: EPSG PCS/GCS code 900913 not found in EPSG support files.  Is this
a valid EPSG coordinate system?
ERROR 6: No translation for an empty SRS to PROJ.4 format is known.
ERROR 6: No translation for an empty SRS to PROJ.4 format is known.
Generating Base Tiles:
ERROR 5: Illegal values for buffer size
ERROR 5: Illegal values for buffer size

From that error, it seems that there is a problem with the EPSG:900913 (google) format. That is why I specifically set -s EPSG:3857.

All of this should work. I did the steps above from scratch, because I have an older server (ubuntu 12.04) where I am currently doing this without error. Once I upgrade to 14.04 the error above shows up. I am trying find the reason for the error. A google search for the the error message results in some old email list traffic that were either not about ubuntu or had other issues.

As a side note, I also tried building from source and there was no change in the errors displayed.


texts: PyQGIS Export table with Replace all selected raw field values by displayed values option

texts:
I have a table with multiple fields with widget types configured (value map, checkbox, date/time).
What I need is to export this table in a geopackage, using this in a PyQGIS script, using the "Replace all selected raw field values by displayed values".
rutaBase = QgsProject.instance().absolutePath()
rutaGpkg = rutaBase + '/geopackage.gpkg'
table = QgsProject.instance().mapLayersByName('tblActuaciones')[0]

opt = QgsVectorFileWriter.SaveVectorOptions()
opt.layerName = 'exported_table'
opt.actionOnExistingFile = QgsVectorFileWriter.CreateOrOverwriteLayer
opt.fieldValueConverter = QgsVectorFileWriter.FieldValueConverter
QgsVectorFileWriter.writeAsVectorFormat(table, rutaGpkg, opt)

But I get this error:
TypeError: sip.wrappertype cannot be converted to qgis._core.FieldValueConverter in this context


texts: Trying to convert to EPSG:27700 (Ordnance Survey) projection, unable to do so

texts:
I'm trying to code an interface that allows a user to draw a vector shape on openstreetmap, which is then converted to a WKT string, processed by the server and sent of as a WFS query to emapsite (which provide the Ordnance Survey data, including a veru detailed vector layer of the UK). However the issue I'm having is that when I draw a polygon around all of the UK (which should obviously return data), I'm not getting data returned (as in the server says there are no vectors there). However when I draw a polygon in the south-west corner of Africa (just off the coast of Nigeria), I'm getting a huge amount of data returned. Ordnance Survey uses EPSG:27700 (OSM uses EPSG:90013), which I don't really understand, but I imagine that I'm not doing projection conversions that are necessary. code below:

http://pastebin.com/wi4CUAmc


texts:
I need to load in a docked QWebView this html:

    < html >
< body >
< iframe width="100%" height="499" frameborder="0"
scrolling="no" marginheight="0"
marginwidth="0"
src="http://dev.virtualearth.net/embeddedMap/v1/ajax/BirdseyeWithLabels?zoomLevel=10&pushpins=45.0610173308_7.64466465993&center=45.0610173308_7.64466465993&heading=123" >
< /iframe >
< /body >< /html >

I'm using this code:

class WebView(QWebView):
def __init__(self, parent=None,framed=False):
    QWebView.__init__(self, parent)
    #QObject.connect(self, SIGNAL("loadFinished(bool)"), self.onLoad)
    #QObject.connect(self.page().networkAccessManager(), SIGNAL("sslErrors(QNetworkReply *, const QList<QSslError> &)"), self.onErrors)
    self.settings().setAttribute(QWebSettings.PluginsEnabled, True)
    self.framed = framed

def load(self,url, params):
    mqurl = QUrl.fromUserInput(QString.fromAscii(url))
    p = []
    QgsMessageLog.logMessage("url="+url, "WebView")
    for k, v in params.iteritems():
        p.append( (k, v) )
        QgsMessageLog.logMessage(str(k)+'=' + str(v), "WebView")

    mqurl.setQueryItems(p)
    if self.framed==True:
        mhtml=HtmlFrame(mqurl.toString(),self.height()-20)
        qmhtml = QString.fromAscii(mhtml)
        self.setHtml(qmhtml,mqurl)
    else:
        QDesktopServices.openUrl(mqurl)

If framed is False it works, if framed is True the webpage is not loaded.
With other webpages it works.
Is it due to Ajax?
How to solve this?
Thanks


texts: OpenLayers: How do I keep a feature fixed while the map layer is moved (like in-car GPS displays)?

texts:
I am aiming to create something simple that works like a GPS navigation display, where the marker for 'you are here' is fixed while the map moves along. I have latitude and longitude coordinates supplied but when I move the map, the marker moves too - even though I tried setting the marker's layer to isFixed: true.

Here's what I'm using to set up the map:

var map = new OpenLayers.Map('map');
var baseLayer = new OpenLayers.Layer.MapServer('base layer', 'url', {map: 'something.map'}, {layers: 'base'}, {transitionEffect: 'resize'}, {gutter:'15'});

var markerLayer = new OpenLayers.Layer.Vector('marker', {isFixed : true});
var feature = new OpenLayers.Feature.Vector(
  new OpenLayers.Geometry.Point(100, 1.5),
  {some:'data'}, {externalGraphic: 'img/marker.png', graphicHeight:25, graphicWidth:30});
markerLayer.addFeatures(feature);

map.addLayer(baseLayer);
mapp.addLayer(markerLayer);
map.setCenter(new OpenLayers.LonLat(lon, lat), zoom);

A script at the bottom of my HTML updates lat and lon with coordinates (in degrees) pushed to it every second via a server-side event. After that, every second I am calling map.setCenter(new OpenLayers.LonLat(lon, lat), zoom) to re-center the map.

At first I thought that if I drew my feature on a new, fixed layer, it wouldn't move, but the feature seems to be being drawn at 100N 1.5E, and staying there despite me panning the map. I want the marker to just stay where it is no matter how you pan the map.

P.S. My map re-centering feels rather inefficient as I am creating a new LonLat every time, but I couldn't figure out from the API docs any way to directly update it e.g.:

center = new OpenLayers.LonLat(lon, lat);
map.setCenter(center, zoom);
FOR every second
  ... update lon and lat with SSE ...
  center = new OpenLayers.LonLat(lon, lat);
  map.setCenter(center, zoom);
END FOR

Any help or suggestions would be much appreciated! :)


texts: Finding differences between two shapefiles in Python

texts: I'm trying to find out whether there were any border changes between two yearly shapefiles. I've heard a lot of praise for Python's open source GIS tools like Shapely/Fiona, and thought I would try those, but I'm at a loss for how to go about implementing this. Is Shapely poorly suited for this task, and, if so, what tool would be better? Any advice you have would be greatly appreciated. Thanks!


texts: Trigger history_ How can I define the last updated column?

texts:
I am trying to implement a trigger function in PostgreSQL. I've written this trigger history function:

CREATE OR REPLACE FUNCTION history_update()
RETURNS trigger AS

$$

    BEGIN

    INSERT INTO data_history (id, geom, origine, descrizione, img, created, created_by)

    VALUES (NEW.id, new.geom, new.origine, new.descrizione, new.img, current_timestamp, current_user);

    RETURN NEW;

END;

$$

LANGUAGE plpgsql;

 CREATE TRIGGER update_trigger
 AFTER UPDATE ON data
 FOR EACH ROW EXECUTE PROCEDURE history_update();

Now I would like to add another column (called last_update) where to set the name of the column in which the value was updated.

For example if the updated value is in the column called "supporto", in the column  "last_updated" I want to see the name of the same column with the updated value (called "supporto").

What if the updated values in the column are more than one?

possible suggestion

1) Insert in the row of the column "last_updated" more values

or

2) Insert as many column as  the number of the column in the original table (for example "last_updated1", "last_updated2"...)

Any advice?


texts:
I am currently using ArcGIS online to create web maps and for mobile data collection; the data is stored/ created in a personal geodatabase, added to a .MXD and symbolized, etc. then pushed up to the web to create a web map.  It is hosted by ArcGIS.

I deal with a lot of users and a lot of data. This set up has had some significant limitations. For example, backing up data; because it is stored in a personal geodatabase in order to back it up I need to over write it each time.  Another limitation is that I can’t get live updates or syncs in the field again because it is stored in a personal geodatabase.

My question, what are alternatives to the system I am using?  I realize this answer may be involved so sending me in the right direction to find information or resources would be great.

This also needs to be within the ESRI environment.


texts: Maintaining transparency of layer during edit sessions of ArcMap?

texts: I'm editing the polygons of one shapefile to more closely match the polygons of another shapefile, essentially tracing the one over the other. I've set the transparency of the layer under edit to 75%, which allows me to more easily recognize errors and discrepancies compared to the other, more accurate shapefile. The problem is that while I have a polygon selected during an edit session, it changes the transparency to something closer to 20%, which makes the process much more difficult. Is there a way to lock transparency during selection/edit sessions?