texts:
I'm writing ArcPy code in ArcGIS trying to calculate a list of rasters.

I get an error message that says:

  TypeError: unsupported operand type(s) for -: 'unicode' and 'unicode'.

Here's my code:

    import arcpy
    from arcpy import env
    import os

    arcpy.env.workspace = "Z:\\comp1\\yitong\\test\\"
    arcpy.env.overwriteOutput = True

    B5 = arcpy.ListRasters("*BAND_05*","")
    B4 = arcpy.ListRasters("*BAND_04*","")

    waterVapor = []

    for rasterB4 in B4:
        for rasterB5 in B5:
            print rasterB4
            print rasterB5
            partName = str(rasterB4[16:30])
            waterVaporOutputName = "waterVapor" + str(partName)
            waterVapor = 0.24+0.503*(rasterB4 - rasterB5)
            waterVapor.save(waterVaporOutputName)


texts: Custom SVG world map usable for react-simple-maps (topojson, d3-geo)?

texts:
A client has supplied a "custom" SVG containing a very detailed world map that includes id's for all countries. I'm not sure how "custom" it is.

Is it at all possible to convert it to topojson, which is needed to use it with react-simple-maps?

https://www.react-simple-maps.io


texts: Using Data Driven Pages to display Attribute table values?

texts:
How can I show ID/page specific attributes on a map book project that has already been made?

I recently took over our GIS mapping duties for a former coworker who retired.  He had already made a map book project using python for Food plot plantings, but there are no Attributes listed on each map page.

I am working with ArcGIS Desktop 10.2.


texts: QGIS OpenCL settings fall back to onboard-graphics

texts:
I am running QGIS 3.4.14 on my Mac 16', with a Radeon 5500 GPU.
I've enabled the OpenCL support in QGIS settings as well as my user profile settings.

It offers me the CPU, the built-in oonboard GPU Nvidia 630 as well as the Radeon 5500.

Unfortunately, when selecting the Radeon and shutting QGIS down and start it again it always selects the NVIDIA and not the more powerful Radeon.

Is it a bug or a feature?


texts:
This is my situation:

I have a table A in schema A which includes personal data. I want to give specific persons the permission to see some data from that table in QGIS, excluding that personal data. Therefore I have a view B of table A in Schema B which doesn't include the columns of the personal data. The problem is that I want the person to edit two columns in that view. So I created a trigger for this view and gave permission to the person to update these two columns in it. However, I also need to give him permission to the table A. My problem is now that the person is able to access table A in Schema A as he needs the permissions to update the two columns. Because that person doesn't know anything about the structure of the database, this might be fine as he would not know that there is a table behind that view. However, when connecting to the database in QGIS he can also see schema A and is able to load the table A.

Perhaps there is a way to limit the acces to table A in the backend. But for now it would be enough when the schema A could be hidden in QGIS. Is that possible? Or does anyone has a better idea for a solution?


texts: How to to validate geometries using ogr2ogr

texts:
I have been trying to dissolve a shapefile with more than 38000 rows. I used query as follows

ogr2ogr -f "ESRI Shapefile" dissolved.shp Barangays.shp -dialect sqlite -sql "select ST_union(Geometry),* from input  GROUP BY concated"

It gave me error that geometry is invalid because of self-intersecting polygons.

I found out there are 86 features which have invalid geometries using QGIS. I tried Make Valid tool as suggested in this question. But this is also not working, giving error undefined symbol, which I couldn't understand. So I decided to validate using ogr2ogr and after searching through google I tried this query

ogr2ogr -f "ESRI Shapefile" valid_barangays.shp Barangays.shp -dialect sqlite -sql "select Geometry.MakeValid(),* from Barangays"

It gave me this error

ERROR 1: In ExecuteSQL(): sqlite3_prepare(select Geometry.MakeValid(),* from Barangays):
  near "(": syntax error

How can I solve this problem?

I also imported my shapefile using ogr2ogr in phpmyadmin mysql. It saved geomtry column as [geometry some_number B]. I am not able to query. I always give error function not exist.


texts: Initializing Google Earth Engine fails in web application

texts:
I am trying to call the ee.data.authenticateViaOAuth() function which then calls ee.initialize().
I follow the instructions line for line but I end up with the following error:
Initialization error: TypeError: Cannot use 'in' operator to search for 'Classifier' in undefined

Has anyone figured out a way around this?


texts: Duplicate any active layer in memory using pyqgis?

texts:
I have seen Duplicating layer in memory using pyqgis?
and I would like to do the same but for any active layer in the TOC whether it is a Polygon, Linestring or Point layer.

Is it possible?

This is the code I tried

 iface = qgis.utils.iface
vl = iface.activeLayer()
feats = [feat for feat in vl.getFeatures()]
mem_layer = QgsVectorLayer(vl.source(), "duplicated_layer", "memory")   mem_layer_data = mem_layer.dataProvider()
attr = layer.dataProvider().fields().toList() mem_layer_data.addAttributes(attr) mem_layer.updateFields() mem_layer_data.addFeatures(feats) QgsMapLayerRegistry.instance().addMapLayer(mem_layer)


texts: QGIS - convert multiple polygons to evenly spaced points within polygon boundary

texts: I need to convert the polygons in a vector layer to evenly spaced points within the polygon boundaries in QGIS. The vector file contains >1000 polygons mostly with overlaps.   The attribute data of the polygon needs to be captured in the new point table. Only extraction of the nodes is not enough.


texts:
As we all know, the attribute table field properties has number format to display number of decimal places to up to maximum 15.
If we have many fields/columns in attribute table, how can we change all fields/columns decimal places roundoff/display at a time ( I mean in single operation. Not changing each and every field properties).

Ex:

ColA:12.987654 ! ColB:34.20987766  ! ColC:89.654309877

In one operation my output should be (Rounding/display to 2 digits)

For Rounding:

ColA:12.99 ! ColB:34.201  ! ColC:89.66

For Display:

ColA:12.98 ! ColB:34.20  ! ColC:89.65

is there any method/tool to do this in Arcgis?


texts: Different trend test result with modifiedmk and zyp package in R

texts:
I am running two sen slope trend tests in R, the first with the modifiedmk package and the second with the zyp package.  I don't understand why with the zyp package the resulting raster raster.mkt2 is empty.

This is my code:

library(Kendall)
library(zyp)
library(gimms)
library(raster)
library(rgdal)
library(modifiedmk)
library(rkt)

#create a fake raster
r = raster(nrow = 50, ncol = 50)
r[] = sample(seq(from = 1, to = 1000, by = 1), size = ncell(r), replace = TRUE)

#stack it
final_stacked = stack(r,r*2,r*3,r*4)

#function with modifiedmk
tsfun <- function(x, na.rm) {
  if(all(is.na(x))){return(NaN)}
  return(pwmk(x)[["Sen's Slope"]])
}

#function with zyp
tsfun2 <- function(x, na.rm) {
  if(all(is.na(x))){return(NaN)}
  return(zyp.trend.vector(x)[["trend"]])
}

raster.mkt1 <- calc(final_stacked, fun=tsfun, na.rm = TRUE)
raster.mkt2 <- calc(final_stacked, fun = tsfun2, na.rm = TRUE)

plot(raster.mkt1)
plot(raster.mkt2)

EDIT:

changing this line final_stacked = stack(r,r*3,r*8,r*9, r*20, r*40, r*55, r*66) so that it isn't linear still leads to different results.

#mean of modifiedmk output
m1 = cellStats(raster.mkt1, stat = 'mean', na.rm = TRUE)

#mean of zyp output
m2 = cellStats(raster.mkt2, stat = 'mean', na.rm = TRUE)

returns 2696 vs 5252

To further clarify lets pull a random vector from final_stacked

t = final_stacked[10,10][,]

returns:

layer.1 layer.2 layer.3 layer.4 layer.5 layer.6 layer.7 layer.8
    841    2523    6728    7569   16820   33640   46255   55506

pwmk(t) returns:

Z-Value.layer.3      Sen's Slope old. Sen's Slope  P-value.layer.3        S.layer.3           Var(S)      Tau.layer.3
    2.703381e+00     4.512682e+03     8.788450e+03     6.863795e-03     1.900000e+01     4.433333e+01     9.047619e-01

and zyp.trend.vector(t) returns:

  lbound         trend        trendp        ubound           tau           sig         nruns       autocor    valid_frac
 3.994750e+03  8.788450e+03  7.030760e+04  1.198425e+04  9.047619e-01  6.863832e-03  1.000000e+00  4.535428e-01  1.000000e+00
       linear     intercept
 8.229786e+03 -1.717743e+04

So they match, but raster.mkt1[10,10] returns 4512 while raster.mkt2[10,10] returns 8788


texts: Efficiently iterating through raster calculations with Rasterio or GDAL

texts:
I've got a very slow and round-about way to currently do this, but I struggle with the general concept of using windows in raster calculations and how to apply an affine or geotransformation between two large rasters. I'm hoping someone can help clarify what the heck i'm doing wrong here.

The following code is aimed at taking a large raster (200mb) and minusing a DTM from the raster values. The main points below:

I have these:

DTM - NED raster tiles merged to bigtiff for the whole USA (77gigs)
Water - Water depth raster over a portion of the US (water + DTM) (~200mb)

I want to perform a simple Water - DTM

Currently have a script that resamples a raster by individual columns where I query 1 point at a time to get DEM value at each point location in the raster. I know this probably isn't ideal but it's the only way I know to get a water depth raster efficiently.

The follow is my attempt at using Rasterio windows and then iterating through the pixels of each window to minus the DEM. At the end of the day i'd like to have another raster, but also a flat file with latitude/longitude/depth.

Most of this code came from here

with rio.open(tst) as src:
    gt = src.get_transform()
    a = src.affine
    prj = Proj(src.crs)
    a1 = a * Affine.translation(0.5, 0.5)
    rc2en = lambda r, c: (c, r) * a1
    nx, ny = src.width, src.height
    print "Width: "+str(nx)+" | Height: "+str(ny)
    for ji, window in src.block_windows(1):
        r = src.read(window=window)
        cols, rows = np.meshgrid(np.arange(r.shape[1]), np.arange(r.shape[0]))
        eastings, northings = np.vectorize(rc2en, otypes=[np.float, np.float])(rows, cols)
        longs, lats = transform(prj, prj, eastings, northings)
        for i, x in enumerate(list(r[0][0])):
            if x > -100:
                lat = lats[0][i]
                lon = longs[0][i]
                dem = xy2ras(lon, lat, dtmgt, bnd)
                print "Lat: "+str(lat) + " | Long: "+str(lon)+" | Flood Height: "+str(x)+" | DEM: "+str(dem)
            break


texts: How do I use Multispec for Landsat data?

texts:
It seems that I only can use the Clustering if I want to cluster only one band data..
If I open up more than 1 band (in one image), the cluster toolbar is disabled?

optional:
Is there any free Satellite Imagery processing tools other than this?


texts:
The Raster to Polygon Conversion results in strange polygons. I must have changed something in the settings, but have no idea what it could be.

Does anyone can help?

The red line is the polygone outline


texts: How to get a "classical" globe projection

texts:
I have a shapefile of the world continents in EPSG 4326.

To what reference system do I have to reproject my file to get a projecting like this globe below?


texts:
I'm using the Sieve command within QGIS (Raster > Analysis > Sieve), but when I run the command with any variation of parameters, it completely changes the histogram of my raster data from between -0.140351 to 0.780933, to a binary output of 0 and 1.

This example was run with a threshold of 2 and pixel connections of 4, on a NDVI image calculated from Landsat 7 imagery.


texts: ArcMap Model Builder Layers and Variables

texts:
I am new to Model Builder in ArcMap.

I have made a model to automatically calculate Land Surface Temperature (LST) by using Landsat 8 imagery.
The issue is that the model doesn't fully run towards a result and I presume that it is due to the fact that when I've made the model, the name of the original file remained in the core of the model, so the name of the file needs to be the same if I want the model to give me results. This means that I have to rename every input data to fit the original name for the files.

How can I make the model to run without renaming the files? I remained stuck at this point.


texts: How to loop through features in a layer, zoom to that feature and export to png

texts:
I have a large number of lines in a shapefile representing sites on a map. Each has a unique 2- or 3-digit Id, but they are not consecutive (i.e. 6, 27, 53,... not 1,2,3,...). I want to zoom to each feature in turn and export the map as a PNG file.I want to save each one as the name of the Id.

I have tried the loop suggested by @Hornbydd but this creates many blank/useless maps for the numbers which are not allocated to sites (i.e. 1,2,3 etc are unneeded), which I had to go through and delete. Note that I cleared the selection before exporting as .png because I don't want the line to show up as highlighted blue.

mxd = arcpy.mapping.MapDocument("CURRENT")
for x in range(6,200):
   arcpy.SelectLayerByAttribute_management("sites","NEW_SELECTION",' "Id" = x ')
   df = arcpy.mapping.ListDataFrames(mxd)[0]
   df.zoomToSelectedFeatures()
   df.scale=5000
   arcpy.RefreshActiveView
   arcpy.SelectLayerByAttribute_management("sites","CLEAR_SELECTION")
   arcpy.mapping.ExportToPNG(mxd, r"X:\site maps\site_" + str(x), "PAGE_LAYOUT", 1280, 960, 160)

This worked for the sites I needed but took too long and saved too many useless files.

So then I tried Data Driven Pages:

mxd = arcpy.mapping.MapDocument("CURRENT")
for pagenum in range(1, mxd.dataDrivenPages.pageCount + 1):
...    mxd.dataDrivenPages.currentPageID = pagenum
...    ID = <-- what goes here?
...    print "Exporting site {0}, page {1} of {2}".format(str(ID), str(mxd.dataDrivenPages.currentPageID), str(mxd.dataDrivenPages.pageCount))
...    arcpy.mapping.ExportToPNG(mxd, r"U:\Mapping\site_" + str(ID), "PAGE_LAYOUT", 1280, 960, 160)

I can't work out what to allocate to ID to get it to be the value in the Id field. I have allocated the Id field as the Name field, but I haven't figured out how to get those values out.

Using

...    ID = mxd.dataDrivenPages.pageNameField.name

sets ID to "Id"


texts: how to configure an arcgis server service to enable multiple editing

texts:
I'm using ArcGIS Server 9.3.1, ArcSDE, VS.NET 2008 and I'm trying to enable a service to be edited multiple times, the service its going to be public and it will be used to report emergencies and other situations that will require attention. I need some pointers on how to accomplish this.

Is it a good idea to have an editable map or its better to save the points coordinates to a database and in a report like app generate all the points at once?


texts:
I am using Raster Calculator in ArcGIS 10.4.1 in order to substract two raster surfaces. However, when I substracted the two surfaces an error 000539 was appeared! I tried to adjust the two surfaces in a specify extent, but clip and extract by mask tools didn't work as I wanted. I tried to fit the first surface into the second one, however, one more column was remained in the former surface.
Any suggestion?


texts: How to select attributes by unicode field value?

texts:
I have an attribute table in QGIS 2.10 with a string field (myfield) and when I want to make a selection based on the unique values of this field I get no selection at all.

>idxmyfield = mylayer.fieldNameIndex("myfield")
>uniquevalues = mylayer.uniqueValues(idxmyfield)
>uniquevalues
>[u'avalue', u'bvalue']
>type(uniquevalues[0])
><type 'unicode'>

>y = len(uniquevalues)
>for x in range(y):
     it = mylayer.getFeatures(QgsFeatureRequest().setFilterExpression  (u' "myfield" = {0}'.format(uniquevalues[x])))
     mylayer.setSelectedFeatures([f.id() for f in it])
     selection = mylayer.selectedFeatures()

If I use another field with integer values the code is working. But if the values of the field are unicode/string it is not selecting anything.

So the type unicode differs from string or is it a type of string? If so, how can I change the type from unicode to string? Or how do I need to change the expression in order to get the selection I want?


texts:
Tables that are created in Oracle using ArcGIS Desktop are given an OBJECTID field by default. The OBJECTID field is assigned a unique not-nullable index.

Does the OBJECTID in an Oracle Geodatabase table have a primary key constraint, by default?

I ask, because there can only be one field with a primary key constraint per table. If the OBJECTID does not have a primary key constraint by default, I plan to add a primary key constraint to a different field -- a field called ID. I will use the ID field in conjunction with a foreign key to enforce referential integrity at the database level. More info here.

Note: A primary key constraint is not to be confused with a unique not-nullable index, a candidate key or a so-called relationship class "primary key".


texts:
I once found a great but simple webpage that has a drop down with several text-based GIS formats such as GeoJSON and WKT, and a text box to copy/paste your data.  The feature is then drawn on a world map.  You can continue to add shapes on the same map.

I've lost this page, and I can't find it.  It was super-useful, but I forgot to bookmark it.

So, does anyone know the page, or any page where I can paste WKT and GeoJSON to see them displayed on a map?


texts: Filtering layers in QGIS atlas while displaying them in the canvas

texts:
I would like to create and atlas that instead of displaying all the features within a layer, it only displayed those related to the atlas.

Certainly, I could do as this excellent answer suggests, which is create a rule-based style for the desired layer and create one rule with an experssion like  'field1' =  @atlas_pagename).

Unfortunately, IMHO, this approach has two problems that I would like to overcome:

No features will be displayed in the canvas outside the atlas
It cannot be combined with graduated styles.

In order to solve that, I tried to copy a layer and use a filter with the expression 'field1' =  @atlas_pagename and then created a theme that I would use in the layout's map. Therefore, I could use the original layer to edit and visualize from the canvas, and the copied and filtered one for the atlas. Unfortunately, although that makes sense to me, it does not work.

So my question is: is it any way to filter features in QGIS Atlas while being able to see them all in the canvas?


texts: Creating a large GeoTIFF from multiple rasters

texts:
I have 1 million hectares of area and 1m raster data in ~200 tiles covering it. I am looking to reclassify the raster data and then create one image of the data in GeoTIFF format.

How can I accomplish this?

I am hoping that with the small raster cell size that the TIFF will still be producible and later usable.


texts: Querying indexed vector layer (finding nearby point)?

texts:
I can find nearby features in a layer by using QgsSpatialIndex: Finding the nearest line to a point in QGIS? Using this method, we need to feed in all existing features and query it with a point to get back the nearby features.

However I wonder what is indexed vector layer, the option that can be enabled when creating a new in-memory layer. I thought that it would perform indexing for all features and would give better performance during querying. Therefore, there must be a way to issue a nearby query on the layer itself without explicitly initialising QgsSpatialIndex, or else why should I turn on the index option?

I can turn on index in Layer Properties > General > Create spatial index or it will be turned on by default if I create a temporary scratch layer. I notice that it improves the speed of rendering.

Now in my PyQGIS script, how can I take advantage of the existing index without creating a QgsSpatialIndex class to query nearby features or doing any nearest neighbor analysis?

spIndex = QgsSpatialIndex()
for pointFeatures in pointsLayer.getFeatures():
    spIndex.insertFeature(pointFeat)
nearbyPointFeatureIds = spIndex.nearestNeighbor(testPoint, 3)

My concern is why I should create my own in-memory index again in the first three lines above. Have I missed that I can get something like a QgsVectorLayer.spatialIndex.nearestNeighbor() from the layer itself?


texts: Circles preview in QGIS for TripAdvisor API

texts: I'm going to use the API provided by TripAdvisor (see here: https://developer-tripadvisor.com/content-api/documentation/) to download hotels in some municipalitys from their servers. Data is requested by specifying coordinates of a point (centre of a circle) and the radius of the circles itself. I wanted to know if QGIS allows to draw a preview of these circles over the region object of studying and, if yes, in which way.


texts: Calculating the area of each cell of a Dataarray (or Raster) using Python?

texts:
I would like to compute the area of each respective cell of a raster.

The pixels are in degree coordinates, so it is not trivial to calculate the area of each one "manually".

In R I know there is a function for this:
https://www.rdocumentation.org/packages/raster/versions/2.6-7/topics/area

Does anyone have an idea of how can I do it in python?


texts: Consuming secured ArcGIS for Server services in Google Earth?

texts:
I am able to easily consume unsecured ArcGIS for Server (AGS) services in Google Earth (GE), but not secured services.

I have created a KML Link in AGS Manager with no luck, and added the service manually as a link in GE with the same result.

The service is added to GE but a red circle is displayed when the service refreshes.

Has anyone been able to do this successfully?


texts: QGIS using spatial query in the Layer Attribute Editor

texts:
I am currently setting up a table to capture polygons into a shapefile to show grounds maintenance areas.
One of the columns in the layer is the parish name. This is contained in another layer

In the grounds maintenance layer I can set the parish column to 'Value Relation' and specify the Parish layer and the relevant Key and Name column.
If I then capture a new record, I can then add the parish name. But I would need to have queried the parish layer first to have established the Parish name.

Is there a way of automating this process so that when I capture the new record into the Grounds Maintenance layer, the parish name is automatically populated.
I know I can do a spatial join but I'd rather not have to do this as it would mean an additional step which is seperate from the main task of capturing the data.


texts: Resample 20m bands in Sentinel-2 image to match 10m bands (GEE)

texts:
Is there a standard way to harmonize the resolutions in Sentinel-2 images? I want to work with the most 10m bands but add two 20m bands in the analysis. However, before running the analysis I would like to resample the 20m bands to match the 10m resolution bands.
I could probably do it by subsetting the 20m bands, resample them and add them to the 10m band subset. However, there might be a better and standardized way?
var image = ee.Image('COPERNICUS/S2/20190630T100031_20190630T100212_T32TQM')
   .select(["B2-4", "B5"])

What's the best way to get B5 to the same spatial 10m resolution (replacing the B5 20m resolution band)?
https://code.earthengine.google.com/70788a7c5f2e49230f590bb528cf55a8


texts:
I have facing a problem on arranging the SAM Classification on ENVI. I applied the Landsat 8 and want to classify the area using SAM on my private libraly, thus the system notice that for band 9 to 11 the bias and gain were lost. I tried to input them but do not know where to do and the system still define that they were still missing.

What can i do for this issue?


texts: QGIS: georeferencing report: tail cutted off

texts:
I'm using QGIS 1.8.0-2 to georeference some old maps and want to report residuals and other information via the PDF report function of the Georeferencer tool. If there's too much information, it's cut off at the bottom of the PDF page.

Is there a way to get all the information (i.e. on a second page)?

Regards,

Thomas


texts: Possible/Feasible to download parcel map and add up to 20 layers

texts:
Novice who just downloaded and installed QGIS a few minutes ago. Currently reading about the functionality of QGIS but thought that I may be able to obtain immediate insight in reference to my objective.

Is it possible to use QGIS to download parcel maps from a internet site and then add up to 20 custom layers?

Specifically, my objective is to download a parcel map for a particular city with "basic" attributes such as:

Parcel Outlines
Zip codes
Hydrography
Highways/Freeways

Then, I would like to have the ability to use the data within my MS Excel worksheet (worksheet 1) to create up to 20 custom layers for display on the parcel map and maybe color code the layers.

Layers could include metrics such as population, average age, median sales, etc.

Finally, I would like to display this map on a worksheet within my MS Excel workbook for presentation.

Is this even possible?


texts: How to get extent of Image layer in OpenLayers?

texts:
Considering the code below (which uses openlayers and Viglino's ol-ext), how do I get the extent of the geoimg?

var map = new ol.Map({
  target: 'map',
  view: new ol.View({
    zoom: 18,
    center: ol.proj.fromLonLat([4.8728,52.3316])
  }),
  layers: [new ol.layer.Tile({
    source: new ol.source.OSM()
  })]
})

var geoimg = new ol.layer.Image({
  name: "Georef",
  opacity: 0.5,
  source: new ol.source.GeoImage({
    url: 'image.jpg',
    imageCenter: ol.proj.fromLonLat([4.8728,52.3316]),
    projection: 'EPSG:3857',
    imageScale: [1, 1]
  })
})
map.addLayer(geoimg);

Here is a jsfiddle, which shows how a 160x160 pixels geoimage, with an imageScale of [1,1] measures 100 by 100 meters on the map. You can measure it yourself, but I placed the picture on the map between 2 streets that I know are 100 meters apart.

Things that do not work:

geoimg.getExtent()

it correctly says undefined because getExtent only gets the extent that was set for the Image Layer, which is a masking extent according to the OpenLayers documentation.

I also tried to look into geoimg.values_.source.canvas_.extent, which obviously returns the extent of the canvas, i.e. the div with the id="map"

EDIT1: After the discussion I had with @TomazicM in the comments, I noticed that my question might not be clear. The picture I am overlaying is not a map picture as in Vigliano's example, but a regular picture. I don't want to find out the coordinates of where the geoimage should land to be georeferenced. This wouldn't make sense in the context of the image not being an actual map image. I only want to overlay an image somewhere on the map (which works perfectly, with the rotation and everything), and afterwards extract the map coordinates of the corners of the image as it was overlaid on the map.

If imageScale =[1,1] would have meant that 1 pixel =  1 map units, I would have managed to do the math myself. Math is not the problem. The problem is that scale seems to create an image that is about 1 pixel ~ 0.625 meters, for the EPSG:3857 CRS. I don't know whether this 0.625 or a close value is hardcoded somewhere in ol-ext, or in openlayers or maybe is calculated based on something. Either way, what I am missing is a reliable way to calculate the map coordinates of the corners of the image, once the image is overlaid on the map.


texts: Can we use existing planet database for Nominatim

texts:
i have imported data from the OSM to PGSQL using osm2pgsql for generating tiles,

and now i need to add Nominatim support also, this Nominatim requires the same OSM data. can we use the existing database? or do i need to create new one?


texts: GoogleMaps URL string to display unique marker

texts:
This is a topic I thought as easy as pie to solve with just a Google search and which is in fact quite complicated.

When I use the following URL string to open a GoogleMaps window on a specific location, I get 2 markers (one green on the exact position of the Lat/Lng coordinates I put in the URL, and one red which is the nearest address).

I would like to get just the green marker, which means opening a new GoogleMaps window with a marker matching the Lat/Lng coordinates from the URL. This is quite simple, but I can't find out on the web if it possible or not.

This is possible with OpenStreetMap using the following URL.


texts: Can't set transparent raster alpha channel in GeoExplorer!

texts:
I am trying to get transparent parts of an 3-channel raster with NoData pixel values in GeoServer.

As suggested in here, I created alpha channel with pixels that have No Data Value (0) as a mask using GDAL
`gdalwarp -srcnodata 0 -dstalpha input.tif output.tif`
As a result I got 4-channel (RGB+alpha) where the pixels that previously were 0 now have 0 values in all four channels and those, that were > 0, in alpha channel have value 255. Also No Data Value is set to 0.

When I upload raster to GeoServer (v2.5.2) using QGIS plugin OpenGeo Explorer, instead of transparent pixels I get black pixels. What could be the problem? I have also tried without setting No Data Value without results.

---------------EDIT----------------
Geoserver version: 2.5.2
Raster style:

   <?xml version="1.0" encoding="UTF-8"?>
   <sld:StyledLayerDescriptor xmlns="http://www.opengis.net/sld" xmlns:sld="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc" xmlns:gml="http://www.opengis.net/gml" version="1.0.0">
     <sld:NamedLayer>
<sld:Name>RGB raster</sld:Name>
<sld:UserStyle>
  <sld:Name>RGB raster</sld:Name>
  <sld:FeatureTypeStyle>
    <sld:Name>name</sld:Name>
    <sld:Rule>
      <sld:Name>Single symbol</sld:Name>
      <sld:RasterSymbolizer>
        <sld:ChannelSelection>
          <sld:RedChannel>
            <sld:SourceChannelName>1</sld:SourceChannelName>
          </sld:RedChannel>
          <sld:GreenChannel>
            <sld:SourceChannelName>2</sld:SourceChannelName>
          </sld:GreenChannel>
          <sld:BlueChannel>
            <sld:SourceChannelName>3</sld:SourceChannelName>
          </sld:BlueChannel>
        </sld:ChannelSelection>
      </sld:RasterSymbolizer>
    </sld:Rule>
  </sld:FeatureTypeStyle>
</sld:UserStyle>

Here is the request:
request=GetMap&layers=jaun:ORO2&styles=&bbox=515813.19564,305943.471,516026.19564,306105.471&width=512&height=389&srs=EPSG:3059&format=application/openlayers&transparent=true


texts: Connecting individual road segments that participate in a dead end

texts:
I am trying to get a road feature class ready to do some analysis, regarding dead end roads. I have a road dataset cleaned up with continuous single segments, split at each road intersection. I have also added a field [DEADEND] noting any segment that dead ends (created through dangle topology). I have also identified any road segment that connect a dead end dangle with a non-dead end road (roads that have intersections on both ends). For example, a road that connects to multiple branching dead ends at one end and a non-dead end on the other, which would not be identified as a dangle in topology (as it has a vertex on each end, but is still part of a dead end).
If I query out the the segments ([DEADEND]="YES"), how do I merge the connected dead ends(sharing vertices) into individual continuous lines.
See picture for example.

I would like all segments labeled in green or purple into one line feature, without using the merge tool in the editor toolbar.


texts: How can I use icons on a my own maptiles klokantech docker container?

texts:
If I take this example:

https://www.mapbox.com/mapbox-gl-js/example/geojson-markers/

an change the mapsource from:

style: 'mapbox://styles/mapbox/light-v9'

to my own docker container:

var map = new mapboxgl.Map({
    container: "map",
    style: "http://localhost:8080/styles/osm-bright/style.json",
    center: [11.7357, 48.8],
    zoom: 8
});

map.on("load", function() {
    map.addLayer({
        id: "points",
        type: "symbol",
        source: {
            type: "geojson",
            data: {
                type: "FeatureCollection",
                features: [
                    {
                        type: "Feature",
                        geometry: {
                            type: "Point",
                            coordinates: [11, 49]
                        },
                        properties: {
                            title: "Mapbox DC",
                            icon: "monument"
                        }
                    },
                    {
                        type: "Feature",
                        geometry: {
                            type: "Point",
                            coordinates: [11, 50]
                        },
                        properties: {
                            title: "Mapbox SF",
                            icon: "harbor"
                        }
                    }
                ]
            }
        },
        layout: {
            "icon-image": "{icon}-15",
            "text-field": "{title}",
            "text-font": ["Open Sans Semibold", "Arial Unicode MS Bold"],
            "text-offset": [0, 0.6],
            "text-anchor": "top"
        }
    });
});

Then disapear the icons on map. How can I use icons on my own container?

I use mapbox-gl-js/v0.41.0


texts: Python find string with two decimal characters in ArcGIS field calculator?

texts:
I have two columns in a personal geodatabase table that I'm trying to add up (to a third field) using a python function in field calculator, however some of the string values are not real numbers.

e.g.

field1     field2
12.24      2
h3.3       4
xyz        1
932.33.32  3

I have a conditional statement that checks for letters, however I'm stuck on how to check for the two decimal character occurrence,

  932.33.32

For the two decimal string example I just want to return the value of field2.  Here is my code so far.

import re
def cal(field1, field2):
  #check for two decimal character occurrence
  elif re.search('[0-9]+', x):
    myvalue = float(field1) + float(field2)
    return myvalue
  else:
    myvaluefinal = float(field1) + float(field2)
    return myvaluefinal


texts:
I'm having problems with SAGA GIS. It segfaults on launch, either from the command line or when called from QGIS.

If I run saga_cmd using gdb, I get this error:-

(gdb) run
Starting program: /usr/bin/saga_cmd
[Thread debugging using libthread_db enabled]
Using host libthread_db library "/lib/i386-linux-gnu/libthread_db.so.1".
____________________________

   #####   ##   #####    ##
  ###     ###  ##       ###
   ###   # ## ##  #### # ##
    ### ##### ##    # #####
 ##### #   ##  ##### #   ##
____________________________

SAGA Version: 2.2.7

[New Thread 0xb1dbbb40 (LWP 14864)]
[New Thread 0xb15bab40 (LWP 14865)]
[New Thread 0xafdb9b40 (LWP 14866)]

Program received signal SIGSEGV, Segmentation fault.
__strlen_sse2_bsf () at ../sysdeps/i386/i686/multiarch/strlen-sse2-bsf.S:50
50  ../sysdeps/i386/i686/multiarch/strlen-sse2-bsf.S: No such file or directory.

My setup:-

Ubuntu Tahr 32bit
QGIS 2.16.3 . This is a fresh install. Using the standard QGIS install as defined here. This is the Stable release, and not UbuntuGIS.
QGIS Running against GDAL/OGR 1.11.2.
QGIS Compiled against GDAL/OGR 1.10.1
SAGA 2.2.7

I found this discussion, which suggests it's the version of GDAL i'm using. This seems to tie in with when this problem first started - when GRASS 7 support started working for me from QGIS.

Is there a way around this? If GRASS and SAGA need different versions of GDAL, I'm worried that attempting to fixing this for SAGA might break GRASS 7 support, which I also need.


texts: One liner to change table style in mapbasic

texts:
Is there a one liner using Set Style to change all the symbology of one tab file?

I want it to be
Global Pen (1,2,0)   Global Brush (1,16777215,16777215)

There's https://groups.google.com/forum/#!topic/mapinfo-l/BZM5ZXfdXSY and http://www.mapbasichelp.com/2014/05/change-object-style-with-mapbasic-alter.html but these seem more complex than necessary.


texts: OpenLayers GeoJSON t.getState is not a function

texts:
I'm trying to make my own example with Openlayer. I have to load a GeoJSON on my map. But it doesn't work and I have a message on the console where is written : Uncaught TypeError: t.getState is not a function

   var geojsonObject = {
     'type': 'FeatureCollection',
     'crs': {
       'type': 'name',
       'properties': {
         'name': 'EPSG:3857'
       }
     },
     'features': [{
     'type': 'Feature',
     'geometry': {
       'type': 'Point',
       'coordinates': [0, 0]
       }
     }]
   };

 var masource = new ol.layer.Vector({
     features: (new ol.format.GeoJSON()).readFeatures(geojsonObject)
   });

 var monlayer = new ol.layer.Vector({
     source: masource
   });

   var map = new ol.Map({
     target: 'map',
     layers: [
       new ol.layer.Tile({
        source: new ol.source.OSM()
       }),monlayer
     ],
     view: new ol.View({
       center: [0, 0],
       zoom: 2
     })
   });

I don't understand what is going wrong.


texts: Using QGIS generated web map in Android app

texts:
I am on Windows 10 (x64 machine). Using QGIS' (version 3.12) qgis2web plugin, I have created a web map of my location of interest. The map contains OSM basemap with overlayed points shapefile (classified into 5 categories). There are more than 2 million point locations.
In the web app, there are index.html and JS and CSS source files in a directory. The app runs locally and displays the map fine. Now, I would like to show the same web app in my Android application (Android Studio, Java).
I thought of using Android's WebView class, or is there any other way to do the same?
I am also concerned about my web app as it is very slow while zooming-in and out- of the map.


texts: How to get attribute of nearest neighbour of polygons from a points layer

texts:
I have 2 layers, one of polygons and one with points. Both have some attributes.
For each polygon, I would like to take one attribute of points layer from the nearest point and put it in a field of polygons layer.
Is there a way to do this with QGIS 3?


texts: R Error in sp::CRS(SRS_string = from$wkt): unused argument (SRS_string = from$wkt)

texts:
A previous script that I had is no longer working after the sf package major change. I am using sf 0.9.6 and raster 3.0.2. I have tried to modify my code but cannot run raster::extract().
My previous code now breaks on the second command and does not generate extent_program. I receive: Error in sp::CRS(SRS_string = from$wkt): unused argument (SRS_string = from$wkt).
# Read in polygons from 'my_shapefile.shp'

sf_program_fallow <- 'my_shapefile.shp' %>% sf::st_read() %>% sf::st_transform(crs = 32611)

# Generate extent of polygons (now generates error)
extent_program <- extent(sf::as_Spatial(sf_program_fallow))

# Calculate zonal statistics (note: 'my_ndvi_RasterLayer' created earlier)
ndvi_fallowed_fields <- raster::extract(x = 'my_ndvi_RasterLayer',
                       y = sf_program_fallow,
                       method = 'simple',
                       fun = mean,
                       na.rm = T,
                       sp = TRUE)

Updating the first line of code to the following produces extent_program but not ndvi_fallowed_fields. I get the same error: Error in sp::CRS(SRS_string = from$wkt): unused argument (SRS_string = from$wkt)
sf_program_fallow <- 'my_shapefile.shp' %>% sf::st_read() %>% sf::st_transform(., sf::st_crs(32611))


texts: Using historical imagery Google Maps API for Business

texts: Is there a way using Google Maps API for Business to see historical aerial imagery, like you can in Google Earth Pro, with the slider bar they have at the top of GEP?


texts: Bent raster after Build Overviews for mosaic dataset using ArcMap?

texts:
I created a mosaic dataset from a set of 12 tiff's spanning the united states.  Before running Build Overviews, the images all join together nicely.  But after, they look like this:

Any idea what is going on?  Did i do something wrong?  Should I have merged the mosaic to a single image first?

Per the comment by @PolyGeo, here is a detailed step by step:

Gather source data from the FAA (The 12 files under High Altitude Enroute (U.S.) GEO-TIFF Files)

Unzip each zip file into a directory
Each of these tif files have an associated html file that shows some metadata, specifically:

Geodetic_Model:
    Horizontal_Datum_Name: North American Datum of 1983
    Ellipsoid_Name: Geodetic Reference System 80
    Semi-major_Axis: 6378137.000000
    Denominator_of_Flattening_Ratio: 298.257222

Open ArcMap (10.3.1.4959 License type Advanced), then create a new Blank Map
Under catalog, right click on any gdb, choose New -> Mosaic Dataset...
On the Create Mosaic Dataset dialog:

Under Coordinate System, choose NAD 1983
Click Environments...
Under Raster Storage set Pyramid resampling technique and Resample to CUBIC
Click OK

Once the dataset is created, right click it and Add Rasters...

Under Input Data choose Dataset then browse for the 12 tif files from 1.1 above
Set Environments... the same as 4.2 above
Click OK

Then, I created feature polygons for each tif by clicking through the map file at all the corners of the useful data and used that to modify the footprints for each raster (note that I also set Always Clip The Image To Its Footprint)

At this point I have a nicely stitched together mosaic dataset.  But when i zoom out past 1:1,250,000 it stopped rendering the image.  I analyzed the dataset which provided the warning that I was missing overviews.  So i went ahead and right clicked the mosaic dataset, then Optimize -> Build Overviews.  The result, is what you see in the attached image...


texts: Why is it important to use feature layers in ModelBuilder?

texts: I am trying to understand creation of a geoprocess with ModelBuilder, but I don't know why it's important to use feature layers instead of feature classes when creating a geoprocess with ModelBuilder. Can someone please explain why?


texts: Point to every other point as dataframe in R

texts:
I understand there are several tools out there that can help me accomplish my goal, but I think I'm running into more of a conceptual roadblock. I've tried to implement a fully reproducible example in seeking assistance.

The situation: I have two point shapefiles, one with 5 points (fewpoints.shp) and the other with many points ~100 (manypoints.shp). For each of the many points, I'd like to find the distance to each of the 5 points, such that the final df looks like:

ID     DistPnt1    DistPnt2    DistPnt3    DistPnt4     DistPnt5
1
2
3
...
100

I've been able to generate a distance matrix and melt it into a dataframe, but I'm having trouble understanding how I'm supposed to retain the original ids such that I can merge them back onto manypoints.shp.

Reproducible example:

library(USAboundaries)

MN <- us_states(resolution="low", states = "MN")

library(sf)

MNspdf <- as_Spatial(MN)

smallpoints <-spsample(MNspdf , 5, type = "random")

manypoints <-spsample(MNspdf , 100, type = "random")

spDists(manypoints, smallpoints)

Is it safe to assume that each row corresponds to the original rows of the manypoints.shp? Am I missing something?


texts:
I have a set a points, and I want to create 10 clusters.

I used attribute based clustering plugin, which worked well for creating the clusters, but the points count isn't perfect, some clusters have two times more points than others.The points are close to each other, but clustering groups the closest points without taking into account the total number of points.

So I was wondering if QGIS has a feature to cluster points and make points count not necessarly equal but at least don't have a big differences.

Hierarchical clustering seems to be an option, but i can't find any implementation of it in QGIS.


texts:
I upgraded to Mint 16 and QGIS 2.4 and unable to get ecw to work following various guides available.

How and why is it so hard?

I run it on a windows machine out of the box.


texts: Tool in QGIS to construct polygons based on network of lines

texts: Is there a tool in QGIS similar to "Contruct polygons" in ArcGIS that creates polygons from existing lines?  This is different from convert lines to polygons that requires closed loops (minimum of three vertices). Instead, I am interested in creating lot polygons from lot boundaries or creating city blocks based on streets.


texts:
Is there a way to add a 'dummy' label to symbols in the legend such as this?

The layer it self has a brown/red outline, as shown in the image, and labels are turned on. Since there are a ton of different numbers and letters for each shape in the layer, I want to show a simple outline with a 'dummy' label.

I've tried the method suggested by csk, but I'm aiming at getting only one symbol - not one for the symbol and one for the label.


texts:
I wrote a python script to configure Weblogic Domain. I can see there is an error when im trying to execute with the below lines. I tried debugging with different options and still it is throwing the same error. Please help me out here.
AttributeError: 'NoneType' object has no attribute 'split'
Exporting the Properties to variables..
Problem invoking WLST - Traceback (innermost last):
File "/test/wls_domain_creation.py", line 304, in ?
File "/test/wls/wls_domain_creation.py", line 61, in export_properties
Here are the lines 61 and 304
machines = _dict.get("machines").split(',')
export_properties()
Thanks in advance


texts: Arc collector: cannot enable offline mode for map

texts:
I am trying to setup arc collector for offline data collection in remote areas. I have followed steps from a few ESRI tutorials.

Problem: when I try to "enable offline mode" for the web map within Arc Online I get a pop-up message that says, "This web map cannot be enabled for offline uses..." (My Content > web map details > edit).

Background: In ArcMap, I created a file geodatabase with domains. I then exported separate polygon, point, and line features to the gdb (features had been previously created). I published the map as a service, changed the capabilities to Feature Access and checked "Create, Update, Sync, and Delete". The map published successfully with a couple "Medium" warnings but zero errors.

Within Arc Online, the features and the web map both say sync is enabled but, again, I am not able to enable offline mode, which is the difference between using collector or not for data collection (given our remote surveys).

What else could be preventing me from enabling offline mode?


texts: How can I get temperature data for each Countries?

texts:
I need maximum and minimum temperature data of each countries.

I searched a lot but could not find the values.

Could you give links if it is possible?


texts:
How do I use ORDER BY, GROUP BY & LIMIT in the following query through column name KPH?

SELECT *
FROM data
WHERE  ST_Within(
    ST_Transform(geom, 2163),
    ST_Buffer(
        ST_Transform(
            ST_GeomFromText('LINESTRING (77.2089421749115 28.5225030428734, 77.2088992595673 28.5223333626164,77.1612417697906 28.5640959203597)', 4326),
            2163),
        500)
    );


texts: How to install PostGIS shp2pgsql-gui on CentOS 6?

texts:
I am installing PostGIS 1.5.3 on a CentOS 6.2 virtual machine to get started web mapping with Geoserver.I have successfully installed PostgreSQL 9.1, PostGIS 1.5.3 and PgAdminIII on the system.  However, I am noticing that there are no PostGIS GUI's installed.

I was under the impression the PostGIS includes the Shapefile and DBF loader plugin for PgAdminIII even in Linux.  When I navigated to /usr/bin shp2pgsql was there but not shp2pgsl-gui.

Is there anyway to install the GUI for PostGIS in CentOS 6.2?  I am very new to Linux and the OSGeo Stack.  I have installed the above software using YUM and the rhel-6 PostgreSQL repository.

Thanks in advance for your help.


texts:
It is tricky to manage the extent of a coordinate system when that extent is not an aligned rectangle in another coordinate system.

For example consider Thursday Island, located at:

142.2°, -10.2° (WGS84, EPSG:4326) or equivalently,
1146km, -1156.5km (Australian Albers projection, EPSG:3577).

According to SpatialReference.org (similar values to epsg.io) this latter projection has:

WGS84 Bounds: 108.0000, -45.0000, 155.0000, -10.0000
Projected Bounds: -2690013.3995, -5098960.4467, 2579169.8548, -1281018.4757

Notice that Thursday Island is well inside of the lat/lon bounds, however it (along with sections of Australian mainland at both Cape York Peninsula and Kakadu) falls outside of the projected bounds.

Say I want to create a web map application using OpenLayers, I need to know the maximum extent of the map in the projected coordinates. The problem is that if I source this information from a standard online reference, I will mistakenly cut off parts of Australia even when choosing a projection that was specifically designed to cover Australia.

Presumably the root cause is that the extent has been represented by an aligned bounding-box and then, rather than transforming the entire boundary of this region precisely (before generating a new bounding box), only two vertices of the box have been properly projected (the bottom-left and top-right corners). Such a method is invalid due to curvature/nonlinearity.

How is this problem normally handled? Is there any authorative source for the extent of a map projection? In particular, is Proj.4 able to report the maximum domain of validity for a transformation? I'd like a solution that works in all cases, including those where the region of validity does not resemble a rectangular coordinate patch (e.g. rHEALPix and QSC).


texts: Getting std::bad_alloc save as NE_100m_land to EPSG3995 for reprojection from QGIS?

texts:
I have loaded ne_110m_land.shp into a new project in QGIS 3.4 running on Mac. I followed what seems to be the sequence in every tutorial to reproject this map. In this case I want to project the map to an Arctic Sterographic projection.

I have done the following (and some more options before I started noting tham)

Right click ne_110m_land
Click export
Click 'save features as' - This brings up the correct dialogue box.
Format: I have tried Geopackage and ESRI shapefile
File name: I give the file name MI2 and folder for it to go to. (I have tried both an empty, new folder and putting this file in with the other ne_110m files)
CRS: I select EPSG3995 WGS84 Arctic Polar Stereographic
Add saved file to map is ticked
Click OK.

I get the dialogue box std::bad_alloc
The layer appears to be added, but the map is the original, except glitchy.

If I do not add saved file to map there is no glitchiness. If I close the current project and open the saved file I also get the dialogue box std::bad_alloc. If I reopen the ne_110m_land file after trying to open MI2, it is glitchy in the same way. I must restart QGIS to get the file to display properly again.

I tried downloading 3.8 and the same happened.

I have tried with version 2 and version 4 of the ne_110m_land files and the same happened.

I tried different coding (alternatives to UTF-8) with no difference

I tried Stereographic World and it seems to have worked as it should in the tutorials.

Triton North Pole Stereographic resulted in an openable file, which does center on the North pole, but it is very strange.

I finally was able to get EPSG 102035 North Pole Orthographic to work.

Is the issue with EPSG3995 WGS84 Arctic Polar Stereographic a bug?

Or is this somehow an inappropriate choice in a way I don't understand?

Is Triton North Pole also a bug or an inappropriate option for reprojection in a way I haven't understood?


texts:
I'm trying to test the code example on the page http://docs.geotools.org/stable/userguide/tutorial/geometry/geometrycrs.html
I copied this example into eclipse

 Transaction transaction = new DefaultTransaction("Reproject");
    try (FeatureWriter<SimpleFeatureType, SimpleFeature> writer =
                    dataStore.getFeatureWriterAppend(createdName, transaction);
            SimpleFeatureIterator iterator = featureCollection.features()) {
        while (iterator.hasNext()) {
            // copy the contents of each feature and transform the geometry
            SimpleFeature feature = iterator.next();
            SimpleFeature copy = writer.next();
            copy.setAttributes(feature.getAttributes());

            Geometry geometry = (Geometry) feature.getDefaultGeometry();
            Geometry geometry2 = JTS.transform(geometry, transform);

            copy.setDefaultGeometry(geometry2);
            writer.write();
        }
        transaction.commit();
        JOptionPane.showMessageDialog(null, "Export to shapefile complete");
    } catch (Exception problem) {
        problem.printStackTrace();
        transaction.rollback();
        JOptionPane.showMessageDialog(null, "Export to shapefile failed");
    } finally {
        transaction.close();
    }
}

But I got some error with transform function. Eclipse "ask" to change Geometry type to Envelope. But when I changed it didn't solve this problem . Eclipse asked to change envelope to geometry again. How can I fix this error?


texts: Strange ArcSDE error inserting null into state_locks?

texts:
I'm working with a heavily customized version of ArcMap 9.3.1. It's heavily customized, but it mostly uses class extensions and special validations.

Environment:

ArcSDE installed in: linux, red-hat, version 9.3.1, no service packs;
Oracle 10g;
Versioned database, moving edits to base, single DEFAULT version;
ArcMap connecting via direct connection;

I have an user, creating simple points using a custom tool we made. This tool basically fires other actions, but nothing too deep. The user, after a while, creating points (simple points!) get's a weird database message, specifying:

    ORA-01400

    Cannot insert null into sde.state_locks.sde_id

    [sde.default][state_id]=6935

After this any point he creates fails, with the same message.

Anyone has a clue of what is happening? Has anyone ever see this?

Just to be clear, the error is intermitent. The tool we made works just fine, with unit tests and everything.


texts: Getting accurate measurements in EPSG: 27700 using QGIS 3.10

texts:
When taking measurements in a project where all layers are EPSG:27700 (with elipsoid Airy 1830) the measurements produced are incorrect.
Is there a better CRS for distance measurements in the UK? Or is there something I can check in the project set up?


texts: Ortho imagery overlapping a top of Hillshade

texts:
I have a hillshade mosaic dataset made from DEM and I also have another mosaic dataset made with a aerial photo.
I will like to overlapping so I can see both layer together, however I don't want to transparency aerial photo mosaic because I lose quality of aerial photo.

I tried https://www.esri.com/arcgis-blog/products/product/imagery/an-alternative-to-overlaying-layer-tints-on-hillshades/

I also tried http://desktop.arcgis.com/en/arcmap/10.3/manage-data/raster-and-images/wkflw-combining-hillshaded-dem-with-a-scanned-map.htm

Both failed. Is there a any suggestion about it? Anyone tried it?


texts:
I am trying to work with Linear Reference system as I found It in ArcGis website: http://resources.arcgis.com/en/help/main/10.1/index.html#//003m00000002000000
But I'm losing something in the backgrounds, maybe in the initial datasets preparing. The metric calculations results shows nothing in the table field.

I am doing this task using the road shape from Open Street Maps shapefile, also I got some Garmin GPS coordinates obtained from specifics and aleatory points in some of these roads in my selected area(two to three points at least per road).

These coordinates were imported to a shapefile and my intention is to use these coordinates as a real reference points to create routes with a linear metric system using ARCGIS tools to create, calibrate(to extrapolate after and before and interpolate between points) to help me locate a thousands of events (that is specified in kilometers in the event table) along those roads instead I have to go in any event place to pick GPS coordinates for each event.

I pasted some screens of the initial datasets I got so far. I pointed the M_field to "shape*" field in calibrate process, is It correctly?

The fisrst screen: The original table from OpenStreet;
The second one: after routes created;
The third: the coordinates table.


texts: Chaining Rendering Transformations: Making isotherms using Barnes Surface + Contour process?

texts:
I'm trying to calculate isotherms following the istructions showed in the post on OSGeo blog (http://blog.opengeo.org/2013/01/10/chaining-rendering-transformations-in-geoserver/)

I'm using GeoServer 2.2.3 on a Windows 7 machine (8 Gbyte RAM, biprocessor)

I've a shapefile with points and temperature  and I'm trying to use the Meteo-1.sld

Using the layer in preview

http://:8080/geoserver/CesareWorkspace/wms?service=WMS&version=1.1.0&request=GetMap&layers=CesareWorkspace:CentralineTemperature-01-01-2001&styles=&bbox=318287.0,4880823.0,515945.0,5142603.0&width=386&height=512&srs=EPSG:32632&format=application/openlayers

it works fine but when I try to use le Meteo-1.sld

http://:8080/geoserver/CesareWorkspace/wms?service=WMS&version=1.1.0&request=GetMap&layers=CesareWorkspace:CentralineTemperature-01-01-2001&styles=&bbox=318287.0,4880823.0,515945.0,5142603.0&width=386&height=512&srs=EPSG:32632&format=application/openlayers&STYLES=Meteo-1

the map is white. It seems that there are no errors in the log.

When I try to use the  Barnes Surface WPS directy (gs.BarnesSurface),  I obtain a black tiff ....

At the following link https://www.dropbox.com/sh/fmbdhzwp5mobpcq/lf43d3SykM there are:

data.zip --> the shapefile I'm using
Meteo-1.sld --> the sld I'm  using
geoserver.log --> geoserver log

Could someone tell me which configurations I've to use in the sld  for my data? Or give me a sample data and a sld that works?


texts: Scanning a raster by moving windows. Output rasters off by a cell. ArcGIS, python

texts:
I'm working in extracting smaller rasters from a bigger one (yellow). I'm using only python, some numpy maybe. I succesfully extracted a patch from a big raster, secuentially, by moving windows. So I ended up with several small rasters, 'quilted'/patching the bigger one. The code works for the first column of patches(red colored) (See graph, sample point 'A'), where the extracted raster and the source raster has the same value at any given cell (extracted raster in red means they match perfectly the source raster). However, all the patches in gray are off by one cell! See the values at sample points 'B' and 'C'. The right value in the patched raster is in the cell besides. And this displacement is kept to all the patches, except the first column. Odd that some rasters are ok and some not. What is wrong with my code?

for rows in range (0,8):
   for columns in range (0,10):
     try:
       small_raster_lower_left  = arcpy.Point(large_raster_extent.XMin + (columns * window_size * large_raster_cellsize_X),
                                              large_raster_extent.YMin + (rows    * window_size * large_raster_cellsize_Y)
                                              )

       numpy_array = arcpy.RasterToNumPyArray(large_raster,
                                              small_raster_lower_left,
                                              window_size,
                                              window_size,
                                              nodata_to_value = -99999
                                              )

       small_raster = arcpy.NumPyArrayToRaster(numpy_array,
                                               small_raster_lower_left,
                                               large_raster_cellsize_X,
                                               large_raster_cellsize_Y,
                                               value_to_nodata = -99999
                                               )

       arcpy.DefineProjection_management(small_raster, large_raster_spatial_ref)
       raster_counter += 1
       small_raster_name = "small_dem_" + str(raster_counter) +"_c" + str(columns) + "_r" + str(rows)
       small_raster.save(small_raster_name)
       del numpy_array
       del small_raster

  except:
       del numpy_array
       del small_raster


texts: raster not displaying below certain scale

texts:
i've clipped a raster in 20 smaller images and i can zoom in until 1:150 in the original raster but, in the smaller ones, when i go below 1:4000, the workspace goes blank.

i'm working with arcgis 9.3 and can't find any configuration/solution to this problem.


texts: Making chart of FIRMS data in Google Earth Engine

texts:
I'm trying to analyze crops burned areas in GEE.
I want to use FIRMS data to make a chart of "T21" and "confidence" bands to see how these variables change in a specific period (a couple of months).
I have tried two different ways, but none seem to work. I get "Error generating chart: No properties to chart"
The two scripts I have tried are
1-
var polygon = shpfile.filterMetadata("Polygon","contains","3187956");
var band = ee.ImageCollection('FIRMS').select('T21')
.filterDate('2020-09-01','2021-01-12');
// Kelvin to celcious
var tempToDegrees = band.map(function(image){
var props = image.toDictionary(image.propertyNames());
return ee.Image(image.multiply(0.02).subtract(273.15).setMulti(props));});
// Add timeband
var withTimeBand = tempToDegrees.map(function(image){
var dateMillis = ee.Image(ee.Number(image.get('system:time_start'))).toInt64();
return image.addBands(dateMillis.rename('dateMillis'));});
print('Temp_chart',ui.Chart.image.series(tempToDegrees, polygon, ee.Reducer.first(), 500));

2- Other script, maybe column chart with mean daily data
var iniDate = ee.Date('2020-09-01')
var endDate = ee.Date(ee.List(FIRMS.get('date_range')).get(1)) //last date is currently 2020-05-17
print('lastDay available', ee.Date(ee.List(FIRMS.get('date_range')).get(1)))

//Add ImageCollection
var filterFIRMS = FIRMS.filterDate(iniDate, endDate).select(['T21']);

//Then initiate the date list with a day less than the availability,
//as the last day available does not contain images:

var difdate = endDate.advance(-1, 'day').difference(iniDate, 'day'); //calculate steps/days in between
var createList = ee.List.sequence(0, difdate) //make list with the needed number of entries
var listdates = createList.map(function(day){ //change list into list of dates
return iniDate.advance(day, 'day')})
print(listdates)

//daily data
//Then filter the image collection using the filterDate, which is (I think) easier than filtering
//on year and doy consecutively:
var firmsImageCollection  = ee.ImageCollection.fromImages(listdates.map(function(summarize_day){
var filterCol = filterFIRMS.filterDate(ee.Date(summarize_day), ee.Date(summarize_day).advance(1, 'day'))
return filterCol.sum().select(['T21'], ['dailyTemp'])
.copyProperties(filterCol.first()).setMulti({
Date: ee.Date(summarize_day), 'system:time_start': ee.Date(summarize_day).millis()});}))
print(firmsImageCollection)
print(ui.Chart.image.series(firmsImageCollection, selecpol, ee.Reducer.mean(), 5000)
.setChartType('ColumnChart'));

I think a lineChart of T21 and a daily mean column chart of confidence band would make me happy.


texts:
I have a script that creates a number of features "in_memory." Once I have them in memory I need my script to create a list of everything "in_memory" - I'm not sure how to achieve this. Is anyone knowledgeable on how to achieve this result?

Below is the portion of my code that creates my features in memory.

for fc in fcs:
    arcpy.Clip_analysis(fc, Buffer, "in_memory" + "\\" + fc)


texts: Issue with overlapping features in PostgreSQL database st_difference not working

texts:
I am working on a problem in a postgreSQL database with the Postgis extension.  My issue is that I have a bunch of polygons in a table and they have some overlap.  These are US census blocks.  Not sure why this is.  So I have been reading what to do and getting some issue.  For the application to work everything needs to be multi_polygons in geography format.

So I need to get this overlap to go away and have the two census blocks butt up to each other.  I have been reading some similar issues and cannot get what they did to work when using st_difference.  Can you compare polygons in the same table with st_difference to remove the overlap area?  I need one of the polygons to take it and the other to not.  Similar to clip in ArcMap/Qgis.  But I cannot use those programs and have to do it all in the database.  If not what would be a good way to do this in Postgis.

Here is what I have been trying and not succeeding with.

update public.census_blocks
set shape_trimmed = st_multi(st_difference(shape_trimmed::geometry,
shape_trimmed::geometry))::geography
where st_intersects(shape_trimmed, shape_trimmed) and
census_blocks.service_area_id = 4017;


texts:
I have an internal map service that serves up GIS records via JSON objects.

For the purpose of this question, instead of using my internal/private
data, we can use the JSON data on ESRI's public Sample Server 6.

https://sampleserver6.arcgisonline.com/arcgis/rest/services/LocalGovernment/Events/MapServer/0/query?where=1%3D1&outFields=*&f=pjson

Question:
I want to select the JSON data and insert it into an Excel table.
In other words, I want to take the JSON data:

And put it into a tab in an Excel spreadsheet as rows:

How can this be done in Excel 2016?


texts: Designing the Layout in opendata.ArcGIS.com

texts:
I am trying to build an open data portal through http://opendata.arcgis.com/, its all the rage right now and pretty exciting development from Esri.

It is fairly new, but it has limited documentation and community support. I would like to completely restyle the site, similar to how http://geo.wa.gov/ or http://data.imap.maryland.gov/ did their sites.

opendata.arcgis allows you to edit the design using their wizard and markdown, or you can customize html (only within the div class="custom-layout" ), the css, and one could insert script tags to extend the html and behaivior:

To create one of the afore mentioned sites it would take a serious amount of hacking inside their little code editor.

My Question: Is there a better way of doing this? Is there a way to download and modify all the code? or do we have to do it all inside this (very) little code editor, using javascript to modify html?

(my rep isn't high enough to create new arcgis-opendata tag or answer my own question)

EDITED The solution might be to delete everything in the Header HTML of the "Design Layout" tab and use an iframe with my own css:


texts: Use Get Field Value tool to calculate field?

texts:
I am running a summary statistic on a clip to get acre counts by area type. I want to then add a second field that has the acre count from the total area. A third field be added so that it shows as the percentage of each acre type as compared to the whole area.

First field = Summary Statistics of area type by acre count

Second field = Total acres of project area <--- This is where I am having trouble

Third field =  (([first field]/[second field]) * 100)

The problem is, I am trying to get the value from the initial project area using the Get Field Value tool. I am able to do this but when using the Calculate Field tool to calculate the third field, I do not know how to get this option.

Can I set the value from the Get Field Value to a variable then in the Field Calculator use Python to call the variable? Will it recognize the variable?


texts: Serving multiple styles from tile server

texts: I installed an OSM Mapnik tile server (according to the switch2osm documentation). Now, i'd like to serve tiles in different styles from the same server. Is this somehow possible?


texts: Can I use VB to search against an ArcServer service in ArcMap?

texts:
I've published an ArcServer service that is going to be consumed by out of the box, standard ArcGIS on the desktop.

Is it possible to write custom tools in VB to query layers from that service in the ArcMap interface?

I'm using ArcGIS 9.3.


texts: label orientation after dxf to shp conversion

texts:
I'm writing a script with python in QGIS that allows me to convert a dxf file to shp. In the dxf file I have oriented text but after the conversion I lose this information. Is it possible to keep it?
this is my script :

import os
from PyQt4.QtCore import *
import datetime
print datetime.datetime.now().minute
uri="C:/Users/Lenovo/Desktop/stage_brome/restitution.dxf|layername=entities"
vlayer = QgsVectorLayer(uri, 'layer', 'ogr')
##For each sublayer, diferent type of geometry, load a layer to map.
subLayers = vlayer.dataProvider().subLayers()
types=["LineString","Point","Polygon"]
for sublayer in subLayers:
    #Extract the geometry type
    geom_type = sublayer.split(':')[-1]
    #Set the path
    url = str(uri)+"|layername=entities|geometrytype="+str(geom_type)
    #Name for sub layer
    dxf_file_name = os.path.splitext(os.path.basename(url))[0]
    layer_name =str(dxf_file_name)+" - "+str(geom_type)
    #Create layer
    sub_vlayer = QgsVectorLayer(url, layer_name, 'ogr')
    #Add layer to map
    QgsMapLayerRegistry.instance().addMapLayer(sub_vlayer)
    if geom_type in types:
        #list of string that contains the name of each layer
        layer_types=[]
        #list of vector layers
        vlayers_by_type=[]
        #list of providers for each vector layer
        vlayers_providers=[]
        geom_types=[]

        for feature in sub_vlayer.getFeatures():
            #fill the lists (vector layer and it's provider and attributes of each layer)
            if (feature.attributes())[0] not in layer_types:
                name=(feature.attributes())[0]
                layer_types.append(name)
                # create vlayer
                vl=QgsVectorLayer(geom_type, "temporary_points", "memory")
                geom_types.append(geom_type)
                pr = vl.dataProvider()
                # add fields
                fields=[field for field in sub_vlayer.pendingFields()]
                pr.addAttributes(fields)
                # tell the vector layer to fetch changes from the provider
                vl.updateFields()
                vlayers_by_type.append(vl)
                vlayers_providers.append(pr)
                error = QgsVectorFileWriter.writeAsVectorFormat(vl, "C:/Users/Lenovo/Desktop/stage_brome/output/layers/"+name+"-"+str(geom_type)+".shp", "CP1250", None, "ESRI Shapefile")
#                if error != QgsVectorFileWriter.NoError:
#                    print "failed  "------------------------------------

        for feature in sub_vlayer.getFeatures():
            name=(feature.attributes())[0]
            if name in layer_types:
                # add features
                vl=vlayers_by_type[layer_types.index(name)]
                pr=vlayers_providers[layer_types.index(name)]
                pr.addFeatures([feature])
                # update layer's extent when new features have been added
                # because change of extent in provider is not propagated to the layer
                geom_type=geom_types[layer_types.index(name)]
                vl.updateExtents()
                error = QgsVectorFileWriter.writeAsVectorFormat(vl, "C:/Users/Lenovo/Desktop/stage_brome/output/layers/"+str(name)+"-"+str(geom_type)+".shp", "CP1250", None, "ESRI Shapefile")
                if error != QgsVectorFileWriter.NoError:
                    print "failed  "+str(layer_types[vlayers_by_type.index(vl)])

print datetime.datetime.now().minute


texts:
I'm trying to use the region grow raster function on a DEM in ArcGIS pro. The goals is to use a single seed point to identify all the surrounding cells that are within a certain elevational range (i.e. +/- 1 5 meters) of that seed.
I created a feature class with my seed point and with fields for the Maximum Growth Radius and Similarity Threshold (following the documentation here) parameters. Regardless of how I set the values in those fields, however, the output raster is always identical to the input - the region around my seed point never gets any bigger! If I inspect the values in the surrounding DEM cells, however, I can see that they are well within the similarity threshold.
Any tips on how to get this tool to work? Also open to using alternative tools (qgis, saga, gdal, etc) that may be able to accomplish the same goal.


texts:
I have two multi-polygon layers, one a map of land uses, the other borders of administrative areas.

I am interested in analyzing the amount of land use in each administrative area, but the polygons in the map cross the boundaries, therefore it need to split these polygons into parts along the administrative boundaries.

I know I can use the Split Tool but that makes each administrative area a separate layer. Is there a way of splitting up the polygons in the map layer but leaving the whole map as a single layer?


texts:
I'm new in PyQgis.
I can create vector layer and lines/polygons on it :

layer = QgsVectorLayer("LineString", "points", "memory")
pr = layer.dataProvider()
pt = QgsFeature()
point1 = QgsPoint(0,0)
point2 = QgsPoint(100,100)
pt.setGeometry(QgsGeometry.fromPolyline([point1,point2]))
pr.addFeatures([pt])
layer.updateExtents()

Is there a simple way, how can I visuall see that line in window?

I tried canvas, but I wasn't able to display anything on it.


texts: Calculate 'Rectangle' Coordinates Given 2 Points and width

texts:
I want to create a rectangular polygon using two points as guides.
So let's say a journey starts in Egypt and ends in London, my polygon should have 4 points:

10 miles further from London than Egypt is, following the line between them (roughly south in this example).
Halfway between the two cities but 50 miles at right angles from the line that joins them.
Like the above point but 50 miles in the other direction.
20 miles further from Egypt than London is, following the line between them (roughly north in this example).

I'll end up with a rough diamond shaped polygon that would completely contain the straight line journey from Egypt to London.
I hope this makes sense; any help for how I can calculate the 4 points is appreciated.
Rough Map Diagram


texts: What is the best way to programmatically convert between WKT and Proj4 string?

texts:
Some shapefiles have a .prj file associated with it, and the .prj file contains the projection info of the shapefile in the format of WKT.
Sometimes I need to convert WKT to proj4 string, and sometimes I need to convert it back.

Is there any ready-made library to do this?


texts: Why has an orthographic projection been used in this NatGeo map?

texts:
I'm wondering why this maps uses an orthographic projection, resulting in 3 maps to view the whole world. What were the requirements or the reasons to choose this projection (why not Mercator ) ?

this map maps the plastic distribution in the oceans


texts: Is there a place where I could download OpenStreetMaps tile packs for the united states? I would like to have offline access for the whole country.


texts:
I'm trying to export a Virtual raster (VRT) to GeoTIFF in Python 2.7 (Windows 10 System), but code runs with some problems.

When i try export the following image (.vrt file) to GeoTIFF i had a problem that i don't know how  to resolve

and this is the image when i export the vrt file

The code used is the following:

from osgeo import gdal
import os, glob, osr, numpy as np
path = 'D:\Tesis\Imagenes\L8 OLI TIRS'
VRT = "C:/DUMP_PYTHON"
pathout = "C:\Users\pp\Documents\LANDSAT\NUEVAS IMGS LANDSAT\Input ETa"

def save_raster ( output_name, dataset, raster_data, driver ,NaN_Value):
    """
    A function to save a 1-band raster using GDAL to the file indicated
    by ``output_name``. It requires a GDAL-accesible dataset to collect
    the projection and geotransform.
    """
    # Open the reference dataset
    g = ( dataset )
    # Get the Geotransform vector
    if raster_data is False:
        raster_data=g.ReadAsArray()
    geo_transform = g.GetGeoTransform ()
    x_size = g.RasterXSize # Raster xsize
    y_size = g.RasterYSize # Raster ysize
    srs = g.GetProjectionRef () # Projection
    NaN_rast = NaN_Value
    # raster_data[raster_data == NaN_rast] = 'NaN'
    raster_data[raster_data == NaN_rast] = np.NaN
    # Need a driver object. By default, we use GeoTIFF
    driver = gdal.GetDriverByName ( driver )
    dataset_out = driver.Create ( output_name, x_size, y_size, 1, \
            gdal.GDT_Float32 )
    dataset_out.SetGeoTransform ( geo_transform )
    dataset_out.SetProjection ( srs )
    dataset_out.GetRasterBand ( 1 ).WriteArray ( \
            raster_data.astype(np.float32) )

#READ DATA
i=672
files= glob.glob(os.path.join(path,'*.tif'))
albedo = gdal.Open(VRT+"/dummy_alb"+str(i)+'.vrt')
save_raster(pathout+'/Albedo_L8_'+files[i][39:46]+'.tif',albedo1,False,"GTiff",99999)

i also tried with the following code, but the result is the same:

i=672
albedo = gdal.Open(VRT+"/dummy_alb"+str(i)+'.vrt')
        src_ds = gdal.Open(VRT+"/dummy_alb"+str(i)+'.vrt')
#Open output format driver, see gdal_translate --formats for list
    format = "GTiff"
    driver = gdal.GetDriverByName( format )
#Output to new format
    dst_ds = driver.CreateCopy( pathout+'/Albedo_L8.tif', src_ds, -10 )

#Properly close the datasets to flush to disk
    dst_ds = None
    src_ds = None

So.. where is the problem?


texts: problem in reading kml data in openlayers

texts:
I am trying to read an KML file (point feature) and show it on a vector layer on OSM as the base map. My code seems to have a problem that does not show the features on the map. I don't now if its the problem with reading the data from KML file or is it related to visualizing it by Openlayers. Any help would be appreciated... Thanks!

My code:

<script src="http://www.openlayers.org/api/OpenLayers.js"></script>
    <script>
      function init() {
        map = new OpenLayers.Map("map");
        var mapnik         = new OpenLayers.Layer.OSM();
        var fromProjection = new OpenLayers.Projection("EPSG:4326");   // Transform from WGS 1984
        var toProjection   = new OpenLayers.Projection("EPSG:900913"); // to Spherical Mercator Projection
        var position       = new OpenLayers.LonLat(48.69,31.345).transform( fromProjection, toProjection);
        var zoom           = 12;
        var layer = new OpenLayers.Layer.Vector("KML", {
            strategies: [new OpenLayers.Strategy.Fixed()],
            protocol: new OpenLayers.Protocol.HTTP({
            url: "points.kml",
            format: new OpenLayers.Format.KML()
            })
        });
        map.addLayers([mapnik, layer]);
        map.addControl(new OpenLayers.Control.LayerSwitcher());
        map.setCenter(position, zoom );
      }

and my kml data:

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://earth.google.com/kml/2.2">
<Document>
<name>Agencies test</name>
<Folder>
<name>BlahBlah!</name>
<Placemark id="2">
<name>Something here!</name>
<description>test</description>
<LookAt>
<longitude>49.69</longitude>
<latitude>32.345</latitude>
<altitude>0</altitude>
<range>24.63686803544318</range>
<tilt>0</tilt>
<heading>1.387289180270979e-005</heading>
</LookAt>
<Point>
<coordinates>49.69,32.345,0</coordinates>
</Point>
</Placemark>
</Folder>
</Document>
</kml>


texts:
I want to use the package tmap to plot a simple shape with a basemap (e.g. OpenStreetMap). According to the help page of tmap this should be able by using tm_layout(basemaps = c('OpenStreetMaps'))

For the dataset of the provinces of the Netherlands this should work

# load package and data

library(tmap)
data("NLD_prov")

# plot  the shape by using tmap

tm_shape(NLD_prov) +
  tm_polygons('population')  +
  tm_layout(basemaps = c('OpenStreetMaps'))

however, no basemap is showed. Also no error message is popping up.


texts:
I am trying to change programmatically using python the grid on the map composer of QGIS.

Does anyone knows the syntax for this operation?


texts: How to create mutually clipping polygons

texts:
Is it possible to create buffers (in my case, but polygons in general also) that would create a boundary at their point of contact ? (so that they don't overlap on each other).
Like this (from left to right)

I am trying to solve this question where I need to join points representing addresses to their nearby parcel without including points related to adjacent parcels or parcels on the other side of the street.
It seems to me that there are two main options: either using voronoi/thiessen diagram but based on polygons, not points, or creating many buffers (making sure they overlap with the address points outside on the street), and then making these new polygons crop each other at the medium of their intersection line (which would be the median line of the street for street facing vertices, and the property boundary line for the vertice between two adjacent parcels/plots of land).
Spatially join attributes from point to polygon side in QGIS3 with numerous contiguous polygons and numerous adjacent points


texts: Change opacity of the layer with range slider

texts:
I want to change the opacity of the layer with a range slider (input type="range") when scrolling
here is my layer code in main.js:
const UnTiled = new ol.layer.Image({
            opacity:0.9,
            title:'UnTiled',
            source: new ol.source.ImageWMS({
                ratio: 1,
                url: 'http://localhost:8080/geoserver/shomal/wms',
                params: {
                         'VERSION': '1.1.1',
                      "STYLES": '',
                      "LAYERS": 'shomal:geotiff_coverage',
                      "exceptions": 'application/vnd.ogc.se_inimage',
                }

              })

            });
map.addLayer(UnTiled)

and here is range slider input code:
<input id="contrast" type="range" value="contrast" max="0.5" min="0" step="0.01"/>

I tried below code in main.js but error is Uncaught TypeError: UnTiled.values_.opacity is not a function :
$('#contrast').on('input', function() {
                UnTiled.values_.opacity( $(this).val());
            });

UnTiled.values_.opacity gives me the value of the layer's opacity
pleas help me to change layer ("UnTiled") opacity with range slider input


texts: UK constituency boundaries that can be used with Google maps

texts:
Is there a comparatively easy way to marry-up UK election boundary data with Google maps?

I am 100% NOT looking for anyone to do the work I should be doing, but with data sets from Open Gov, Ordnance Survey and a couple of other places I have not seen any ways that that data maps on to my (very outdated) knowledge of Google maps.

I hope someone clever might already know the answer and be able to just say a "use this (data set) and then search Google for info on ????". Then I do the work, join them up and have constituency map and then have to decide what to do with it.


texts: Building graph from database structure to implement A* algorithm?

texts:
I imported an osm file into a postgresql database using osm2pgsql.

Now, I want to implement a route/path generating algorithm (routing algo) based on two locations entered by a user.

It is not really clear for me where in my database structure I should search for these two points. For example if a user enters as starting point "Airport A" and as end point "Airport B" where should I search for these two locations so I can start implementing my routing algorithm.

I need some clear information on how the database is structured and how to create a graph so that i can implement A* algorithm between the 2 points the user has enetered.


texts: Overlaping lines differently in different locations

texts:
I have two lines (two different Shapefiles) that are part of a network and cross each other multiple times.
In some parts of my map I want one of them (let say line A) overlap the other one (B), but in other location I want the opposite, B over A.

Is there a way to do so?

I've reached a workaround that works but it is not elegant at all:

I put A over B
I copy B and put it over A
I cut the part of the "new B"(like me) that I don't want to be under A

QGIS 2.14 on Windows 7


texts: Leaflet map scrolled when dragging the page on mobile

texts:
I have a Leaflet map embedded in a div with vertical scroll. The scroll works well on desktop but on firefox mobile I get a strange behavior: when I scroll the page (touching outside the map) the map get scrolled as well.

When I finish the drag, the map goes back to its right position.

I couldn't manage to reproduce the error in a simple scenario so I made a gif:

Note that there is a marker (the image does not load but this is not relevant) that moves at a different speed than the background OSM layer.

All this happens even if I set the map options with dragging:false.

Does anybody know what can be affecting the map this way? Why it reacts at all to events that take place in the parent div?

UDPDATE: Added GestureHandler as suggested and the panel it installs gets also affected by the page scroll.


texts: Multiply a raster list with a field value list in a FOR loop in ArcPy

texts:
I'm trying to create a simple map algebra calculation not managed by ModelBuilder, but writing it in Python is trickier than I thought.
I have a list of raster in a GDB and I must multiply them by each value contained in a population field.
I created my raster list, the a search cursor for the population value, and I use everything in 2 nested FOR loop to do my multiplication.

# Import arcpy module
import arcpy
from arcpy import env
from arcpy.sa import *

#Check Spatial Analyst extension
arcpy.CheckOutExtension("Spatial")

# Set the current workspace to iterate on rasters
arcpy.env.workspace = "D:/PROJECTS/03_2014_AB_HYDRO_MODELLING/02_TRANSFORM_DATA/WPWFS/WaterSupplyVulnerability.gdb"

#Variables
PiS = "D:/PROJECTS/03_2014_AB_HYDRO_MODELLING/02_TRANSFORM_DATA/WPWFS/SurfaceWaterUtilities.gdb/SurfaceWaterNetwork/SurfUtsAB_Proj_corr25m" # Water treatment station feature class
Pi = ["POPN_CURRE"] # Get the population value in the field "POPN_CURRE" for each water treatment station
WiS = arcpy.ListRasters("WEI*", "All") # Create a list of rasters
Qn = arcpy.Raster('WaterYield_CubicMeters_SqKm_Year') # Water yield raster

#  For each Station's population Pi in All stations PiS, multiply POPN_CURRE with rasts (119 rasters)
with arcpy.da.SearchCursor(PiS, Pi) as cursor:
    for row in cursor:
        row.next()
        for Wi in WiS:
            print Wi
            WiPi = Wi * row
            WiPi.save("D:/PROJECTS/03_2014_AB_HYDRO_MODELLING/02_TRANSFORM_DATA/WPWFS/WaterSupplyVulnerability.gdb/WiPi" + row)
            WiPiS = arcpy.ListRasters("Wi*", "All") # Create a list to iterate through WiPi calculation rasters
            SWiPi = CellStatistics("WiPiS", "SUM", "DATA")
            SWiPi.save("D:/PROJECTS/03_2014_AB_HYDRO_MODELLING/02_TRANSFORM_DATA/WPWFS/WaterSupplyVulnerability.gdb/SWiPi" + row)
            IMPn = Qn * SWiPi
            IMPn.save("D:/PROJECTS/03_2014_AB_HYDRO_MODELLING/02_TRANSFORM_DATA/WPWFS/WaterSupplyVulnerability.gdb/IMPn" + row)
        arcpy.Delete_management(WiPiS) #Delete intermediate raster from Times processing

It returns the following error :
Runtime error

Traceback (most recent call last):
  File "<string>", line 4, in <module>
TypeError: can't multiply sequence by non-int of type 'tuple'

I understand that my population field is not a integer, and so ?
When I try the 2 for loop with just a print rather than the calculation, it works.
So what is the function or method to call ?

Thanks


texts: Extract corner coordinates of polygon extent in PyQGIS / GDAL

texts:
My goal is to extract corner coordinates (i.e. UpperLeft, UpperRight, LowerRight, LowerLeft) from polygon extent using PyQGIS or GDAL.

In ArcPy there is a handy property of an extent class (i.e. feature['SHAPE@'].extent.upperLeft),
but is there something similar in PyQGIS / GDAL?


texts: Seeking book or site for help with ArcGIS Desktop?

texts:
I am new to the field of GIS and use Trimble Pathfinder and TerraSync.

I just started with ArcGIS 10.3 and was wondering if anyone can recommend a entry level "Q&A" book or site for ArcGIS 10?


texts: Google Earth Engine: memory error when calculating monthly average precipitation for a small region using GLDAS data

texts:
I am trying to calculate the total monthly precipitation rate using the Rainf_f_tavg band of the GLDAS-2.1 image collection on Google Earth Engine (https://developers.google.com/earth-engine/datasets/catalog/NASA_GLDAS_V021_NOAH_G025_T3H#bands). This is the code I have so far. However I keep getting the error user memory limit exceeded. I assume this happens because I am trying to analyse images over such a long timeframe (2001-2020), but I still get the error even when I am analysing a much shorter time frame. I can only analyse 1 single year without getting an error. I used the code in the answer to this question to try and reduce the images that need to be analysed by first computing the daily means Download GLDAS precipitation monthly 1 km resolution data by Earth-engine but it still seems to be too much. Is there any way around this or a way that I can reduce the number of images that need to be analysed further? Or is there a similar dataset I could use which is less computationally intensive?
var startDate = ee.Date('2001-01-01');
var endDate = ee.Date('2020-12-31');
var startYear = startDate.get('year');
var endYear = endDate.get('year');

var years = ee.List.sequence(startYear, endYear);
var months = ee.List.sequence(1,12);

var lst = ee.ImageCollection("NASA/GLDAS/V021/NOAH/G025/T3H")
  .select('Rainf_f_tavg')
  .filterDate(startDate, endDate);

print(lst.size());

var days = endDate.difference(startDate, 'days');
var daysStep = 1;
var timeSeries = ee.ImageCollection(
  ee.List.sequence(0, days.subtract(1), daysStep)
    .map(function (offsetDays) {
      var start = startDate.advance(offsetDays, 'days');
      var end = start.advance(daysStep, 'days');
      var composite = lst
        .filterDate(start, end)
        .mean(); // You need to decide how to combine the images
      return composite
        .set('system:time_start', start.millis())
        .set('empty', composite.bandNames().size().eq(0));
    })
  ).filterMetadata('empty', 'equals', 0);

print(timeSeries.size());

var transformed = timeSeries.map(function (image) {
  var date = image.get('system:time_start');
  var img = image.multiply(86400).set('system:time_start', date);
  return img;
});

//print(transformed);

// make monthly mean mosaics
// loop over the years and months to get summed monthly images
var byMonth = ee.ImageCollection(ee.FeatureCollection(years.map(function(y){
  var yearCollection = transformed.filter(ee.Filter.calendarRange(y, y, 'year'));
  var byYear = ee.ImageCollection.fromImages(
    months.map(function(m) {
      var maxImage = yearCollection.filter(ee.Filter.calendarRange(m, m, 'month'))
        .reduce(ee.Reducer.sum());
      var date = ee.Date.fromYMD(y, m, 1).format("MM_dd_YYYY");
      return maxImage.set('system:time_start', ee.Date.fromYMD(y, m, 1)).rename(date);
      //.set('month', m).set('year', y); // eventually set year and month
  }));
  return byYear;
})).flatten());

// filter the empty one out
var outputMonthly = byMonth.filter(ee.Filter.listContains('system:band_names', 'constant').not())
                    .sort('system:time_start').toBands();

var features = outputMonthly.reduceRegions(donetsk, ee.Reducer.first(), 30);
//print(features);

Export.table.toDrive(features,
"Donetsk_monthly_precip_totals_2000_2005_v21");


texts: Script runs in ArcMap Python window but not in standalone PythonWin or IDLE?

texts:
I am attempting to automate some basic geoprocessing for a web map service. I am using ArcGIS 10.1 for Desktop.

This code was built in ModelBuilder and exported to python.

This runs perfectly in ArcGIS for Desktop.

However it returns error 000732 for both the "Service_Group" dataset and the Excel sheet that are to be joined when ran Outside of Arc.

I have changed the #Local Variables to a UNC path with no luck.

    import arcpy
mxd = arcpy.mapping.MapDocument(r"X:\Mikes_Workspaces\Outage\Outage.mxd")
df = arcpy.mapping.ListDataFrames(mxd,"Layers")[0]
# Script arguments
OutCurrent = arcpy.GetParameterAsText(0)
if OutCurrent == '#' or not OutCurrent:
    OutCurrent = "X:\\Geodatabases\\WebData\\Water_Service.gdb\\OutCurrent" # provide a default value if unspecified

# Local variables:
Service_Group = "Service_Group"
Update_ = "S:\\Pump and Canal Patrol\\Water On\\4-7-15\\Outage_Board.xls\\Update$"
Group_Out = "Service_Group"

# Process: Add Join
arcpy.AddJoin_management(Service_Group, "Group_", Update_, "Group_Out", "KEEP_COMMON")

# Process: Copy Features
arcpy.CopyFeatures_management(Group_Out, OutCurrent, "", "0", "0", "0")
# Process: Symbology
arcpy.ApplySymbologyFromLayer_management("OutCurrent", "X:\Mikes_Workspaces\Online Shapefiles\Outage_today.lyr")
# Process: Remove Join
arcpy.RemoveJoin_management(Service_Group, "")

This is the error received:

    Traceback (most recent call last):
  File "C:\Python27\ArcGISx6410.1\Lib\site-packages\Pythonwin\pywin\framework\scriptutils.py", line 326, in RunScript
    exec codeObject in __main__.__dict__
  File "E:\GIS\Mikes_Workspaces\Automation\OutageUpdate.py", line 15, in <module>
    arcpy.AddJoin_management(Service_Group, "Group_", Update_, "Group_Out", "KEEP_COMMON")
  File "C:\Program Files\ArcGIS\Server\arcpy\arcpy\management.py", line 5325, in AddJoin
    raise e
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Layer Name or Table View: Dataset Service_Group does not exist or is not supported
ERROR 000732: Join Table: Dataset S:\Pump and Canal Patrol\Water On\4-7-15\Outage_Board.xls\Update$ does not exist or is not supported
Failed to execute (AddJoin).


texts:
I recently upgraded to Ubuntu 16.0 and installed QGIS 2.18. As I tried it, this error came:

ERROR 1: NUMPY driver was compiled against GDAL 2.1 but current library version is 1.11

Segmentation fault (core dumped)

I have basic ubuntu/linux capability.


texts:
I am working on spatial point pattern analysis using R.

I have two sets of coordinates, both defined inside a 32*20 region.
I first convert two data sets to ppp objects:

Tumor_CD8 <- ppp(Tumor_CD8[,1],Tumor_CD8[,2], c(0,32), c(0,20))
Tumor_CD3 <- ppp(Tumor_CD3[,1],Tumor_CD3[,2], c(0,32), c(0,20))

And now I am confused about how to combine them to a multitype point pattern object, in order to perform functions like: Gcross, Kcross and Lcross, etc.


texts: Calculating real length and slope of roads

texts:
I have collected points alongside road in form of csv files with z-values. When imported in QGIS I used Set Z value and Points to path to make line. When I use right button in Identify Features tool I can see that all vertices have appropriate z values.
However, when I use $length function in Field Calculator I get exactly same values as if I have draw line but hand without z values.
Also, how can I calculate average slope for same line?


texts: combining st_contains with st_buffer gives very slow query

texts:
I have a layer with roads (4000 features) and a layer with routes (80 features). I want to select every road located on a route, combined with the route-type.

Normaly I would use select road.geom, road.id, route.type from road, route where st_contains(route.geom, road.geom).

Because the routes are not properly digitized I need to buffer the routes to select all roads contained by a route.  So when I do select road.geom, road.id, route.type from road, route where st_contains(st_buffer(route.geom, 3), road.geom) the query keeps running for a long time without result.

How can I make this query running faster?


texts: Toggling feature visibility in Openlayers + GeoExt

texts:
I'm setting up a small web mapping application using PostGIS + mod_tile + FeatureServer + OpenLayers + GeoExt, each of which I'm learning as I go. In addition to the base layer, I've pulled in a few vector layers from FeatureServer, each with many features. For example, a request to my "states" layers returns a GeoJSON feature for each state, and "Congressional districts" returns a feature for each Congressional District. I can render these layers to the map and toggle them in their entirety in a GeoEXT layer tree.

What I'd like to do (seems) straightforward: For each layer, display a grouped list of all features in a layer tree with checkboxes to toggle the visibility/rendering of each feature, individually. I'm struggling to find a straightforward way to do this.

There was a question on the GeoEXT mailing list that included code to set the renderIntent to 'delete' on uncheck:

var layerList = new GeoExt.tree.LayerNode({
    text: 'The United States',
    expanded: true,
    layer: usa, // refers to a vector layer using HTTP protocol to get GeoJSON.
    leaf: false,
    children : [{
        text: "Alabama",
        leaf: true,
        fid:1,
        checked: false
    },{
        text: "Alaska",
        leaf: true,
        fid:2,
        checked: false
    }, {
    //...  load other fifty states
    }]
});

var layerTree = new Ext.tree.TreePanel({
    title: 'Map Layers',
    root: layerList,
    autoScroll:true,
    listeners: {
        'checkchange': function(node, checked){
            var tempfid = node.attributes.fid;
            if(checked){
                node.parentNode.attributes.layer.getFeatureByFid(tempfid).renderIntent = 'default';
                // I am setting visibility with renderIntent because styles are reset when layer is redrawn form panning or zooming.
                node.parentNode.attributes.layer.redraw();
            }else{
                node.parentNode.attributes.layer.getFeatureByFid(tempfid).renderIntent = 'delete';
                node.parentNode.attributes.layer.redraw();
            }
        }
    }
});

I'm thinking there's a better way to do this (perhaps with Strategy.Filter?), but can't seem to wrap my head around it. Could someone point me in the right direction and/or to some examples that include this functionality? Thanks in advance for your help!