texts: Implementing WaterML2 output in my web service

texts:
I'm trying to implement WaterML2 output for my sensor network webservice. Currently I report sensor data where the data includes IDs for each sensor, time of sensor readings, units and, optionally, GPS data. Typical output looks like the following:

PackageID   TimeStamp   TimeZone    RelayID SensorID    SensorType  Units   Data    Longitude   Latitude    Altitude    Speed   Climb   Track   LongitudeError  LatitudeError   AltitudeError   SpeedError  ClimbError  TrackError
210 2014-04-23 20:36:45.728 UTC SumitTest01 FFF34E141400    Temperature C   24.125  11.574138333    50.924998333    228.2   0.0 0.0 10.57   17.626  16.966  48.22   33.27
209 2014-04-23 20:36:28.978 UTC SumitTest01 FFF34E141400    Temperature C   24.1875
208 2014-04-23 20:29:33.004 UTC SumitTest01 FFF34E141400    Temperature C   23.625
207 2014-04-23 20:29:14.503 UTC SumitTest01 FFF34E141400    Temperature C   23.625

SensorType can be Temperature, Humidity, pH with more to come as we add more sensors.

A few years ago we looked at implementing WaterML1, but the specification was confusing and we decided against it. I've been at talks since where researchers have really pushed WaterML2 and one person even told me how it fixed a lot of problems with the confusion in WaterML1.

I've been looking at the specification and examples. It looks pretty straightforward, but I can't seem to find any libraries that can handle the models and XML generation for me in Java (my web service is written in Scala and I can use any Java library).

I've found the following Maven artefact (sbt include):

"org.jvnet.ogc" % "waterml-v_2_0" % "2.2.0"

...but no documentation on it. I've looked through the jar file and it just seems to be automatically generated jaxb stubs.

So my question is, are there any Java/Scala libraries out there that can help me create WaterML2 using a simple API so I don't have to write my own generator in scala.xml/jdom/etc?

If not, are there better examples out there for simple data that only deals with readings (like Temperature) and fixed GPS data points? How would I deal with the accuracy measures in my sample set I get from the GPS? Can those be represented or should I just discard them?

TL;DR Are there any good Scala/Java libraries to help generate WaterML2 and, if not, where are more examples of creating WaterML2 for point/sensor data?


texts: Formatting date inside ArcGIS Online Indicator widget

texts:
In my data I have a field in "date" format. The data source is a survey123 layer.
I am not able to format the appearance to just year, month, day. In the widget the Time and pm or am is always shown.
Is there no way to format the date in the indicator widget?


texts: How to change the case of input text (xtype: "textfield")

texts:
I have an input box which asks users to input a UK postcode but which specifies that this must be in lower case, for example so16 0as.  Some users find this frustrating as they would just like to copy and paste a post code into the textfield; the copied postcode may contain upper case letters ie SO16 0AS.  This will return an error.

The input code is as follows:

items: [{
            xtype: "textfield",
            width: 220,
            name: "search_lc",
            value: "",
            emptyText:'Enter Postcode (e.g. so16 0as)',
            fieldLabel: "Find",
            allowblank: false
        }],

I have tried to use fieldStyle: 'text-transform:lowercase' both below the fieldLabel and above it but this doesn't work.

I'mm thinking that this should be a fairly simple thing to resolve but I am experiencing a 'can't see the wood for the trees moment'!  Any help gratefully received!

Many thanks,

Red.


texts: Fill polygons via data-defined pie charts

texts:
I am trying to summarize my polygons by each's distribution of some variable.

To couch the discussion, let's say we wanted to summarize the demographics of each county (of a small state) by it's demographic makeup--say, black/white/other.

What came to mind would be to treat the fill of each county's polygon as a pie chart, and slice it accordingly.

Is there a way to do this in base QGIS (or SAGA)?

If not, what about through scripting? I was thinking one might be able to split the polygon layer into three new layers--white, black & other--with one slice of each polygon corresponding to the right wedge size.

This sounds vaguely doable to me, but I've got nothing like this GIS chops it would take. Any pointers?

Here's a rough sketch of what I've got in mind for two counties (with made up data).

The give-up alternative is to use the pie charts labeling feature.
Which may be more readable in the end anyway... just wondering whether this is possible anyway.

Note that the areas of the resulting sub-polygons should correspond (in proportion to the overall area) to the data points--i.e., the B sub-polygon in Providence should be 10% of the Providence's area, and the B sub-polygon in Kent should be 50% of Kent's area.


texts: I found these peninsula very interesting.  Tracker showed ISS over the bottom portion of South America and the features resemble portions of the Falkland Islands but I'd love to explore more.  Any ideas?


texts: Extract the filename from a chunk when processing a LAScatalog

texts:
I'm using the R packages lidR and TreeLS to process lidar data and I made a function that takes in a single LAS file and splits out a shapefile, it looks like this:

lasfile = 'tile_1.las'
las = readLAS(lasfile)

stemdetection = function(las, ...)
{
  UseMethod("stemdetection", las)
}

stemdetection.LAS = function(las, vox, dens, n, rad)
{
  thin = tlsSample(las, voxelize(vox))
  map = treeMap(thin, map.hough(hmin = 0.25, hmax = 16, hstep = 0.5, pixel_size = 0.05, max_radius = rad, min_density = dens))
  tls = stemPoints(las, map)#classify stem points
  seg = stemSegmentation(tls, sgmt.ransac.circle(n))#extract measures
  cs <- crs(tls)
  thin <- NULL
  map <- NULL
  coordinates(seg) = cbind(seg$X, seg$Y)
  output <- writeOGR(seg, td, paste(tools::file_path_sans_ext(basename(lasfile)), "_stems", sep=""), driver="ESRI Shapefile", overwrite_layer=TRUE)
  return(output)
}

I'm trying to change the function to work on a LAScatalog but having some issues with getting the name of each las file (chunk). Here is one example of what I tried:

stemdetection.LAScluster = function(las, vox, dens, n, rad)
{
  las <- readLAS(las, select = 'xyzi')# Read the LAScluster
  if (is.empty(las)) return(NULL)              # Exit early (see documentation)

  thin = tlsSample(las, voxelize(vox))
  map = treeMap(thin, map.hough(hmin = 0.25, hmax = 16, hstep = 0.5, pixel_size = 0.05, max_radius = rad, min_density = dens))
  tls = stemPoints(las, map)#classify stem points
  seg = stemSegmentation(tls, sgmt.ransac.circle(n))#extract measures
  cs <- crs(tls)
  thin <- NULL; map <- NULL ; tls = NULL
  coordinates(seg) = cbind(seg$X, seg$Y)
  output <- writeOGR(seg, td, paste(tools::file_path_sans_ext(basename(ctg@data$File)), "_stems", sep=""), driver="ESRI Shapefile", overwrite_layer=TRUE)
  return(output)
}

stemdetection.LAScatalog = function(las, vox, dens, n, rad)
{
  catalog_apply(las, stemdetection, vox = vox, dens = dens, n = n, rad = rad)
}

Is there a way to append the filename to the LAS object while (or after) performing the readLAS function?

Thanks in advance!


texts: Get JSON WebMap from ArcGIS JavaScript API Map object

texts:
I'm trying to get a WebMap object (as JSON) from a JavaScript Map object in the ArcGIS JavaScript API. Is there any way to do this within the API, without using ArcGIS.com? Ideally something like:

var webMapAsJSON = map.toWebMap();

From the "Export Web Map Task" documentation in the REST API, there's this line that suggests it should exist:

  "The ArcGIS web APIs (for JavaScript, Flex, Silverlight, etc.) allow
  developers to easily get this JSON string from the map."

However, I don't see anything in the Map object or elsewhere in the API that would do this. Any ideas?


texts:
MY apologizes about the next doubt, but I can't understand what is an orphan point. I am using QGIS with a module to apply the Linear Reference System (the module is called "LRS"). During the Calibration step, I select the Linear layer and Point layer. After I go to "Errors" step and the module identifies the next errors:

  Not enough points
  Orphan point

The question is: What is an orphan point?

Solved: The problem was solved. The id's of each layer must have the same ID.


texts: I have a 3D building done using Revit , i have exported the file using FME extension into IFC file in order to make the integration with ArcGIS ! i didn't find the way to make this happen , the only way i have discovered is to export the IFC file to KML. so any one got an idea the process of the integration BIM & ArcGIS.


texts: Currently, the formula I'm using to convert MGRS to Lat/Long and vice versa is giving me slightly different results than most other conversion utilities.  It's usually between 10 to 30 feet off.  This isn't a very big difference, but it matters for the particular application of data.  I'm wondering if anyone knows of a simple library that does conversions between MGRS and Lat/Long or UTM.


texts: How do you get the coordinates from a drag event in the Google Maps API V3?

texts: I have made a Google Version 3 Geocoder , I want to be able to pick up the coordinates of the route when it is dragged. I know how to pick up the coordinates giving the start and end points, but, after change the route, I don't know how to do!


texts: Automatically creating themes in QGIS 3.x

texts: I have to make around 1000 maps that represent the same area (it is reclassified satellite image - so raster). Style is the same, only date changes. I figured out how to style all layers in no time and how to print all images in Atlas. I would like to automatize also the process of making Theme (it was called preset in QGIS 2.x) in the layer panel, so I would not need to manually make all the themes I would use later in Atlas. I would like to name them after the name of the layer or raster image (it represents the date of the satellite image).


texts: QGIS Print Layout - Rotation of Grid Coordinates to match map rotation

texts:
I have a layout in QGIS where 3 map views are rotated to 35 deg. The grids automatically rotate to match the maps rotation, however the labels do not.

Does anyone know a way to make the coordinate labels match the rotation of the map?


texts: Checking if the point lies within the polygon in R using sf package

texts:
I have a point's coordinates as well as an object with the class of sfc_polygon in R. I would like to check if this point lies inside that polygon or not.
I created an object with the class of sfc_point by applying st_point() function on the point's coordinates and I used the st_intersects function to check what I mentioned above. The code is as follows:
pnts = st_sfc(st_point(c(long,lat)), crs = 3006)
iso = ABC$geom
(lst1 = st_intersects(pnts, iso))
(mat1 = st_intersects(pnts, iso, sparse = FALSE))
which points fall inside a polygon?
apply(mat1, 1, any)
lengths(lst1) > 0
I plot the polygon and the point and I am sure that the point lies inside the polygon, but the result of the code above is false!
What is the problem with the code?


texts: What is the most appropriate journal for publishing a paper about data collection with QGIS and mobile equipment?

texts:
I am a QGIS user, and a few months ago, we used QGIS on a PC Tablet to collect data on nesting birds on a cliff directly on the field. The principle is to have QGIS installed on a PC-Tablet with an SQL database, and with easy way to write data.

We are still work on it, and to share this example and cite this work, we would like to write a manuscript, but we don't know in which journal.

Does someone have advice as to which journal we could submit this kind of work?


texts: Transform a lat/lon coordinate to x/y pixel on a Google Earth screen

texts:
For a given Google Earth URL
(for example https://earth.google.com/web/@33.83405747,-117.89835255,51.77456494a,140.84823168d,35y,0h,0t,0r)
I want to find the lat/lng coordinates of the corners of the screen, so that I can calculate any pixel location on the screen given the lat/lon. The screenshot of this Google Earth page would look like this:

For example in the above link the screen size is 1280 × 720. The coordinates 33.83405747,-117.89835255 would be the exact center, so the pixel coordinates would be (640,360). Now I want to calculate these coordinates for any given lat/lon(Of course, assuming they will lie on this field of view.)
The breakdown of a Google Earth URL would be this:

33.83405747 is latitude
-117.89835255 is longitude
51.77456494a is the altitude of the center from sea level(I think it includes the building height too, if there is any at the center. It is in meters).
140.84823168d is apparently the distance of the viewer from the center(not sure what the units are and apparently it is not from the sea level).
35y is the field of view.
0h is the heading
0t is the tilt
0r is the rotation of the view

I'm assuming that I'll keep heading, tilt and rotation always 0.


texts:
On this website http://www.gisday.com/, if you scroll down a little bit, it shows a map of locations where GIS workshops are being held.

As you zoom in (or out), the number of points changes and they snap away from each other to their more appropriate location.

I looked around and thought the Integrate tool did that, but I was incorrect.

Does anyone know of a tool(s) in ArcGIS Pro that I can achieve that, or can provide insights?


texts:
I am trying to add CRS  to a shapefile that I am exporting:

shp1= gpd.read_file(i[1],encoding='iso8859_7',crs_wkt='2100')
shp1.to_file(out,encoding='iso8859_7',crs_wkt='2100')

When I import it to a GIS later, it doesn't have the CRS that says in crs_wkt.

I tried to do this:

shp1=shp1.to_crs(epsg=2100)

but says:

RuntimeError: b'no arguments in initialization list'

If check for the crs it has

shp1.crs

it says:

{}


texts: Script not taking Input Parameter from ModelBuilder?

texts:
I have a script that reads user-queried values from a file, then writes just those values to a text file. The text file is then further manipulated in ModelBuilder. I want the user to be able to pick the original file. I added a parameter called INPUTFILE to the script, with type equal to "File" and Direction equal to "Input". In the script, that variable is called.  This is the first several lines of the script:

import struct
import numpy as np
import matplotlib.pyplot as plt
fid = open(INPUTFILE,'rb')

When i run the model, it fails at the line "fid = open..." and gives this error message:
"NameError: name 'INPUTFILE' is not defined"

But it is an input to the script in the model I built. Its a variable, when you run the Model, it is one of the things you set. you set it by browsing to the file you want the script to read.

What do i have to do to get the script to accept the parameter?


texts: How to draw a route line between more than one lat long in leaflet map

texts:
now I am using Leaflet map in my application but the problem is am not able to bring the line between lat long markers. Poly line just drawing the lines but not taking along with the road..
Currently using the recent version of leaflet route machine.js file but getting the error like

leaflet-routing-machine.js:15868 Routing error:

message: "HTTP request failed: error"

url: "https://router.project-osrm.org/route/v1/driving/-74.75774,40.225438;-73.998023,40.305747;-74.223136,40.410324;-74.297394,40.287739;-74.274553,40.469638;-74.229082,40.823

target: ProgressEvent
isTrusted: true
lengthComputable: false
loaded: 0
total: 0
type: "error"

Failed to load resource: net::ERR_CERT_DATE_INVALID

Please help me anyone to achieve this.. thanks in advance..


texts: How to use CQL_FILTER for the highest value in a layer field?

texts:
I need to query/filter my wms in openlayers to only see the polygons classified with the highest values in a given field.

Would be something like:

wms_layer.mergeNewParams({'CQL_FILTER': "ND = max(ND)"});

But this is not working.

Any ideas?


texts: Statistics of Image Neighborhood in Google Earth Engine

texts:
When using the reduceNeighborhood function with a kernel, the results differ according to the zoom level. When the results are zoomed in at the maximum (you see individual pixels clearly), the kernel is computed only for the edges of a pixel. The differences are also evident when interrogating the same pixels over multiple zoom levels. See below the results from the GEE example from Statistics of Image Neighborhood.

Are there some parameters to setup in order to have the same results no matter what the zoom level is?

Zoom out image

Zoom in image


texts: Should I always choose Sentinel-2 atmospheric-corrected imagery?

texts:
Sentinel-2 imagery is available in the default L1C form, or in the atmospheric-corrected L2A form.
I haven't been able to find a L1C - versus - L2A comparison that would provide the pros and cons of each.  In the absence of such a comparison why wouldn't I always want the corrected L2A form?  Are there situations where L1C would be preferable?


texts:
I have 2 ESRI feature layers on my map; points and polylines. I can select them just fine either by rectangle or polygon. What I am missing is how to highlight (larger size and different color) the selected points and polylines for display. Do I need to create new layers with just the selected features and use a renderer or place a graphic container on the map with those features highlighted? The ultimate goal is to have the selected and highlighted features available for deletion.

Using ArcEngine and ArcSDE. C# .NET.

Thank you.


texts: gdal2tiles "ERROR 5: Illegal values for buffer size" error for a single zoom level only

texts:
I'm experiencing what seems to be a very particular issue when running gdal2tiles on a virtual raster at zoom levels 10 to 17. The function runs without issue for all zoom levels except for zoom level 16, where I receive an "ERROR 5: Illegal values for buffer size" error.

I'm using gdal2tiles version 2.1.3, and the command I'm using is nothing out of the ordinary:
gdal2tiles -s EPSG:3857 -z 16 georef\output.vrt tiles

I've run the function at various levels using a verbose command, and I believe the issue is related to the Native Extent being set to (0,0,0,611) for the first tile, which (I think) gives the tile a width of zero. The verbose output for zoom level 16 is shown below:

('Options:', <Values at 0x7f77e7661e60: {'profile': 'mercator', 'kml': False, 's_srs': 'EPSG:3857', 'copyright': '', 'resume': False, 'url': '', 'googlekey': 'INSERT_YOUR_KEY_HERE', 'bingkey': 'INSERT_YOUR_KEY_HERE', 'webviewer': 'all', 'zoom': '16', 'quiet': None, 'resampling': 'average', 'tmscompatible': None, 'title': 'output.vrt', 'srcnodata': None, 'verbose': True}>)
('Input:', 'output.vrt')
('Output:', 'tiles')
Cache: 395 MB

('Input file:', '( 91664P x 68298L - 3 bands)')
NODATA: []
('Preprocessed file:', '( 91664P x 68298L - 3 bands)')
('Bounds (output srs):', 0.0, -68298.0, 91664.0, 0.0)
('Bounds (latlong):', (-0.6135196480581518, 0.0), (0.0, 0.823431722035318))
('MinZoomLevel:', 16)
('MaxZoomLevel:', 16, '(', 2.388657133911758, ')')
Generating Base Tiles:

Tiles generated from the max zoom level:
----------------------------------------

('dataBandsCount: ', 3)
('tilebands: ', 4)
(1, '/', 16912, 'tiles/16/32767/32767.png')
('\tNative Extent (querysize', 611, '): ', (0, 0, 0, 611), (611, 0, 0, 611))
('\tReadRaster Extent: ', (0, 0, 0, 611), (1024, 0, 0, 1024))
ERROR 5: Illegal values for buffer size
ERROR 5: Illegal values for buffer size
Traceback (most recent call last):
  File "/usr/bin/gdal2tiles.py", line 2436, in <module>
    main()
  File "/usr/bin/gdal2tiles.py", line 2433, in main
    gdal2tiles.process()
  File "/usr/bin/gdal2tiles.py", line 481, in process
    self.generate_base_tiles()
  File "/usr/bin/gdal2tiles.py", line 1309, in generate_base_tiles
    dsquery.WriteRaster(wx, wy, wxsize, wysize, data, band_list=list(range(1,self.dataBandsCount+1)))
  File "/usr/lib/python2.7/dist-packages/osgeo/gdal.py", line 2000, in WriteRaster
    buf_pixel_space, buf_line_space, buf_band_space )
TypeError: not a string

As mentioned above, there are no errors experienced with other zoom levels, and verbose output for other levels shows a non-zero third argument to Native Extent and ReadRaster Extent.

For example for zoom level 15:

('dataBandsCount: ', 3)
('tilebands: ', 4)
(1, '/', 4256, 'tiles/15/16383/16383.png')
('\tNative Extent (querysize', 1223, '): ', (0, 0, 1, 1223), (1222, 0, 1, 1223))
('\tReadRaster Extent: ', (0, 0, 1, 1223), (1023, 0, 1, 1024))

I'm uncertain if this is a bug, or if I'm perhaps missing something with my configuration.


texts: Modify NAIP Images for seamless mosaic with gdal

texts:
I'm processing lots of USGS NAIP imagery which overlaps.
if I load multiple files into a QGIS project I can clearly see the edges of each tile being that some tiles have a darker or lighter contrast than neighboring tiles.

by using GlobalMapper I'm able to put a fade around all of the edges of the tiles to basically fade them together. but I'm getting away from that method to create a full script to handle my proccess much easier. is there a way to "fade" or "blend" the edges together using gdal or maybe do some type of histogram correction?


texts: How can I delete duplicate LiDAR points? I incorrectly read the same las twice and saved as a one file resulting in duplicate XY points.  How should I delete these duplicate points?


texts:
Can anyone show how to calculate the RMSE (root mean square error) between following two rasters step-by-step and discuss on the obtained results' min and max values, and how to interpret them.

 First raster (original, 2 by 2):
 1 2
 3 4

 Second raster (obtained, 2 by 2):
 2 2
 4 1


texts: Problem with snap layer when in name is space

texts:
I'm trying to make a option where user can choose layer to snap. All is okay when layers names are without space, but if there is a space in the name (like 'something layer name') snapping does not work.

How can I solve it?


texts: Error occured during an attempt to connect 10.1 SDE geodatabase on SQL Server from ArcMap 10.5

texts:
I have an SDE geodatabase created and running in version of 10.1 on SQL Server 2008. I was trying to connect to it via an SDE connection file from ArcMap 10.5.

An error was thrown out as below while it attempted to show databases on the Database Connection dialogue. I am sure the username and password are correct.

Any assistance?


texts: What do the COUNT values in the output table of a viewshed RASTER (from multiple locations) actually indicate?

texts:
In ArcGIS (10.1), I ran a VIEWSHED analysis using as input observer locations a POINT shapefile with 6 locations. When I open the output RASTER's attribute table, I seen that it contains a number of rows equal to the input locations, and a COUNT field.

Now, I seem having problem grasping the meaning of the latter field, and the documentation I found seems not clarify my doubt. My best guess is that those values represent the number of cells that are visible from each locations. Am I right, or am I missing something?


texts: Overlap feature's style with layer style

texts:
I have a Vector layer, where each feature is created using a image with different rotations.

var point = this.__getFeaturePoint();

featureStyle = {
    externalGraphic: this.__getImageURL(),
    graphicHeight: 50,
    graphicWidth: 50,
    graphicXOffset: -25,
    graphicYOffset: -50,
    rotation: this.getAzimuth()
};

feature = new OpenLayers.Feature.Vector(point, {}, featureStyle);

I would like to have a 'select' style so the image URL changes, but the rotation remains.

style = new OpenLayers.Style({
    graphicZIndex: 1
});
selectStyle = new OpenLayers.Style({
    graphicHeight: 50,
    graphicWidth: 50,
    graphicXOffset: -25,
    graphicYOffset: -50,
    fillOpacity: 1,
    graphicZIndex: 2
});

this.__layer = new OpenLayers.Layer.Vector("Antennas", {
    renderers: ["SVG"],
    rendererOptions: { zIndexing: true },
    styleMap: new OpenLayers.StyleMap({
        "default": style,
        "select": selectStyle,
        "highlight": highlightStyle
    })
});

It seems that only the feature's style is considered. How can I overlap both styles? Or maybe implement such a styleMap for each feature?


texts:
I have a set of files made for use with ArcGIS Desktop that I would like to use with ArcGIS Pro. One of these is an .AVL symbology file, which on first light is not recognized by ArcGIS Pro.
Is there a way to use these files in ArcGIS Pro for setting symbology?


texts: Bug using Python 3.7 (from QGIS) and trying to install package with pip

texts:
I use QGIS 3.6 and Python 3.7 and I try to install a package with pip in command line tool under Windows 7 but I have the following message :

Fatal Python error: initfsencoding: unable to load the file system codec
ModuleNotFoundError: No module named 'encodings'

Current thread 0x0000188c (most recent call first):

My path and PYTHONPATH are well configured with link to QGIS' Python folders. The problem appeared recently with use of Python 3.7.


texts: Problem with numeric variable in i.topo.corr

texts:
I have a script to do topographic correction on landsat images.

Following the manual, I first have to create the illumination model. I do it with the following line:

i.topo.corr -i --overwrite output=temp basemap=DTM zenith=$zenith azimuth=$azim

When I launch it through the script the output is:

ERROR: Sorry <127.95340672> is not a valid option
ERROR: Missing value for parameter <azimuth>

What makes me crazy is that if I launch it interactively with the same numbers it works!

Any idea?


texts: Applying symbology from layer using ArcGIS ModelBuilder?

texts: I am trying to automatically apply symbology to my polygons based on a .lyr file by using "apply symbology from layer" on a .shp file. It seems to work if the generated .shp file from the module is manually selected (blue bubble) but not the very same .shp file produced from the module within modelbuilder. See image below for clarification.


texts:
Many recent maps I have seen have these neat two-way color gradients that can be used to color correlations between continuous estimates.

Here is an example is taken from a recent publication.

Christian Levers, Pieter J. Verkerk, Daniel Müller, Peter H. Verburg, Van Butsic, Pedro J. Leitão, Marcus Lindner, Tobias Kuemmerle - Drivers of forest harvesting intensity patterns in Europe. Elsevier. 1 March 2014 (access: http://dx.doi.org/10.1016/j.foreco.2013.12.030).

I believe those gradients have been generated with ArcGIS 10.x and I was wondering how I could reproduce such color gradients with freely available GIS tools (or at least those available on Linux)?

Happy to hear suggestions or code snippets.


texts:
I am looking for a solution for this with no success using:

result = arcpy.mapping.ConvertWebMapToMapDocument(Web_Map_as_JSON, templateMxd)

AND

arcpy.mapping.ExportToPDF(mxd, Output_File, georef_info=True)

All my services are in WGS84 web mercator. My only successful export was in arcmap with all data in WGS84 and dataframe in UTM, weird.


texts: Trying to combine nearby points into a polygon

texts:
We have collected millions of cell phone strength readings. I am trying to group points of the same strenght that are within a distance of each other in to a polygon. I asked on the
postgis mailing list but didn't get a response.
To make the process simpler, I already extracted all the point of the same strength to one table. So I have a table of 1.5 million rows. It has a gist index on the geom field and btree indexes on the uniqueid of each row and on the 'clusterid' which is supposed to be the same for all points that are in the same 'cluster' (ie point a which is within distance of point c which is within distance of point z would all be the same cluster).
How would you accomplish this task???
I did wonder whether it might be worth trying to load the geometries in to memory and use JTS to process them? Although I doubt that would be any better.
Thanks - Bryan.


texts: What is the best method of visualizing a network on a map?

texts:
I have several hundred geo-referenced data points, and the relationships from that point, to other points. I'm trying to figure out the best way of visualizing this on an interactive map (possibly using google maps).

One idea I had was that when a user clicks on a point, it then displays all the links from that point, to the related points.

I would appreciate any suggestions or examples of how best to do this.

I have experience using ArcGIS, QGIS, Python and a small amount of javascript.


texts: In ArcGIS raster, how to calculate the average value of each region group

texts:
I am doing a ArcGIS raster project. Now I am at the step where I already created a scoring system for my targeted area and I assigned the score to each of the cells in side the targeted area, as shown in the picture below

Now if I want to find the average value of the cells of each region group.

What tool shall I use?


texts: How to get fieldnames in a value table depending on the input value

texts:
I would like to add a parameter as value table with one of the fields being a value list. The value table consist of 3 columns, namely "table","field","description", whereas "field" is a value list. The values of the value list shall depend on the input "table, thus change with each new entry. E.g. adding a table xyz to the value list should lead to a list of fields which can than be selected.

My code for now:

        def getParameterInfo(self):
            in_tables = arcpy.Parameter(
                displayName = "test",
                name = "test",
                datatype = "GPValueTable",
                parameterType = "Optional",
                direction = "Input")
            in_tables.columns = [["DETable", "table"],["GPString","field"],["GPString","description"]]
            in_tables.filters[1].type = 'ValueList'
            in_tables.filters[1].list = ["<marker spalte>"] #dummy list which is updated in updateParamters

        ...
        def updateParameters(self, parameters):
            if parameters[0].value and parameters[0].altered:
                fields_list = [field.name for field in arcpy.Describe(parameters[0].value[0][0]).fields]
                parameters[0].filters[1].list = fields_list #update of dummy list

This actually works but I always get the field list of the very first entry. I thought about using len in order to always get the last entry (which is a bit more useable but still not satisfying). What I would prefer is the identification of the selected row in order to dynamically update the value list with the related field names.

Is there a way to do that?

ArcGIS 10.3.1

Update:
The problem with my parameters[0].value[0][0] solution is, that it is not working if the any of the other tables has other fields than the first one. the problem when working with len() is that an error occurs if one of the former selected fields does not exist in the last field list. Thus I changed my updateParameters method.

        if parameters[0].value and parameters[0].altered:
            all_tables_list = parameters[0].value
            number_of_tables = len(all_tables_list)
            i = 0
            fields_list_tmp = []
            for in_table in all_tables_list:
                fields_list_tmp.append([field.name for field in arcpy.Describe(parameters[0].value[i][0]).fields])
                i += 1
            fields_list = list(set([item for sublist in fields_list_tmp for item in sublist]))
            parameters[0].filters[1].list = fields_list

which basically creates a list with all existing fields in all selected tables.


texts: Creating a Time Aware Population Layer with Point Features

texts:
I am interested in using ArcGIS Pro's time aware capabilities to create a time aware map of the populated places for the state of Wyoming from the 1870 Census to the 2010 Census.  In the end, I would like to be able to adjust the time slider to display the points using graduated symbols based on that population value which use the same method for every year (i.e. symbol size classifications remain constant over the entire range of time).  Then I will publish to ArcGIS Online.

Currently, I have all of the places in a point feature layer in a geodatabase with fields containing population counts for every census during those years.  It is my understanding that I need to add a time field and then replicate each feature for each census year which I have data for that specific feature and then add the population for that specific year.

Is this the correct approach and, if so, is there a simple way to produce this table?  While there are not many cities in Wyoming there are still 15 decennial census over that time period so I would like to speed up this process by not manually moving these data.


texts: How to use Mapbox geocoder in OpenLayers application?

texts: I want to use Mapbox geocoding services in OpenLayers 5 project. I'm using this method https://github.com/jonataswalker/ol-geocoder  but it doesn't support Mapbox API


texts: Simplified polygons have overlapping lines

texts: I am working with QGIS. I am creating a contour world map, I started my project with a raster image of the world and then extracted contour lines much too complex for my purposes. I simplified the polygons using a line generalization plugin. However I noticed that after simplifying some of my polygons are lapping into each other, how can I fix this.


texts: Removing events map.on("click"...) from mapbox gl?

texts:
I remove the layers with map.removeLayer("clusters"); but got this error

evented.js:111 The layer 'clusters' does not exist in the map's style
and cannot be queried for features.

as I have map.on('click', 'clusters', function (e) { ... event. So how to remove this?


texts: Query for tables that are not registered with option to move edits to base

texts:
I have an Oracle enterprise geodatabase. In theory, all 500+ tables in the database should be registered as versioned, with the option to move edits to base.

However, it has come to my attention that a few of the tables were not registered with the option to move edits to base.

I would like to find all tables that were not registered with base (and then manually re-register them correctly).

Is it possible to write an SQL query for this purpose?


texts: Editing 3D object in ArcScene and exporting as a VRML file

texts:
This might be easy, but I'm a newbie in using ArcScene and I can't figure out how to do it. I have a 3D terrain model - srtm dem (*.tiff) file visualized in ArcScene. I need to add to this model some 3D objects like a sphere, cube etc. and export it all as one *.wrl file.

I can't find the right way. It works when exporting a terrain model alone. Then I tried to add some point features (*.shp) and visualized it as a 3D simple marker - sphere. But when I try to export it, it always crashes down. Probably because of the incompatibility of tiff and shapefile?

So I made a sphere in SketchUp and exported it as a vrml file and tried to add it to the terrain model in ArcScene. But I dont know how to put it right on the model surface and how to change the size and colour, if possible. And the exportation fails as well. So I can't export more than one layer to the *.wrl file?

How to do it then?


texts: Connecting points to polygon edges using angular bearing

texts:
I am trying to map some caves located along a coastline (using QGIS 2.18.3 with GRASS 7.2.0). The survey took place on a boat, so the GPS points were taken at sea.

I have information on the angular bearing (ranging from 0-360 degrees, where 0 is North) of the caves from the boat and the lat-long position (in decimal degrees) of the boat. I plotted this information as shown below - [purple is land, yellow is position of the boat, and "rumo" is the angular bearing]

I need to find the corresponding point on the coastline based on the angular bearing and distance to shore (edge of polygon). I have tried looking at various options to do this, but as far as I can see there are ways to find the nearest distance to polygon, but not along any angular direction. The Azimuth and Distance plugin also seemed promising, but from what I understand I have to manually enter each angular bearing. Since I have >1000 such points, it is not very feasible! Seems like there is a tool in ArcMap to do this - Bearing Distance to Line. But there doesn't seem to be a direct equivalent in QGIS.

I am far from an expert in QGIS, but could be pushed into writing some code.

A sample dataset and related polygon can be found here


texts: Setting the correct zoom level from Sector bounds for flat map using WorldWind SDK

texts:
This is not coming out correct for me. Setting the zoom level based on some other code I got on here using the EarthFlat model, which does not set the correct zoom for the placemarks within a sector bounds.

private void invoke() {
    SwingUtilities.invokeLater(new Runnable() {

        public void run() {

            ArrayList<Position> pos = new ArrayList<Position>();
            pos.add(Position.fromDegrees(32, -117));
            pos.add(Position.fromDegrees(48, -72));
            pos.add(Position.fromDegrees(40, -70));

            final Sector bounds = Sector.boundingSector(pos);

            setWwPanel(new WWPanel());
            final WWPanel wwPanel = getWwPanel();
            wwPanel.setPreferredSize(new Dimension((int) getSize().getWidth(), (int) getSize().getHeight()));
            wwPanel.getWwd().addGLEventListener(new GLEventListener() {

                public void display(GLAutoDrawable e) {
                    //System.out.println("Map displayed");
                }

                public void dispose(GLAutoDrawable e) {
                    System.out.println("Map disposed");
                }

                public void init(GLAutoDrawable e) {
                    System.out.println("Map GL initialized");
                }

                public void reshape(GLAutoDrawable e, int arg1,
                        int arg2, int arg3, int arg4) {
                    System.out.println("Map reshaped");
                }

            });

            LayerList layerList = new LayerList();
            BMNGWMSLayer bmng = new BMNGWMSLayer();
            layerList.add(bmng);
            ViewControlsLayer vcl = new ViewControlsLayer();
            vcl.setShowPitchControls(false);
            vcl.setShowVeControls(false);
            vcl.setShowHeadingControls(false);
            layerList.add(vcl);
            Model flatModel = new BasicModel(new EarthFlat(), layerList);
            ((EarthFlat) flatModel.getGlobe()).setProjection(FlatGlobe.PROJECTION_MERCATOR);
            wwPanel.getWwd().setModel(flatModel);
            wwPanel.getWwd().addSelectListener(new ViewControlsSelectListener(wwPanel.getWwd(), vcl));

            RenderableLayer rl = new RenderableLayer();
            for(Position curPos : pos) {
                PointPlacemark pm = new PointPlacemark(curPos);
                pm.setAltitudeMode(WorldWind.RELATIVE_TO_GROUND);
                pm.setLabelText("Test");
                PointPlacemarkAttributes pma = new PointPlacemarkAttributes();

                pma.setImageColor(Color.GREEN);
                pma.setLabelFont(Font.decode("Verdana-Bold-11"));
                pma.setLabelMaterial(Material.GREEN);
                pm.setAttributes(pma);
                pm.setValue(AVKey.DISPLAY_NAME, "Text displayed on mouse over");

                rl.addRenderable(pm);
            }

            layerList.add(rl);

            BasicOrbitView fov = (BasicOrbitView) wwPanel.getWwd().getView();
            fov.setCenterPosition(Position.fromDegrees(bounds.getCentroid().getLatitude().degrees, bounds.getCentroid().getLongitude().degrees));

// this is what i am using to set zoom ??????
            double radius = wwPanel.getWwd().getModel().getExtent().getRadius();
            double viewRange = Math.sqrt(Math.pow(bounds.getMaxLatitude().degrees-bounds.getMinLatitude().degrees, 2)+Math.pow(bounds.getMaxLongitude().degrees-bounds.getMinLongitude().degrees, 2));
            double altitude = radius * Math.sin(Math.toRadians(viewRange * 0.5));

            fov.setZoom(altitude);

            add(getWwPanel(), BorderLayout.CENTER);
        }

    });
}

It sets the zoom fine, but the placemarks are off the view bounds.

Expected result:

Undesired/actual result:


texts:
I have two 3.6 MB GeoTIFF compressed with LZW. When I merge them with:

gdal_merge.py *.tif -o ~/out.tif  -co COMPRESS=LZW -co BIGTIFF=YES

The files is 10x larger at 31 MB.

Is there some reason the filesize would increase when merging? They are a continuous region (all touching each other), in EPSG 3413 projection.

This question was asked here (How to merge raster scenes while maintaining small file size?) but not answered. That OP selected the "use compression" answer, but my compressed file is larger. I'm posting gdalinfo for one input file (first) and the final file (below).

NOTE Oddly, when I do this for >2 files, or the whole folder which is 800 MB, it is 8.4 GB uncompress, and LARGER, or 11 GB with LZW compression.

Driver: GTiff/GeoTIFF
Files: gimpdem0_0_v1.1.tif
Size is 8310, 15000
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Polar_Stereographic"],
    PARAMETER["latitude_of_origin",70],
    PARAMETER["central_meridian",-45],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",1],
    PARAMETER["false_northing",1],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-640000.000000000000000,-2905550.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -640000.000,-2905550.000) ( 57d25'19.50"W, 63d 1'28.29"N)
Lower Left  ( -640000.000,-3355550.000) ( 55d47'53.80"W, 59d11'57.96"N)
Upper Right ( -390700.000,-2905550.000) ( 52d39'30.45"W, 63d24'19.24"N)
Lower Right ( -390700.000,-3355550.000) ( 51d38'28.63"W, 59d31'30.22"N)
Center      ( -515350.000,-3130550.000) ( 54d20'53.46"W, 61d18'11.32"N)
Band 1 Block=8310x1 Type=Int16, ColorInterp=Gray

And the merged file:

Driver: GTiff/GeoTIFF
Files: /Users/kdm/Desktop/20171011-204436.tif
Size is 8310, 30000
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Polar_Stereographic"],
    PARAMETER["latitude_of_origin",70],
    PARAMETER["central_meridian",-45],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",1],
    PARAMETER["false_northing",1],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-640000.000000000000000,-2455550.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -640000.000,-2455550.000) ( 59d36'29.72"W, 66d52'51.74"N)
Lower Left  ( -640000.000,-3355550.000) ( 55d47'53.80"W, 59d11'57.96"N)
Upper Right ( -390700.000,-2455550.000) ( 54d 2'25.84"W, 67d20' 7.02"N)
Lower Right ( -390700.000,-3355550.000) ( 51d38'28.63"W, 59d31'30.22"N)
Center      ( -515350.000,-2905550.000) ( 55d 3'28.16"W, 63d14'13.92"N)
Band 1 Block=8310x1 Type=Int16, ColorInterp=Gray


texts:
Recently installed the GDAL software on my Windows 10 computer for the first time. Process indicated the install was successful.
When I did the command gdalinfo --version I get gdal 3.0.0 with the release date

I am not sure what the syntax of the pip install command should be
I see a suggestion of "pygdal>=2.1.2<2.1.3" for older versions.
Tried pip install "pygal>=3.0.0" and get a message
"Could not find a version that satisfices the requirement 3.0.0.
Second message states no matching distribution found. Last number of versions shown is 2.4.0


texts: How do I add a .vrt file to a QGIS project?

texts: I've created a vrt file through a text editor. Trying to add it to a project has me baffled.


texts:
Is it the case that when I use Arc GIS to compute nearest distance from a point to a poly-line/polygon, I need to project them both in an equidistant projection? Note: I am using the near tool.

What happens, for example, if I do not project them at all and so the everything is just in the Geographic Coordinate System (WGS1984 in my case)?


texts: Determining what side of line point is on using PostGIS?

texts:
I want this in order to generate  house numbers, odd or even depending on wether a feature is on one side of a street line or not .  Thanks

In the end i think the  most portable way is to use this( BTW line is actually 'linestring' not segment):

ST_LineCrossingDirection( linefromAtoB,streetline) <0  where A is the point , B is  ST_ShortestPoint(thepoint,streetline); if <0 then linefromAtoB enters through the left of The streetline; if >0 the it enters from right side of streetline). You have to determine what street line is closest to each point first.

I made a mistake:
C is the symmetric of A in relation to B=  ST_ShortestPoint(thepoint,streetline).
C's coordinates are XC=2*(XB-XA) -XA and for y   YC=2*(YB-YA) -YA.
ST_LineCrossingDirection( linefromAtoC,streetline) <0  where A is the point , ; if <0 then linefromAtoC enters through the left of The streetline; if >0 the it enters from right side of streetline). You have to determine what street line is closest to each point first.


texts: How to get multiple SLD files of the same shapefile?

texts: I need to get multiple SLD files for the same shapefile. Geometry is polygons, and they need to be colored by one attribute column.  What is the best solution? I'm working with ArcMap2SLD or ArcMap to SLD Converter.


texts: Adding shapefile to multiple MXDs using ArcPy?

texts:
My ArcMap verison is 10.6. The purpose is to add a wetland shapefile to multiple MXDs. It doesn't work and the error said the addLayer is not defined. My script is below.

Does anyone have an idea?

import arcpy,os
arcpy.env.workspace = ws = r"C:\Users\Rachel\Desktop\Example\Wetland.shp"
wetland = r"C:\Users\Rachel\Desktop"

mxd_list = arcpy.ListFiles("*.mxd")

for mxd in mxd_list:

current_mxd = arcpy.mapping.MapDocument(os.path.join(ws, mxd))

dFrame = arcpy.mapping.ListDataFrames(current_mxd, "Main Map")[0]
addLayer = arcpy.mapping.Layer(wetland)
arcpy.mapping.AddLayer(dFrame, addLayer)

current_mxd.save()
print("done")

del mxd_list


texts: Generating Interactive Tiles using Nodejs, Mapnik, and PostGIS for Leaflet

texts:
I am looking to do something like this but where I have all the data on the same server.

I've been following through the set of tutorials @ sparkgeo.com/labs/big/ and I'm having some difficulty. Ultimately I'd like to get it working with Leaflet like the tutorial @ sparkgeo.com/labs/bigleaflet/ but I can't get the version on the prior tutorial to work either. So the final version will be generating R Shiny interactive maps and layers from tiles accessed through Leaflet generated by Postgis by means of node.js using mapnik and mapbox.

The Wax libraries used in that tutorial have been removed so I installed this. I pointed the tiles in the HTML file to http://localhost:8000/{z}/{x}/{y}.png and pointed to my location of wax.g.min.js. No census tracts are showing up with I run the cencon.js code. Similarly with the Leaflet version of the HTML file, I pointed to the localhost for the layers and they are not displaying. They are successfully loaded to the PostGIS database and changing the database login info.

OK, so I'm new to all of this and I'm wondering if I'm missing something obvious. I did not alter the code other than what I previously mentioned. Is there something else I needed to do to generate the tiles @ http://localhost:8000/{z}/{x}/{y}.png? Also, now that Wax is superseded, do you have any links for more suitable tutorials using mapbox?


texts: I'm new to least cost analysis. I have a DEM-layer of a region and I want to know what's the least cost path from a to b. Do I have to use the DEM layer to calculate that? I also have a hillshade layer and a slope layer of that region. I think the DEM-layer only has the heigth of the landscape, but doesn't slope affect the costs of a specific route as well? Do I have to combine these layers? I have QGIS 3.8 and 3.6


texts: Where to obtain free 15m ASTER image for New Mexico?

texts:
I am looking for a free Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) image that will let me to perform NDVI for the pilot crop health/soil moisture that can go as far back ten years or so.

It needs to have:

15 m (meters)
temporal
Ten years back or so
Perform NDVI

This will be in a small area somewhere near Albuquerque, NM.


texts: The Add Feature Button is not present on toolbar in QGIS 3.6

texts: This may sound crazy, but I do not have an Add Feature button on my toolbar in QGIS 3.6.  Am I missing something?


texts: Why do features in QGIS drawn on same datum as Google imagery from OpenLayers move?

texts:
I have just started using QGIS. I've loaded the "openlayers" plugin and added a google hybrid layer, which is based off the WGS84 datum to my understanding. Now I create a shapefile and add a new line vector layer based on the same datum. I draw a line and then zoom out. Now, the line has moved with regards to the google map (i'm drawing a road and now the road no longer lines up with the original road I'm using for a guide).

My end goal is to make a georeferenced shapefile with some roads. Is my methodology flawed or is their a problem with my plugin (assuming the former of course!)?


texts: ArcGIS 10.5 Spatial Analyst Extract-by-Mask work around?

texts: I have 11 different shapefiles all using NAD27 Lambert Conformal Conic with which I use ArcGIS 10.5 Extract By Mask to extract a chunk of a forest map of the same projection parameters. But each time I extract with these shapefiles, I end up with an extraneous portion of a tif to the lower right of the original shape/tif file that did not exist before hand, that only comes about as a result of the Extract By Mask. It is real. I can see it, I can see the values in it. I can probably fill in the bogus values with 0's. But how to tell ESRI this is happening? How to eliminate this? What is the workaround, so this does not reoccur? The extraneous portions are different.


texts: Count the number of points in a (multi)polygon in shapely?

texts:
I'm using Python, Shapely and Fiona to process some shapefiles.

Is there an easy way to count the number of points in a polygon/multipolygon in Shapely and Fiona?


texts: Using one to many join to *.xlsx in ArcGIS Desktop?

texts:
I'm working with ArcMap 10 and have some trouble joining a shapefile and a .xlsx file.

The .xlsx contains sensor data (CO2, NO2 etc.) without location information. The shapefile contains all the sensor locations. Sounds easy but here is the problem:

The .xlsx file contains data from multiple points in time (every 5 seconds of a single day). It looks like this in ArcMap:

And the Shapefile attribute table looks very simple:

I Want to join the two based on ID, which is the matching sensor number. A basic join obviously won't do the trick. I tried to mess around with creating a Query Table but could not make it work. How do I do this?


texts: Create Postgis GEOM from OSM bounding box

texts:
I'm trying to use the data that is retuned from the nomatim API to perform a simple geospatial search. The users search for a place, i.e Leeds, and this returns simple info, including the following bounding box:

"boundingbox": [
      "53.6983747",
      "53.9458558",
      "-1.8004213",
      "-1.2903452"
 ],

Using this, I should be able to create a bounding box using the Postgis function:

geometry ST_MakeEnvelope(double precision xmin, double precision ymin, double precision xmax, double precision ymax, integer srid=unknown);

However, when I do this, I am unable to get correct geometries that I can then use, and plug into another function wherein I check if a point is within that geometry.

As far as I can tell, the SRID for the OSM data is 3857. When I perform the following:

ST_MakeEnvelope(53.6983747,-1.8004213,53.9458558,-1.2903452, 3857);

And then perform a search on this, I have to add 500000 metres to the ST_DWithin to be able to get a point which is right in the centre.

Does anyone know If there is a better way to approach this, or spot any mistakes in this? I have tried a couple different approaches and I cannot find any more accurate than this one.

For reference, this is the search I've been performing:

Select * from geometries
where ST_DWithin(ST_Transform(geom,27700),ST_Transform(ST_SetSRID(ST_MakePoint(1.546033, 53.798075),3857),27700), 70)

at 75 I get the Newcastle box, and at 80 I get the Leeds box, and that point is in Leeds. I am very confused.


texts:
I use openlayers 2.12.

I created a form where the user can insert a bounding box, to insert to the database a new borough.

To check the EPSG of the borough I tried the following

//ba bb bc bd came from the form
var thisIsNew=new OpenLayers.Bounds(ba, bb,bc, bd).getProjection();
alert (thisIsNew);

and also

var thisIsNew=new OpenLayers.Bounds(ba, bb,bc, bd);
var newest = thisIsNew.getProjection();
alert(newest);

I never get an alert and on the console I see that "the object does not support that method", for both of the above ways.

What am I missing?

Thanks in advance


texts: How to create bounding boxes that are not orientated to the north?

texts:
I want to create bounding boxes around polygons, the bounding boxes will then be used in a ST_Hausdorffdistance to determain the squarness of every polygon. My first choice is using the ST_Envelope function but the problem with this is that all the envelopes are orientated to the north.
The result is the following image.

The following image is how I would like it to look.

This is created with the Minimum bounding geometry tool in ArcGIS. Is there any possibilty to do this in PostGIS?
I tried doing it with the Box2D function in PostGIS, but I can't seem to find the way to save the geometry to the geom column?


texts:
I am composing an Arcpy tool to help users cache their current mxd's webmaps (our own services). I want to know a way to call the ArcGIS Web Map's scales with ArcPy. Once I have the web maps scales, I can apply policies and pass the desired scales to arcpy.gp.ExportMapServerCache for extraction.

Notably, the geoprocessing tool "Export Map Server Cache" calls the maps scales into the tools dialogue, a user then ticks the scales they want the cache. I figure if the tool can do it then some code can do it too...

My end goal is a tool that gets a web maps scales, applies an appropriate extraction extent polygon for the scale range, and extracts appropriate scales into a cache.

##### For Many Web Cache
arcpy.CreateFolder_management(UserProjectPathCache,lyr.name)
arcpy.AddMessage("    New Cache Folder:  "+ UserProjectPathCache +"\\"+lyr.name)
ServiceCache = (UserProjectPathCache +"\\"+lyr.name)

##### Passing To Publisher Connection String:
ExportService = UserProjectPath +"\\Temp_Pub_Connection\\"+lyr.name+".MapServer"
arcpy.AddMessage("    Temp Publisher Con: "+ExportService)

#TAKEN FROM THE GLOBAL IMAGERY SERVICE
S1 = 591657527.591555
S2 = 295828763.795777
S3 = 147914381.897889
S4 = 73957190.948944
S5 = 36978595.474472
S6 = 18489297.737236
S7 = 9244648.868618
S8 = 4622324.434309
S9 = 2311162.217155
S10 = 1155581.108577
S11 = 577790.554289
S12 = 288895.277144
S13 = 144447.638572
S14 = 72223.819286
S15 = 36111.909643
S16 = 18055.954822
S17 = 9027.977411
S18 = 4513.988705
S19 = 2256.994353
S20 = 1128.497176

DefaultScales = "{0};{1};{2};{3};{4};{5};{6};{7};{8};{9};{10};{11};{12};{13};{14};{15};{16};{17};{18};{19}".format(S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12,S13,S14,S15,S16,S17,S18,S19,S20)

target_cache_path = ServiceCache
export_cache_type = "TILE_PACKAGE"  #"CACHE_DATASET" # "TILE_PACKAGE"
copy_data_from_server = "COPY_DATA" # Because the server might not have write access
storage_format_type ="COMPACT" # Easier To Move
num_of_caching_service_instances = 1 # I pulled this number from the air
area_of_interest = CacheAOI
export_extent = ""
overwrite ="MERGE"

arcpy.AddMessage("    --------------")
arcpy.AddMessage("        I: CACHING WEB MAP TO: "+ServiceCache)

################## EXPORT MAP CACHE HERE
try:
        arcpy.AddMessage("        I: CACHING WEB MAP: PLEASE WAIT...")
        arcpy.gp.ExportMapServerCache(ExportService, target_cache_path, export_cache_type, copy_data_from_server, storage_format_type, DefaultScales, num_of_caching_service_instances, CacheAOI, export_extent, "MERGE")
        arcpy.AddMessage("        I: CACHING WEB MAP: Completed")
        del ExportService
except Exception as e:
        arcpy.AddError("       "+e.message)


texts: Pink Tiled Map with Openlayers and Geoserver

texts:
I'm having a problem with the pink tiles when trying to display a map.
I'm using Geoserver 2.6. In the main html page, I reference to Openlayers with:

<script src="openlayers/lib/OpenLayers.js"></script>

This is how I define my layer:

new OpenLayers.Layer.WMS("dim_province",
                    "http://localhost:8080/geoserver/wms", {
                        layers: "dw1:dim_province",
                        transparent: true,
                        format: "image/gif"
                    }, {
                        isBaseLayer: false,
                        buffer: 0,
                        visibility: false
                    }
                ),

The layer is displayed with no problem on Geoserver>preview as OpenLayers, however It shows a pink tile image when I'm running the html page.


texts: Cartodb not loading layers on IE8-10 on longer requests

texts:
Check the demo here: http://unirazz.com/cartodb/cartodb1.html
After map is loaded, click on GO button to load cartodb layers. It works in latest browsers excepts IE8-10. If I reduce the length of request then it loads fine otherwise it failed. As per documentation I found that cartodb automatically sends POST request if the length of GET is large. In case of IE8-10 its always using the GET request no matter whats the size of request is.

Also I Found that as on IE8-10 there is limited support of CORS that could be also the cause that layers are not loading correctly.

How can fix it?

The code is simple, just loading layers like this

layer.getSubLayer(0).setSQL("SELECT * FROM cartodb2 where zcta5ce10 in ('" + blueZips.join("','") + "')");
        layer.getSubLayer(1).setSQL("SELECT * FROM cartodb2 where zcta5ce10 in ('" + redZips.join("','") + "')");
        layer.getSubLayer(2).setSQL("SELECT * FROM cartodb2 where zcta5ce10 in ('" + greenZips.join("','") + "')");


texts: Postgis ST_NearestValue() returns nothing

texts:
My query in  (Dutch SRID, in meter):
 SELECT ST_NearestValue(ahn.rast, 1, ST_SetSRID(ST_Point(135742.6965165309,453864.2795337387), 28992)
   FROM ahn2.ahn
  WHERE ST_Intersects(ahn.rast, ST_SetSRID(ST_Point(135742.6965165309,453864.2795337387), 28992)::geometry, 1))

returns nothing in pgAdmin4:

I am sure the point is located in the middle of a raster. Slightly changing the coordinates results in:

135742.6965165309,  453864.2795337387 (original)
135740.6965165309,  453864.2795337387 (returns value)
135744.6965165309,  453864.2795337387 (returns value)
135742.6965165309,  453865.2795337387 (returns value)
135742.6965165309,  453863.2795337387 (returns value)
135741.6965165309,  453864.2795337387 (also nothing)

Why is my query returning nothing and how can it be fixed?
Background:
I'm working on a script that needs to calculate whether or not a certain locations has a shadow. To do that a line of sight (LOS) between the ground and the sun is calculated. This line is in fact a list of coordinates (x, y, and z). To determine whether or not there is shadow the height of the terrain needs to be determined for each x and y coordinate. If it is below the z (from the line) than the LOS is not blocked. The above query serves to extract the terrain height.
Version info
POSTGIS="2.4.4 r16526" PGSQL="100" GEOS="3.6.2-CAPI-1.10.2 4d2925d" PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 2.2.4, released 2018/03/19" LIBXML="2.7.8" LIBJSON="0.12" LIBPROTOBUF="1.2.1" RASTER


texts: Examining each connected pixel to select its neighbors to separate clusters using ArcGIS for Desktop?

texts: I am trying to develop a methodology for separating pixels from one another to determine whether they belong to a city cluster or not. Starting at the center pixel of the city I want to select all pixels that fall under a certain pixel neighborhood. I want to examine a few different ones 3x3, 4x4, etc.. The selection would expand outward, examining each connected pixel to select its neighbors, and so on, until there are no neighbors for any of the pixels. I have tried focal statistics in arcgis, but I cannot find out if there is a way to use the tool on just one pixel. The picture below shows the issue. I cannot tell if the pixels circled in black are actually connected to the center pixel. This is the output of focal statistics showing number of neighbors in 3x3 rectangular neighborhood. I would prefer a python solution if possible.


texts: Hiding Map Elements Outside of Neatline in ArcMap

texts:
I have a map layout in ArcMap 10.1 that incorporates a neatline that is offset from the edges of the layout by about an inch on the inside. Certain elements of my map extend beyond the neatline and layout view.

Does anyone know of a way to clip or hide the elements that extend between the neatline and edge of the layout when exporting the map?

Thanks!


texts: Looking for a tool that will place the cursor to a specific x,y coordinate

texts:
I am a new user to Quantum GIS. I am looking for a tool/function that will allow a user to type in coordinates, and place the cursor on that exact position. This is for navigation in the field...

Here is a bit of context...

I need to provide a laptop with an interactive map to a colleague. He will hike around with a gps and the laptop. When he stops, he will check his position on the gps, and (hopefully) be able to type in the coords in QGIS on the laptop to see were he is on the map.

I know how to do this with ArcGIS (its the "go to x,y" tool).

Hopefully there is something similar in QGIS.

Thanks a million!


texts: Using onchange and onclick together in Google Earth Engine

texts:
In this code

var select = ui.Select({

  items: ['1','2'],

  onChange: function() {
    Map.onClick(test)

  }
});

Map.add(select)

var test = function() {

  print ("123")

}

When I change the dropdown second time after I clicked on the map the test function runs 2 times , for the third time I change dropdown it runs 3 times after click. I want for every change in the ui.Select after the click the function will run only one time.

What should I do?


texts: QGIS Server - Docker - Multiple Projects

texts:
I am currently trying to improve my companies server architecture with the help of Docker. I have successfully managed to run a docker container of this image: https://hub.docker.com/r/camptocamp/qgis-server. However the container only allows to serve one project by default.

Does anybody know how to serve multliple projects with a qgis server docker container?

I have tried to run a docker container of this image: https://github.com/kartoza/docker-qgis-server. However it didn't work. I guess the best approach is probably to change some configuration files in the docker container of the first mentioned container but I have no idea which files.


texts: Using Google Earth to figure out the coordinate reference system from only lon, lat information

texts:
Lately, I've been trying to use lat/long information to plot on QGIS without any good results because I lack the correct EPSG code for the coordinate reference system. This is not the case with Google Earth, which always plots in the correct places using lat/long information.

I concluded that Google has an automatic way to discern them an select the correct coordinate reference system, and other posts have also suggested it. Is there a way to know which one are they using for my data?

45.801417207831314,9.092656559575259
45.801417207831314,9.092656559575259
45.801417207831314,9.092656559575259
45.801417207831314,9.092656559575259
45.801421972667725,9.092675509431874
45.80142711235822,9.092700934521718
45.8014346957649,9.09272071353453
45.80144571424553,9.092746763378287

Note: my data is from Como, Italy.


texts: Raster cell size changing with new projection

texts:
I have a high-resolution World View 2 satellite image with a cell size of 0.5 m. However the image is in a WGS84 coordinate system, so to get meter distance values instead of degrees I have converted it to a WGS84 UTM projection.
The raster in UTM projection has a cell size of 0.35.
How can I convert the image into a UTM projection and maintain the 0.5 m cell size?
I am using ArcGIS Pro and converted the coordinate system from WGS84 to UTM using "Project Raster" tool.


texts: Calculating distances and relationships between point land line vector layers using QGIS?

texts:
I have two vector layers. One lines and the other points.

I am trying to find  how far dots are apart from lanes, on average?

I am thinking about setting a buffer around the lines, and finding how many points would be in that buffer, but may be there are other simpler functions in QGIS


texts: imposm3 on PostGIS to mapnik conversions

texts:
I am working with osm2psql and imposm3 to load a postgres db with OSM data.

The performance of imposm3 is seemingly far in advance of osm2psql (to a factor of 4 i am finding).

However I cannot find a supported conversion to mapnik xml/json anywhere.

Has anyone any advice as to how to generate this, or if there's a good online source?


texts:
I am using R and the current 3x version of the lidR package to process point clouds.  Sections of my point cloud are 'unclassified' and I need them to be classified as 'ground'.  (These mis-classifications are happening at interstate overpasses.)  I would like to manually assign these points their 'ground' classification value.  If I were using a GIS, I would select the points using my road polygon and manually reassign their classification value.  Is there a workflow in v3x of the lidR package to do this - to change a point's classification based on, in my case, a polygon shapefile?
In v1.6, there was a command called lasclassify that would allow me to create a user defined attribute in the point cloud and assign my own classifications (as in the lasclassify example shown here https://www.rdocumentation.org/packages/lidR/versions/1.6.1/topics/lasclassify).  I think it has been long deprecated, though.
There is not much to see, but I provide my code below for anybody who is interested:
library(lidR)
library(rgdal)
mypath<- "c:/myPath/"
tdat<- readLAS(file.path(myPath,"myLas.las"))

shppath<- myPath
shpname<- "myRoadPolygon"
shape_spdf<- readOGR(dsn = shppath, layer = shpname)
shppoly<- shape_spdf@polygons[[1]]@Polygons[[1]]

###here I would like to select the points in *tdat* by the polygon shapes in *shppoly*...


texts: How to get unique feature count from a field w.r.t another field using ogrinfo?

texts:
We have some 40 folders which has two - three tab files in each folder. Table structure will be the same for all the tab files. What i would like to extract is Feature Count, LineNo (will be 0,1,2,3,etc..), StationNo (Will be from 0 - 85). LineNo and StationNo are related. For LineNo "0" = StationNo will be "0,1,2 upto 85" and again for LineNo "1" = StationNo will be "0,1,2,3 upto 85" (see the image).

So how to extract these informations using ogrinfo ?


texts: Why does QGIS 2.14 install via apt-get return 403 Forbidden?

texts:
After add the repository http://qgis.org/debian and run

  $ apt-get update

, before try QGIS install, I get the following message:
 http://qgis.org/debian/dists/vivid/main/binary-i386/Packages  403  Forbidden [IP: 104.27.138.220 80].

Obviously, this mean that I can not connect to the repository, but
I can access it through web browser after answer a reCaptcha process. I figure out that QGIS servers are behind CloudFlare.

How can I proceed to install QGIS 2.14 via apt-get ?


texts: How to calculate the Longitude and Latitude of rotated polygons in Openlayers3?

texts:
We are using external Transform library with Openlayers for the polygon scale,Translate and rotate actions.
After rotating the rectangle/square/polygon, when we scale it, the shape is getting stretched and it loses it's shape.
Can anyone share the general formula for calculating the Longitude and Latitude which I can use inside a function using JavaScript?


texts: How not prompt the user to save changes of layers that are on the edit mode?

texts:
When a layer is on the edit mode and I click on the exit button of the QGIS project, it asks me if I want to save changes that I made (automatically). So, how can I disable it using PyQGIS? I've tried QgsProject.instance().setDirty(False) but it keeps asking if I want to save my changes. The image below shows the message that pops up.

How to automatically discard changes and never prompt the user to save changes?
I saw simillar questions like this one: Is there is a signal in QGIS when a project is about to be closed?
But it was also an open question.


texts: Displaying trailing zeros to right of decimal in label using QGIS?

texts:
is there a way I can get zeros to show up in my label when they come to the right of a decimal in a number such a s 270.00? I'm trying to display two decimals in all my numerical labels, but QGIS automatically removes trailing zeros.

I've included a snapshot of my attribute table showing a column that i specified the precision to 2, and it shows 680.00, like I want it. But when labeling, these zeros aren't included.


texts: Converting geopackage to feature class (shapefile or geodatabase) using ArcGIS Pro

texts:
I work with ArcGIS Pro and got a geopackage (.gpkg) from a client. It seems that one cannot directly work with the geopackage data in ArcGIS, so I would like to convert it to a shapefile or feature class.
How can a geopackage be converted to a shapefile or a feature class in ArcGIS Pro?
I neither can use an online converter (file is too big) nor another GIS program.


texts:
I am new to coding and need some help debugging...I am trying to classify a region in Peru using Sentinel 2. The error I am receiving is "'image' is not defined in this scope." I believe my issues is on line 31 and the fact that I am working within a collection not an image...but I am still completely lost. Any tips?
var growth = ee.FeatureCollection('users/gtnatia/just_trees_training_2');
growth = growth.geometry();
Map.centerObject(growth);
Map.addLayer(growth, {color: 'red'}, 'growth');

//

var all_other = ee.FeatureCollection('users/gtnatia/all_other_training');
all_other = all_other.geometry();
Map.centerObject(all_other);
Map.addLayer(all_other, {color: 'black'}, 'all_other');

// Make a cloud-free Landsat 8 TOA composite (from raw imagery).
var sentinel2 = ee.ImageCollection('COPERNICUS/S2_SR')
                    .filterDate('2018-01-01', '2018-12-31')
                    .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
                    .filterBounds(aoi);

// Use these bands for prediction.
var bands = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B10', 'B11'];

//

var polygons = ee.FeatureCollection([
  ee.Feature(growth, {'class': 0}),
  ee.Feature(all_other, {'class': 1}),
]);

// Get the values for all pixels in each polygon in the training.
var training = image.sampleRegions({
  // Get the sample from the polygons FeatureCollection.
  collection: polygons,
  // Keep this list of properties from the polygons.
  properties: ['class'],
  // Set the scale to get Landsat pixels in the polygons.
  scale: 30
});

// Create an SVM classifier with custom parameters.
var classifier = ee.Classifier.libsvm({
  kernelType: 'RBF',
  gamma: 0.5,
  cost: 10
});

// Train the classifier.
var trained = classifier.train(training, 'class', bands);

// Classify the image.
var classified = image.classify(trained);

// Display the classification result and the input image.
Map.setCenter(-72.355,-13.880, 9);
Map.addLayer(image, {bands: ['B4', 'B3', 'B2'], max: 0.5, gamma: 2});
Map.addLayer(polygons, {}, 'training polygons');
Map.addLayer(classified,
             {min: 0, max: 1, palette: ['red', 'green']},
             'supervised_classification');


texts: PostGIS - Count adjacent polygons connected by line

texts:
I have a layer with municipal limits.
With the following code I calculated the number of adjacent municipals:
SELECT distinct a.id, count (b.*), a.geom, a.codistat, a.name
FROM municipal a
JOIN municipal b
ON ST_Intersects(b.geom, a.geom)
where a.id != b.id
GROUP BY a.id, a.geom, a.codistat, a.name

Now I am looking for the solution to calculate the number of adjacent municipalities that are connected by a road (layer way with only id attribute) to understand the value of the connections between municipalities.
In some cases this value will be lower than the one I previously calculated.
I have some problems figuring out how to integrate the code.
Any help is welcome


texts:
I have a DXF with text and I have a shapefile with some squares.

Is it possible, in QGIS, to add the text from the DXF to a new column in the shapefile?


texts: What characters are not acceptable in an ArcGIS text field (File Geodatabase)?

texts:
Trying to join a table to my feature and it keeps leaving out a field. The field is a legal description of a property, it has all sorts of funky characters:

:
;
/
\
()
'
"
and more

The field is also long, from 100 to 1500 characters, but based on this I should be ok for field length: Maximum length of text fields in shapefile and geodatabase formats?

I can't find any documentation regarding what can or cannot be entered here. I assume its only Unicode or ASCII values, but I've removed pretty much everything and am still having the join error. The join error, doesn't show up in the join, it just excludes the field in the join.


texts: Sort a feature table by geographic location

texts:
Here's in interesting challenge: I need to embed a table in a map layout so that the records in the table are sorted according to their geographic position on the ground. For example, in a north-to-top map, the northern-most feature (points in this case) should be first record, and the southern-most the last.

The map composition is in ArcGIS, but I'm willing to use gdal/ogr, qgis, etc. for the table sorting if need be.

To expand this question beyond my immediate use case and thereby make it more useful to the community, a robust solution would allow choosing other geographic sorts, such as west to east or inner to outer.


texts: Link the edges and the nodes with a relationship class table in a geometric network

texts:
I have a bus network in arcgis composed by:

The Bus stops (points)
The Bus lines (multilines)
A relashionship table who link the stops and the line

So the relationship (Many-to-Many) contains a table who looks like that:

[ID_NODE, ID_LINE]
[1      , 1      ]
[2      , 1      ]
[3      , 1      ]
[2      , 2      ]
[4      , 2      ]
...

I'm now trying to create a geometric network based on those features.

When I create the geometric network I can add:

The lines
The stops

But not the many to many relashionship table !

It happens that a line go throught a stop but the bus don't serve this stop, so the line and the stop should be unconnected.

How can I add the information contained in my relationship table in the geometric network ?


texts: Converting multiple CSV automatically to shapefile using ArcPy?

texts:
I have more than 300 CSV files. Every CSV file contains different Longitude, Latitude and other parameters (e.g. Date, Longitude, Latitude, Elevation, Max Temperature, Min Temperature, Precipitation, Wind, Relative Humidity, Solar Radiation).

I need to convert CSV to shapefile using ArcGIS Python to automatically loop through every CSV. I have made one code but it takes only a single file, and does not loop through every CSV file.

Below I copied my script which is creating only single shapefile taking from last CSV file. Not creating shapefile from all CSV files

    import arcpy,os

... shpworkspace = r"D:\SWAT-WEATHER-DATA2"
... arcpy.env.workspace = shpworkspace
... arcpy.env.overwriteOutput = True
...
... csvlist = arcpy.ListFiles("*.csv")
...
... for csvfile in csvlist:
...     print csvfile
...     outlayer = "CSVEventLayer"
...     spatialreference = "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]];-400 -400 1000000000;-100000 10000;-100000 10000;8.98315284119522E-09;0.001;0.001;IsHighPrecision"
...     arcpy.MakeXYEventLayer_management(csvfile,"Longitude","Latitude",outlayer,spatialreference,"#")
...
...     shpfile = os.path.splitext(csvfile)[0].replace('-','_')
... arcpy.CopyFeatures_management(outlayer,shpfile)
... del outlayer
...
weatherdata-233688.csv
weatherdata-233691.csv
weatherdata-233694.csv
weatherdata-233697.csv
weatherdata-233700.csv
weatherdata-233703.csv
weatherdata-233706.csv
weatherdata-236688.csv
weatherdata-236691.csv
weatherdata-236694.csv
weatherdata-236697.csv
weatherdata-236700.csv
weatherdata-236703.csv
weatherdata-236706.csv
weatherdata-239688.csv
weatherdata-239691.csv
weatherdata-239694.csv
weatherdata-239697.csv
weatherdata-239700.csv
weatherdata-239703.csv
weatherdata-239706.csv


texts: How to find locations for hydro power using raster analysis?

texts:
I have a stream network and DEM file both in raster format. I need to extract values from DEM to stream network. For the analysis in stream network after that, I need to find out the places for hydro power based on various criteria of head, stream order.
So how can i start my work?

I am planning to convert both raster file into text file using numpy array and extracting head for all the cells of stream and after that i will start processing. Is that a good idea to start or is there any other way?

I am quite new to python.


texts:
import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, "*")[0]
addLayer = arcpy.mapping.Layer(r"F:\project\Final\Ground floor\steps.lyr")
arcpy.mapping.AddLayer(df, addLayer, "BOTTOM")

this code is for adding a layer in arcmap. when i use the same code to add a 3D layer in arcscene i get the following error

cannot open map document

Failed to execute (addlayer).


texts: Visualization of 3D GIS daten with 3D glasses

texts:
A colleague of mine wants to inspect 3D aerial photos in 3D with 3d glasses. He thinks he just needs the right 3d Graphics card and monitor with more Ghz - something along the lines of nvidia 3d Gaming cards I imagine. The 3d Visualisation could then be produced automagically without any need for extentions in ArcGIS.

Has anyone any experience with this kind of hardware? what works, what doesn´t what should be avoided at all cost etc?


texts:
I'd like to show place labels (cities and towns) only for the one country. Is it possible to achieve this in Maputnik editor? I'm looking for something like Mapbox "iso_3166_1" parameter. I've tried to add iso_a2 to the layer filter settings in such way:
"id": "place_city",
...
"filter": [
  "all",
   ["==", "$type", "Point"],
   ["==", "class", "city"],
   ["==", "iso_a2", "FR"]

but it didn't work.


texts:
I am working with 4 band GeoTIFF images. I used GDAL to open these images. Each layer needs pre-processing and after completion of processing, I need to write back into a TIFF file. I used the below code, but I couldn't get all the bands together. When I check the written file, it contains only a single band.

for data_file in os.listdir(destSAR):
    data = data_file.split('_')
    joinStr = [data[6],data[7],data[8],data[9]]
    joinedName = '_'.join(joinStr)
    imgSAR = joinSAR+joinedName
    imgEO = joinEO+joinedName
    SARt = gdal.Open(imgSAR)
    opref = cv2.imread(imgEO)
    SAR=SARt.ReadAsArray()

    output=getSARmask(SAR,opref) #processing the SAR images
    os.chdir(outputPath)
    filename = 'post_process_'+joinedName
#     print(len(output)) # have 4 bands in output
    for i in range(len(output)):
        ds = SARt #original SAR image with 4 bands, using for Geo information
        arr = output[i]
#         arr = band.ReadAsArray()
        [cols, rows] = arr.shape
        arr_min = arr.min()
        arr_max = arr.max()
        arr_mean = int(arr.mean())
        arr_out = np.where((arr < arr_mean), 10000, arr)
        driver = gdal.GetDriverByName("GTiff")
        outdata = driver.Create(filename, rows, cols, 1, gdal.GDT_UInt16)
        outdata.SetGeoTransform(ds.GetGeoTransform())##sets same geotransform as input
        outdata.SetProjection(ds.GetProjection())##sets same projection as input
        outdata.GetRasterBand(1).WriteArray(arr_out)
        outdata.GetRasterBand(1).SetNoDataValue(10000)##if you want these values transparent
        outdata.FlushCache() ##saves to disk!!
        outdata = None
        band=None
        ds=None


texts: Using Spatial References with ArcPy after *.prj files no longer installed by ArcGIS 10.1 for Desktop?

texts:
I'm writing a python script that uses the CreateFeatureclass_management tool, without a template. One of the required parameters is a spatial reference. I used to be able to point to the system *.prj file in 9.3, but that's no longer an option, and I don't want to point to an existing shapefile's *.prj because that will cause issues when I share the script.

Now that Esri doesn't install all the *.prj files with ArcGIS, how can I set the projection in a script without referencing an existing shapefile?


texts: Displaying degree symbol in popup template of ArcGIS API for JavaScript?

texts:
I am editing javascript code of ArcGIS API for JavaScript for online map.

the piece of code is here:

var lvEsriPopupT = new esri.dijit.PopupTemplate(
                   {
                      "title": "Weather",
                      "outFields": "*",
                      "fieldInfos": [{
                         "fieldName": "Temperature",
                         "label": "Temperature  deg& C",
                         "visible": true},
                         ... ],
                      ...
                  });

The code works fine. the issue that I am trying to fix, is displaying the degree symbol inside double quotes.
anything I tried (° ℃ ℃) comes ups without being decoded to the actual symbol. I suspect it is because it is inside double quotes.

Is there a solution for that?


texts: Create a collection of geometries in JTS

texts:
I'd like to create a program that a user inputs coordinates as points or lines and gets a geojson file.
The problem is I don't know how to store all the geometries in one file.
For example I want to create a point array, then a line array and a polygon array. Then all to be combined in one and written in a Json file.
 A sample code will be like that

import java.util.ArrayList;
import java.util.Scanner;
import com.vividsolutions.jts.geom.*;
public class new {
public static void main(String args[])
 {

 int b;
 Scanner in = new Scanner(System.in);
 System.out.println("Enter number of points");
 b = in.nextInt();
 int arr[] = new int[b*2];

 System.out.println("Enter "+b+" pairs of coordinates");
 for(int i=0;i<arr.length;i++)
     arr[i]=in.nextInt();
 Coordinate[] coordinates = new Coordinate[arr.length/2];
 int i = 0;
  // Go through the arr and add each point as a Coordinate object to the coordinates array.

     while (i < arr.length) {
    // transform string arguments into double values
    double x = arr[i];
    double y = arr[i+1];
    // create a new Coordinate object and add it to the coordinates array
    Coordinate newCoord = new Coordinate(x,y);
    coordinates[i/2] = newCoord;
    Point point = new GeometryFactory().createPoint(newCoord);

    GeometryCollection geome = new GeometryFactory().createGeometryCollection(null);
 // Here is the problem. How to write every point in a geometry collection?

    i=i+2;

  //The same will apply for line and polygon.
  //Then all combined geometries could be written in a Json
  //  GeometryJSON g = new GeometryJSON();
    try {
        g.writeGeometryCollection(GeometryCollection, "final.json");
    } catch (IOException e) {
        // TODO Auto-generated catch block
        e.printStackTrace();
    }


texts: To get attribute data of polygon layer by giving Lat/Lon

texts:
Is anyone having an idea about how to get attribute data of a polygon layer by giving Coordinates with in that polygon. Here i have tried with these codes..

SELECT * FROM kar_bbmp_198 WHERE ST_Contains(geom, ST_GeomFromEWKT('SRID=4326;POINT(779322 1440590)'))

SELECT * FROM public.kar_bbmp_198 WHERE ST_Within(geom,ST_PointFromText('POINT(779322 1440590)',32643));

But I am getting an error message like

ERROR: Operation on mixed SRID geometries ********** Error ********** ERROR: Operation on mixed SRID geometries SQL state: XX000

Any idea... I am using PostGIS 2.0 and PostgreSQL 9.1


texts: How to edit using a WFS service, without showing it on the map in OpenLayers?

texts:
I have a service with about thousand points. I wish to enable the user to edit the attributes of the features.

Since there are so many points, I don't want to display them as a vector layer (using WFS Protocol) on the map. I am using a WMS layer to display it on the map.

Using a WMSGetFeatureInfo control, I am getting the attributes of a selected point and showing them in a form to the user.

How can I now save the edited attributes to the database?  Is there some other way of achieving it?


texts: How to join shapefile to CSV (not other way around)?

texts: I have a CSV with multiple rows per county (it contains data on number of coronavirus cases per week). I want to make a choropleth map showing the change over time, but when I join this data to a county shapefile in QGIS, it only preserves the most recent datapoint. How do I execute a join such that all of the rows in the time series CSV are preserved, but have a geometry attribute? Ideal output is a geojson.


texts: Draping surface on custom data rather than built in elevation data for 3d web scene in CityEngine webviewer app?

texts: When I publish a web scene from ArcGIS Pro to ArcGIS Online for use in a cityengine 3d web scene, all of my surfaces that I have custom elevations assigned to come out flat. As far as I can tell the only elevation data that the 3d webscene will read or apply, is the built in esri elevation data.   I'm trying to display subsurface geologic data.


texts:
Last week I tried everything to build up a tile server but the Tile Server I was using (tileserver-gl) doesn't work for the 2012 r2 WS Server.

After that i tried it with leaflet but I couldn't figure it out. And I also tried it with docker but docker requires the WS Server 2016.

I need a tileserver there you can edit the style and using your own .osm, where the map is stored in. What i tried: tileserver-gl, leaflet (can't use my own style and map), brouter, mapbox, mapnik.

How can I edit the style from the map so that I can make the water red or green?

I need a tileserver using my own .osm files and also to configure my own style.


texts: So I have numbers that range from 5 to above 50, and when I sort these from least to greatest, the numbers less than 10 are treated as higher than anything else.  It also does the same when I try to symbolize by category, and treats the less than 10 numbers as higher than the actual higher integers.  I am a total noob an have no idea what I am doing in QGIS or how to begin to solve this problem.  Thanks for reading.


texts: OSRM Table service returns empty array when I use a HTTP Server

texts:
I am currently trying to use OSRM's Table service to get the durations of the fastest routes between two points, to be able to build a time matrix using Python and later use it as a parameter in an optimization problem.

My problem is that whenever I run my own server through my jupyter notebook (using five different locations) the 'durations' array returns empty, like in the following example:

str_locations = '-30.59464,-71.19366;-30.60386,-71.21342;-30.58654,-71.18429;-30.60298,-71.20075;-30.58603,-71.19196'
r=requests.get('http://localhost:5000/table/v1/driving/'+str_locations+'?sources=0')
r.json()['durations'][0]

returns:

[0, 0, 0, 0, 0]

I've already checked that those points exist and are connected using the two first points in map.project-osrm.org (here's the link). Plus, whenever I replace the local server with the online demo (http://router.project-osrm.org/) the array is not empty, but that doesn't work with more than 100 coordinates and I need at least 250.

I think it could have something to do with the algorithm that I am using, MLD rather than CH?


texts:
In QGIS I used the Georeferencer on a jpg image I had prepared in Gimp.  But when I load the resulting tiff image into my QGIS project it is invisible, even though an image viewing program can show it perfectly well.  I used exactly the same method to georeference a jpg image of an OS sheet covering the same area, and that shows up perfectly well in QGIS.

Why is my second layer invisible?


texts: How to find the destination point given a start point, distance and bearing?

texts:
If I have a starting point as well as distance and a bearing to the destination point, how can I find the coordinates of the destination point using PostGIS?

I'm looking for a PostGIS equivalent to this function:
http://www.movable-type.co.uk/scripts/latlong.html#destPoint


texts: Correct syntax for retrieving features from GML format in OpenLayers

texts:
Yesterday I made with qgis the shapefiles of my town and now I wanted to add them in a web application I'm making with OpenLayers. The problem is that, while it loads them perfectly with the wms service, the vector layers fail to get the feature from the xml/gml file and I don't know why. I don't even get any error. Here is a some part of my code:

OpenLayers.ProxyHost = 'proxy.php?url=';

    var map = new OpenLayers.Map("map", {
        //projection: new OpenLayers.Projection('EPSG:900913'),
        projection: new OpenLayers.Projection('EPSG:3857'),
        units: 'm',
        maxResolution: 4.777314266967774,//156543.0339,
        maxExtent: new OpenLayers.Bounds(-1718.3524101391, -3513.9879361762, 3058.9618568286, 785.5949040948),
        numZoomLevels: 4
    });

var mroads = new OpenLayers.Layer.WMS(
                 'Major Roads',
                 'http://andri.ns0.it:8080/geoserver/Grottole/wms',
             {
                 layers: 'Grottole:MajorRoads'
             },
             {
                 isBaseLayer: false,
                 opacity: .9
             });

var mRoadsVector = new OpenLayers.Layer.Vector(
            'Major Roads Vector',
        {
            projection: new OpenLayers.Projection('EPSG:3857'),
            protocol: new OpenLayers.Protocol({
                url: 'mroads.xml',
                format: new OpenLayers.Format.GML()
            }),
            strategies: [new OpenLayers.Strategy.Fixed()]
        }
        );

I believe I made everything good so far, so I really hope you can help me.


texts: Linear fishnet grid over layer with latlong coordinates in QGIS

texts: In QGIS I have a city boundary layer which I want to transform into a simple fishnet grid with 30m spacing. The boundary layer is in geographic degrees and the only solution I found so far is choosing a projected coordinate system which produces a rotated and uneven grid. What could I do to overcome this?


texts: Exporting multiple feature classes with joined fields using ArcPy?

texts:
I have multiple FGBs (10) with many feature classes (about 30). Most of the feature classes are joined to a table which is specific to each feature class (i.e. 30 feature classes to 30 tables). I would like to make the joined fields permanent so that I can merge the datasets as required.

The manual process to do this would involve r-click > 'export data' for each feature class in the TOC but this is not feasible given the number of datasets.

Based on esri help pages and other SE threads I started with this code to copy features:

# Copy features from a FGD to FGD

import arcpy
import os

# Set environment settings
arcpy.env.workspace = "C:/GIS/Groundwater_Inventory/Data/GW_Theme_Recharge/Study_Areas/Recharge_Summary_areas.gdb"
arcpy.env.qualifiedFieldNames = False

# Set local variables
out_workspace = "C:/GIS/Workspace/Inventory_data/Recharge_simple.gdb"

# generate a list of features in the workspace
fc_list = arcpy.ListFeatureClasses()

# Execute CopyFeatures for each input feature
for feature in fc_list:
    out_featureclass = os.path.join(out_workspace, os.path.splitext(feature)[0])
    arcpy.CopyFeatures_management(feature, out_featureclass)

Images of tables and feature classes attached. I suspect I need to recreate the join in this script but cannot find any threads which describe this process.

I am using ArcGIS, with entry level python/programming experience.

EDIT: I resolved this issue by copying the required data in the tables into a single csv file using a VBA script based on this help page: https://www.thespreadsheetguru.com/the-code-vault/2014/4/23/loop-through-all-excel-files-in-a-given-folder

I then merged all feature classes and from there I could perform a simple attribute join to the new table.


texts:
I have a QGIS 3.12 project which covers a large area and has 25 layers of different types.  Each layer has its own saved style.
Occasionally I am required to create a "sub-project" which forms a smaller area of the large project.  After creating the clipping polygon to define the boundaries of the smaller project, I use the Clip Multiple Layers plugin to create 25 clipped layers which come into existence with some type of default styling.  Is there some way to avoid the time-consuming task of applying the saved styles individually to each of the 25 clipped layers?


texts: Choropleth in R: custom breaks and plotting

texts:
I have a question with regard to plotting a choropleth map.
What is recommended with regard to plotting a map with custom breaks? And make it look a bit attractive.

In my script I have used the plot and spplot command to visualize some data but I am not entirely happy with the result.
It is easy to see why. Using spplot I get a colour ramp that repeats colours for different values which is a bit weird.
Using plot with the custom breaks I have not been able to make better intervals for the values, such that for instance zero-values are white on the map. The skewed proportions between legend and map are I think probably the result of Rstudio.

So I wondered whether someone has better suggestions for plotting choropleth maps in R using custom breaks?

EDIT 20-02-2013

As suggested I have used ggplot2 for creating a choropleth map. I've used the code below. I encountered two problems. One is that after fortifying the data not all data from the shapefile is transferred to the data frame. So I have to rename the "id" variable such that I can merge the data. Minor issue.

A bigger issue is that despite the fact that the data in the dataframe is correct, the ggplot shows a map where the data is incorrect. The colour value according to the legend does not correspond with the value of the variable in the data. Am I overlooking something here?

# Plotting polygon shapefiles (try-out)
# Load the packages
require(rgdal)
require(maptools)
require(ggplot2)
gpclibPermit() # required for the fortify method
# Load the shapefile
africa=readOGR("/home/GIS",layer="africaII")
names(africa)
# Load the data on conflict
conflict<-read.csv("africa_conflict.csv", header=TRUE)
names(conflict)
# Merge the data together
africadat<-merge(africa, conflict, by="ISO3")
# Changing the data in the shapefile
africa@data <- africadat
# Fortify so that ggplot can plot
africa.points = fortify(africa,region="ISO3")
# Change "id" to "ISO3"such that the datasets can be merged again
names(africa.points)[names(africa.points)=="id"]<-"ISO3"
africa.df=merge(africa.points,africadat,by="ISO3")
# Plot the data
ggplot(africa.df) +
aes(long,lat,group=group) +
geom_polygon(aes(fill =onset))+
geom_path(color="white") +
coord_equal()

Link to replication data

Update 13-11-2015:
Fixed link to shape file.
See answer below for solution on how to plot the choropleth.


texts: Anti-aliased rasterization of vectors in QGIS?

texts:
I have some polygons that I want to rasterize, but in a fairly low-resolution grid (cell size 250x250m). The raster values of the area inside the polygons should normally have the same value (100), but most of the polygons are fairly small (say max 50 m wide), which means that right now they often don't even show up in the final raster (using gdal:rasterize or saga:rasterize). What I would have liked is for the raster value to have been interpolated, in the sense that the final pixel value would have been lower than 100 (say 30) but still present in the output raster. Is there any stock algorithm in QGIS capable of this?
I guess that another way would be to first do a much higher-resolution bi-valued rasterization, and then downsample the resulting raster with an algorithm that hopefully handles this. This would probably be slower and require more memory.


texts: How do I build a new raster based on raster stats in QGIS?

texts:
I want to create a new raster that uses the mean value of the original raster as the base class of the image and then uses a +- 1 and 2 standard deviations as the other four classes of the image.  Is there a simple way to do this in QGIS?

Here is the equation I used before posting my question. What am I doing wrong? 2018712@5 is the band I'm working with. 194.26 and 9.8663 are the mean and standard deviation respectively.

("2018712@5" >= (194.26+(9.8663*1)))*1 + ("2018712@5" >= (194.26+(9.8663*1.5)))*1.5 + ("2018712@5" >= (194.26+(9.8663*2)))*2 + ("2018712@5" <= (194.26+(9.8663*-1)))*-1 + ("2018712@5" <= (194.26+(9.8663*-1.5)))*-1.5 + ("2018712@5" <= (194.26+(9.8663*-2)))*-2