texts: PointStacker - Random locations at different zoom level

texts:
I'm applying PointStacker transformation on one of my layers using Geoserver 2.5.2 and OpenLayers 2.

I know, that I have 7 datapoints located on the exact same location.
However, when I change the zoom level the stacked point is moving around randomly.
I would've thought that the stacked point's location is calculated by the mean of the vector points being stacked. In other words, if they are located at the same location, it shouldn't move around changing zoom levels.

What am i doing wrong?

zoom level 7:

zoom level 8:

zoom level 9:

zoom level 10:

zoom level 11:


texts: Mapbox-gl: flyTo feature with certain property

texts:
I have been playing around with mapbox-gl.  I am trying to get the map to flyTo() a specific polygon in one of my layers.

After looking extensively through the API I can only find methods to return features after clicking the map or selecting with a bounding box -- no method to parse the features of a layer and its attributes.

I want the latlon (or bounding box or whatever) of a polygon in the dataset, where "NBHD_NAME" equals the variable I am searching for.  I am using a polygon PBF hosted on Mapbox.  I cannot seem to access the features of the vector layer through map.getLayer() or map.getSource().

Am I missing something obvious, or is this functionality not really present in mapbox-gl?


texts: I want to assign a specific long/lat position on a map to the elevation from SRTM3 data files, but have no idea how to find the specific value. So I want some example of how I can find in N50E14.hgt elevation to 50°24'58.888"N, 14°55'11.377"E.


texts: Filling column of selected attributes with same number using QGIS

texts:
In QGIS, I have selected multiple point features by using "Select by Location" in Vector - Research tool. In the attribute table, and I want fill a particular number (1) for those selected features.
How do I do it?
I tried with

replace('grid_id','NULL','1')


texts:
I'm a student engineer and I'm preparing my final project study about the problems of exchanging Data between CAD and Gis in general - particularly Microstation and ArcGIS  - and the changes that affect the files during their transfer: color's changing, thickness, position...

Moreover, the causes of these changes and how we can proceed to reduce the loss of information when transferring files.

Have you any information about this?


texts: Mapbox GL JS vector tiles show in higher zoom level without enough vector tiles

texts:
If I have a set of vector tile files that only available to max zoom at 16, how can I show it on the map at zoom 17 or even 18?

Because of the cause of the data, or the amount of the vector tiles in a layer at a zoom level , my vector tiles generating speed is slow or generated vector tiles is unavailable because it may effect the display of the other layers. That is, I may unable to generated enough vector tiles when I needed them.

My map setting:

map = new mapboxgl.Map({
                container: "gmap",
                zoom: 8,
                center: [120.685296,24.136738],
                style: {
                    "version": 8,
                    "layers": [
                        {
                            "id": "TempBaseMap",
                            "type": "background",
                            "paint": {
                                "background-color": "#C0C0C0"
                            },
                            "minzoom": 8
                        }
                    ],
                },
                maxZoom: 18,
                localIdeographFontFamily: "'Microsoft Jhenghei'",
                logoPosition: "top-left"
            });

And I want to add this vector tiles to a layer...

map.addLayer({
                    "id": "Road",
                    "source": {
"type":"vector",
"tiles":["http://(my local machine)/Road/{z}/{x}/{y}.pbf"]
},
                    "source-layer": "Road",
                    "type": "line",
                    "maxNativeZoom": 18,
                    "layout": {"visibility":"visible"},
                    "paint": {},
                    "minzoom": 8
                });

maxNativeZoom seems works with mapbox.js but I feel that it didn't work with mapbox GL JS.


texts: Reinitialize ArcGIS10.0 Python installation

texts:
I am developing a Python script for a Toolbox that required a few libraries which I installed manually into C:\Python26\ArcGIS10.0\ and it worked just fine. To make the distribution easier, I dived into the Python installation and libraries files and ended up deleting the wrong files and messing up my installation. I fixed a few things but I don't trust that setup anymore. For example, it runs fine outside of ArcMap (from command line or PyScripter) but raises weird errors when run from ArcMap (with the same parameters) .I know that was a mistake and I should have been more careful.

FYI: The fix that I have done was to follow the advice on Esri's blog article A simple approach for including 3rd party Python libraries with your scripts : "The secret to most Python modules and packages: they’re just files. You can copy any Python library  into the same folder as the .py file that imports them and it will just work". I copied the necessary library files to the same folder as my main script. Result: it's still missing libraries or gives me a "bad argument to internal function error" when run from ArcMap but works fine when run separately.

My question is: is there a way to clean up the Python configuration without uninstalling/reinstalling everything? I would like to restart with a clean config.

This may seem a bit dumb but I've been struggling with that since I messed things up...Also, if you have advice or good practices concerning this kind of things, I'd gladly hear them as I am quite new to developing on Windows...

Thanks!


texts: How to interpolate daily climate data using ArcGIS for Desktop?

texts: I want to interpolate daily climate data available from 3 locations. I am not sure how to interpolate the data for daily values for 30 years from the available 3 locations.


texts: What types of GIS problems does functional programming address?

texts:
After reading @ray vernagus post about using ArcObjects with F#, I have been wondering ... what types of GIS problems are better handled by Functional Programming, as opposed to more traditional forms of programming?

Perhaps more specifically, when would I be better off using F# as opposed to C# with ArcObjects?


texts:
I'm searching for a way to determine the location of wet patches, after rainfall, in a DEM. By this, I mean puddles.

We are investigating the best locations to locate paths for recreational use. Paths that remain wet for a longer time need a lot more maintainance.

I would like to visualise the wet areas, to take into consideration, while planning a new pathstructure.

Does anyone know how I can achieve this?

I would like it to look something like this:
http://www.powersim.com/sitefiles/site4053/img/RainFallMadla3.jpg


texts:
I am trying to create a decent cloudfree mosaic of Landsat 7 Tier 2 TOA reflectance around the Amery Iceshelf in Antarctica. I did find a good script for landsat 8 TOA (Apply a cloud mask to a Landsat8 collection in Google Earth Engine - time series) and those work for me, however when I try to change this to landsat 7 with the correct parameters, it doesn't.

Next to the couldmask issue, the result is a collection, whereas I would like to have 1 (mosaiced) image. Any idea?

This is the code I am using at the moment:

var geometry = /* color: #00ffff */ee.Geometry.Polygon(
        [[[92.57694067490729, -70.98282239796863],
          [57.06912817490718, -70.57777747363254],
          [51.35623754990718, -76.61497412799037],
          [91.96170629990729, -77.28946992862909],
          [92.66483129990729, -71.18228294956876]]]);

// TOA L7 dataset
var mosaic =  ee.ImageCollection("LANDSAT/LE07/C01/T2_TOA")
                  .filterDate('1999-12-01', '2003-03-01')
//                  .filterMetadata('CLOUD_COVER', 'less_than',10)
                  .filterBounds(geometry)
 //                 .sort('CLOUD_COVER',false)
//                  .first()

 var getQABits = function(image, start, end, newName) {
    // Compute the bits we need to extract.
    var pattern = 0;
    for (var i = start; i <= end; i++) {
       pattern += Math.pow(2, i);
    }
    // Return a single band image of the extracted QA bits, giving the band
    // a new name.
    return image.select([0], [newName])
                  .bitwiseAnd(pattern)
                  .rightShift(start);
};

// A function to mask out cloudy pixels.
var cloud_shadows = function(image) {
  // Select the QA band.
  var QA = image.select(['BQA']);
  // Get the internal_cloud_algorithm_flag bit.
  return getQABits(QA, 7,8, 'Cloud_shadows').eq(1);
  // Return an image masking out cloudy areas.
};

// A function to mask out cloudy pixels.
var clouds = function(image) {
  // Select the QA band.
  var QA = image.select(['BQA']);
  // Get the internal_cloud_algorithm_flag bit.
  return getQABits(QA, 4,4, 'Cloud').eq(0);
  // Return an image masking out cloudy areas.
};

var maskClouds = function(image) {
  var cs = cloud_shadows(image);
  var c = clouds(image);
  image = image.updateMask(cs);
  return image.updateMask(c);
};

var mosaic_free = mosaic.map(maskClouds);
var visParamsPanSharp = {bands: ['red', 'green', 'blue'],min: [0,0,0],max: [1, 1, 1]};

Map.addLayer(mosaic, visParams, 'With clouds');
Map.addLayer(mosaic_free, visParams, 'Cloud free');


texts: Add "x" in attribute table if raster is georeferenced (ArcGIS and Python)

texts:
I have a new question and I hope somebody can help me with this.

I have got a shapefile which shows the extent of every raster I have. Some of these rasters are georeferenced but some not. Now I would like to add a x to the attribute table where the raster is georeferenced. I know that there must be a way with python GetRasterProperties but I'm not sure how to write the script. Maybe someone of you did the same and can give me some help.

I have to do this with ArcGIS 10 and it must be easy to update when new rasters are georeferenced.


texts:
I am creating shapefile of buildings based on QGIS background images and need to add information  on the final shapefile about the source and date at which these images were taken.

Can anyone tell me how to know production date of images used as background in QGIS?


texts: Work with projects created in old QGIS version after update

texts: If I install the new version of QGIS, can I work on the projects that I created in the old version ?


texts:
I have been struggling trying to use Tabulate Area in Arc 9.3 for a large raster dataset.

I have one raster with 8 land uses and another vector layer with zones (which I converted to a raster with the same resolution as the raster with 8 land uses, I also used the snap raster function when I did this).

When I use tabulate area separately for the two datasets, I get different areas.
To cross check, I tabulated area using one of the polygons in the vector dataset for tabulate area, and got quite a significant difference in area as well.

Any suggestions on what would be the best way to tabulate area for a large dataset?


texts: Exporting offsets from Trimble Rangefinder?

texts:
I already asked this question at http://www.rpls.com/forums/viewtopic/8/1276880, but didn't get any answer so far.

We have the following Trimble devices:
GPS Pathfinder Pro, Handheld with TerraSync and a TruPulse 360b Laser Rangefinder.

In a (dense) forest, we want to map all trees with 13 hectares (~3000 trees). Some gaps in this forest allow measurements with the dGPS, but mostly there is no reception. By measuring distances, azimuths and slopes between neighboring trees, we could calculate coordinates for all trees, based on the few GPS measurements we take in the gaps.

For this, I can create a feature (without position --> no reception) for each measurement and take an offset measurement from one tree to another with the Rangefinder. I can than export the data file to GPS Pathfinder Office software (V4.2). In the software, I can see the offset for each feature (distance, azimuth, slope) by clicking on it.
But then there's a problem: In order to use the offset for the calculations, I'd have to export it. In this isn't working. It is possible to export a lot of attributes for each feature, but apparently not the offset.

So here are my questions:

Is it somehow possible to export the offset measurements with GPS Pathfinder Office (V4.2)?
If not: Is this possible with the newest version of this software?
If not: Could any other software accomplish this?


texts: Modify circle but don't allow to change center in OpenLayers

texts:
I have a Circle on my map with the following modify code:
var modify = new Modify({
            source: this.getRadiusLayer().getSource(),
            wrapX: true,
            style: new Style({
                stroke: new Stroke({color:'#AA0000', width:2}),
                fill: new Fill({color: '#AA000059'}),
                image: new CircleStyle({radius:5, fill: new Fill({color: '#AA0000'})})
            }),
        })

I want my users to be able to change the radius of the circle, but I don't want them to be able to change the center point of the circle. Is there any way I can disable moving the circle to another center point in version 6.4.3?


texts:
I am very new to mapping and I have loaded a set polygons from a shapefile into ArcGIS Pro. They appear correctly however they have clearly lost (or didn't have in the first place), their longitudinal and latitudinal information as they are currently placed in the middle of the ocean, which I believe is the 0 lat, 0 long position.

I can manually move these polygons but it won't let me save the moves, so when I click 'select', it moves back to the original location. Both of the 'save' and 'discard' buttons are greyed out and cannot be clicked. If I click 'finish' on the move tool then it will snap back to the original location too.

Does anybody have any ideas?


texts: Confusing lines in shapefile in ArcGIS Desktop?

texts:
I have SHP of lines extracted from AutoCAD drawing, but the problem is that the end of the cables are not “Confused” (see picture below).

The result should be something like this [2 separated line but they have the same end (confused end)].

Since I have lot of lines to do I can’t do it manually.

Is there a way to do this automatically in ArcGIS Desktop?


texts: 60GB .gdbtable file in my file geodatabase – what is it, and can I delete or make smaller?

texts: The title says it all really – I have a single file named 'a0000000.gdbtable’ in my default file geodatabase that is 60 GB in size. I have a number of other files with .gdbtable extension that aren't nearly as big. Just wondering what it is and whether I can delete it or make it smaller? It's taking up too much space!


texts: How to define a threshold in Principal Component Image for Urban mapping?

texts:
I am working on an Urban mapping and i would like to know how i could define a threshold value for extracting built-up areas from a PC (Principal Component) Image ?

I am interesting on how to determine a threshold value?

What are the most common methods of identifying thresholds in satellite Image Processing?


texts:
Needs to have from port, to port, and time stamps for each (arrival, departure times)

Is there a dataset with all of this, or will I need to combine some other ones?


texts: Large-scale cloud masking, grouping and mosaic in Descartes Labs platform

texts:
I am looking to create a stack of cloud-free sentinel-2 imagery by day for a large portion of California for all of 2020 in the Descartes Lab platform. I realize that this is a very large amount of imagery and will require significant memory to store and so I'm looking for methods to mosaic imagery by day. I'm new to Descartes and I think using the groupby and mosaic methods could work but not sure how to implement the groupby function on a masked image collection.
First, I establish a small area of interest in Northern California and dates from January to March of 2020 (I am using a small area and small time range to test out the workflow):
import descarteslabs as dl
import numpy as np

# Define a bounding box in N CA, in a GeoJSON
aoi = {
    "type": "Polygon",
    "coordinates": [
        [
            [-121.8430, 38.9676],
            [-121.6991, 38.9676],
            [-121.6991, 38.8590],
            [-121.8430, 38.8590],
        ]
    ],
}

start_datetime = "2019-01-01"
end_datetime = "2019-03-30"

Next, I create a Sentinel-2 scene stack with desired bands:
sentinel_scenes, sentinel_ctx = dl.scenes.search(
    aoi,
    products="sentinel-2:L1C",
    start_datetime = start_datetime,
    end_datetime = end_datetime,
#    cloud_fraction=0.7,
    limit = 500,
)

sentinel_all_stack = sentinel_scenes.stack('red green blue red-edge nir swir1 derived:ndvi',
                                           sentinel_ctx,
                                           processing_level='surface')

Then I create a Descartes Labs Sentinel-2 cloudmask stack for masking:
dlcloud_scenes, dlcloud_ctx = dl.scenes.search(
    aoi,
    products=["sentinel-2:L1C:dlcloud:v1"],
    start_datetime=start_datetime, end_datetime=end_datetime,
    limit=None
)

dlcloud_valid_cloudfree_stack = dlcloud_scenes.stack('valid_cloudfree',
                                                     sentinel_ctx ,
                                                     data_type='Byte')

Next, I mask the Sentinel-2 collection with the Descartes Labs cloudmask:
dlcloud_valid_cloudfree_stack_rep = np.repeat(a=dlcloud_valid_cloudfree_stack,
                                          repeats=7,
                                          axis=1)

# Mask the pixel values where valid_cloud_free is 0
cloudfree_stack = np.ma.masked_where(dlcloud_valid_cloudfree_stack_rep==0, sentinel_all_stack)

From here, I would like to group the imagery by day. I tried:
for (year, month, day), day_scenes in imagery.groupby(
    "properties.date.year", "properties.date.month", "properties.date.day"):
    print("{}: {} scenes".format(month, np.array(day_scenes.stack("blue green red red-edge nir swir1 derived:ndvi", ctx)).shape))

And I get the error: 'MaskedArray' object has no attribute 'groupby'
I would like to mosaic all images grouped by same day afterwards. All ideas are very welcome!


texts: Save as output layer from QGIS processing algorithm to PostGIS database using PyQGIS

texts:
I have two layers(layer1,layer2) in PostGIS database where I have connect in QGIS project.

My database have specific precision. Now I want to use from QGIS processing some algorithms like GISDIFFERENCE(for example I know tew postgis have query for this). But I don't know how to import the algorithm output in PostGIS database automatic to keep my database precision.

For example, if I try to run graphical this tool show me choose to save as to PostGIS database. Because I want to create a QGIS plugin.

Some way is to run algorithms and I save manual the output to shapefile and after to import this shapefile to PostGIS database like this :

outputs_QGISDIFFERENCE_1=processing.runalg('qgis:difference', layer1,layer2,False,None)
load= QgsVectorLayer(outputs_QGISDIFFERENCE_1['OUTPUT'],'myname', 'ogr')
 QgsMapLayerRegistry.instance().addMapLayer(load)

and save to PostGIS database.

But that way create it me topology error because I loose the database precision if I use the graphical way to save as to PostGIS database then that layer is correct.

But how to do this using PyQGIS?


texts: In MapBasic, how do I register a table in XLS format that has a different range from one execution to another?

texts: I need my application to be able to register an XLS table that does not have a set range.  One time it may be A2:AL10000, the next it may be A2:AL2000, etc.  I have not seen this anywhere in the reference manual.  Any ideas?


texts: GRASS GIS - 'cat' directory showing 0 categories for reclassified DEM?

texts:
I am trying to use an online model (Titan2D) which requires a material map (a DEM that has been reclassified so that each group can be assigned a value of constant bed friction angle once it has been uploaded into the model) which must be in a GRASS GIS format, however, I'm having problems... When I upload the material map the model says that there are no 'material names' which are created in the 'cat' folder during the GRASS GIS part. So I think that must be where I'm going wrong?

Within ArcMap, I have reclassified a DEM so that ranges of elevation have been assigned an integer value. I have ended up with the DEM split into 11 groups.
Next, I have used the 'raster to ASCII' tool and imported the new text file into GRASS GIS using the 'r.in.gdal' option.
Next, I go into the 'cat' directory and edit the file so that each of my reclassified groups has a material name, however, to begin with, the file says that there are no categories...

'# 0 categories

0.00 0.00 0.00 0.00'

Could this indicate that the reclassified groups I previously made haven't been picked up?
I've tried editing this file so it reads as:
'# 0 categories
Vector map: Slope

0.00 0.00 0.00 0.00
1:SlopeA
2:SlopeB
3:SlopeC
4:SlopeD
5:SlopeE
6:SlopeF
7:SlopeG
8:SlopeH
9:SlopeI
10:SlopeJ
11:SlopeK'

But this is what causes me the error.

Any ideas?


texts:
Is there an ArcGIS extension or utility that will view 360 degree panoramic images?  I've found utilities that are part of larger software packages for managing mobile scanning data, but all I'm after is the spherical viewing capability.

If there isn't a tool out there, does anyone know of a .NET control that I could incorporate into an ArcMap extension that will handle spherical display?


texts: Use two shapefiles in one Leaflet map in R

texts:
My query along with a minimal working example is already posted here:

https://stackoverflow.com/questions/48924514/using-data-from-2-shapefiles-with-leaflet-in-r?noredirect=1#comment84855712_48924514

I did not receive a reply on stackoverflow so I am posting again.

Basically I have a shapefile for a region and subregions in a country. I want to make a Leaflet map such that each time I hover over a point, I see the boundary of the REGION and  the name of the SUBREGION and REGION.

I am not able to do this in one leaflet map. I can do this in two separate maps.

Can someone help me?


texts:
I'm using the ArcGIS REST API to create an OD cost matrix by using this link https://logistics.arcgis.com/arcgis/rest/services/World/OriginDestinationCostMatrix/GPServer/GenerateOriginDestinationCostMatrix/submitJob?.
However, when using a feature service URL to create the OD cost matrix, using Postman I always get the same error.
When I have origins and destinations equal to an array as specified here https://developers.arcgis.com/rest/network/api-reference/origin-destination-cost-matrix-service.htm, it all works well. However, when I use a feature service URL this happens.
origins = {"url":"http://sampleserver3.arcgisonline.com/ArcGIS/rest/services/Network/USA/MapServer/1/query?where=Name=%27Indian%20Consulate,%20San%20Francisco%27&outFields=Name,RouteName&f=json"}

destinations = {"url":"http://sampleserver3.arcgisonline.com/ArcGIS/rest/services/Network/USA/MapServer/1/query?where=Name=%27Ghirardelli%20Square%27&outFields=Name,RouteName&f=json"}

I always get this error
{
"type": "esriJobMessageTypeInformative",
"description": "Submitted."
},
{
"type": "esriJobMessageTypeInformative",
"description": "Executing..."
},
{
"type": "esriJobMessageTypeError",
"description": "Invalid value for parameter Origins - Details : invalid URL json "
},
{
"type": "esriJobMessageTypeError",
"description": "Failed."
}
]
}


texts:
When I use the WGS 84 CRS the Chukotka peninsula east of the 180 meridian is shown on the left of the map. If I project (on-the-fly) the map using a Lambert Azimuth projection, the Chukotka peninsula is shown as a whole, but the 180 meridian is shown as a line through the peninsula.

Can someone please tell me how I can either:

project my map so that the 180 meridian is not shown, or
tell me how I can hide the meridian in the vector layer, or
merge the two polygons (there does not seem to be a Merge option in QGIS, unlike ArcGIS)

The syntax I am using is:

+proj=laea +lat_0=55 +lon_0=40 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

I am using QGIS 2.14.1 on Windows 10.


texts:
I have a method where I open a geotiff file and print its dimensions

    public void parseAndSaveFile(String pathFileToParse) throws IOException, Exception {

        GridCoverage2DReader reader = null;
        GridCoverage2D image        = null;
        BufferedImage img           = null;

        try {

            File f = new File(pathFileToParse);

            ParameterValue<OverviewPolicy> policy = AbstractGridFormat.OVERVIEW_POLICY.createValue();
            policy.setValue(OverviewPolicy.IGNORE);

            // this will basically read 4 tiles worth of data at once from the disk...
            ParameterValue<String> gridsize = AbstractGridFormat.SUGGESTED_TILE_SIZE.createValue();
            //gridsize.setValue(512 * 4 + "," + 512);

            // Setting read type: use JAI ImageRead (true) or ImageReaders read methods (false)
            ParameterValue<Boolean> useJaiRead = AbstractGridFormat.USE_JAI_IMAGEREAD.createValue();
            useJaiRead.setValue(true);

            reader = new GeoTiffReader(f);
            //reader.read(new GeneralParameterValue[] { policy, gridsize, useJaiRead });
            image = reader.read(new GeneralParameterValue[]{policy, gridsize, useJaiRead});
            //Rectangle2D bounds2D = image.getEnvelope2D().getBounds2D();
            // calculate zoom level for the image
            GridGeometry2D geometry = image.getGridGeometry();

            img = ImageIO.read(f);
            //WritableRaster raster = img.getRaster();
            //int numBands = raster.getNumBands();

            int width   = img.getWidth();
            int height  = img.getHeight();

            logger.info("Width: {}, Height: {}", width, height);

        } catch (Exception e) {
            e.printStackTrace();
        } finally {
            reader.dispose();
            image.dispose(true);
        }
    }

In the calling method, I try to delete the file

Files.deleteIfExists(Paths.get(pathFileToParse));

but I get the error:

Unable to access the file. The file is used by another process.

Because? Did I correctly release the geotiff file?


texts: How to install Python libs for XYtools to read Excel files

texts:
I need to use the XYtools plugin to read an Excel file. After choosing the XYtools incon and indicating an Excel file to read in, a message appears "there is a problem with loading a python module which is needed to read/write Excel files. ... see documentation/help on how to install the xlw and xlwd libraries." I found the site where the Python items reside but I cannot understand at all what to do next. Apparently the needed Python libraries are not installed with XYtools.
I would appreciate being pointed to somewhere that might give some help with exactly what to do to get XYtools to work with Excel files.
Thanks -


texts: Is there an efficient way of sending a large raster image?

texts:
Rather simple question: I have two large raster images that are stored in ArcSDE. I need to send these images to someone else, and post it on their ftp site. What is the most efficient way to send these rasters?

I tried creating a geodatabase and copying to images and pasting into the database, but that took forever. Is there a faster way?


texts: Converting VB codeblock to Python in ArcGIS Field Calculator?

texts:
I have created a tool in Model Builder to plot telecoms data - the user clicks a point on the map, enters a radius and enters an azimuth. The model then creates the buffer and cuts a segment out of the buffer based on the azimuth angle.

In order to calculate the bearing of each line in the segment, I have added some code to the Field Calculator to ensure it takes into account there is 360 degrees in a circle...see below.

Expression:

val

Expression Type:

VB

Code Block:

dim val
If [Azimuth] < 60 Then
val = 360 - (60-[Azimuth])
Else
val = [Azimuth] - 60
End If

I want to publish this tool as a geoprocessing service to ArcGIS Portal, to then use as a widget within WAB.

I cannot publish the tool with the codeblock written in VB and therefore need to convert to Python, can anyone help me with the syntax please?

I'm using ArcMap 10.5.

Python I'm testing now but getting an error of invalid syntax in line 1...

Expression:

calc(!Bearing_1!)

Expression Type:

Python

Code Block:

def calc(!Bearing_1!):
if !Azimuth! < 60
return 360 - (60 - !Azimuth!)
else:
return !Azimuth! - 60


texts: Calculating field from variable using PyQGIS

texts:
using python editor in QGIS, I want to update a field based on a variable calculated previously in my script (x for this example).
I have used the code with success when using an expression based on other fields
e.g:
expression = QgsExpression("field1||'_'||field2")
but I can't get it to work when the expression is based on a variable, not a field.
x = 5
id_new_col= layer.dataProvider().fieldNameIndex("field")
expression = QgsExpression("x")
context = QgsExpressionContext()
context.appendScopes(QgsExpressionContextUtils.globalProjectLayerScopes(layer))
for f in layer.getFeatures():
        with edit(layer):
            context.setFeature(f)
            f['field'] = expression.evaluate(context)
            layer.updateFeature(f)

This results in an empty field.
How can I call the variable into the field calculator expression?


texts: Create a local tile package from online world imagery

texts:
Is it possible to create a local tile package (.tpk file) containing data from an online source - for example something like world_imagery.mxd? ArcView 10.2 allows me to set up such a map document, but when I choose to Share As | Tile Package it says: 00005 Layer Type is not supported. I understand that it would need to download the images.

The point of this is to create an off-line base map. I've successfully created base maps from the EDN data that I downloaded from ESRI, but I don't have satellite imagery. An alternative would be a terrain map that shows shaded hills and mountains, but I don't see that on the edn data either.


texts: Python Zonal Statistics Problems in QGIS

texts:
I am trying to run Zonal Stats for Block Group shape files against NLCD Percent Impervious raster. I have shapefiles for 48 contiguous states all projected to the same as the raster.

I've written this Python script to try and do each state shapefile in a loop. At a Linux command line I run the script, it prints out all the shapefile names to the console and does not an error but the Zonal Stats are never added to the shapefiles. It lists all the shapefile names in a second.

If I try running it from the QGIS Python Console it appears to run the script using execfile() and no error or print statements are given but the Zonal Stats are also not added to the shapefiles. Does anyone know what I am doing wrong? Zonal Stats is running fine from the QGIS desktop.

I see there are a number of other topics on this subject but I haven't found a response in them that solves this issue.

import os, sys, traceback
from PyQt4.QtGui import *
from PyQt4.QtCore import *
from qgis.core import *
from qgis.utils import iface
from qgis.analysis import QgsZonalStatistics

def main():

    try:

        shp_fldr = "/media/psf/Home/Documents/Projects/NLCD/blocks_albers"

        raster = "/media/psf/Home/Documents/Projects/data/2016/NLCD_2016_Impervious_L48_20190405.img"
        f_lst = os.listdir(shp_fldr)

        for fl in f_lst:

            if fl.find(".shp") > 0 and fl.find(".xml") < 0:

                print "start: " + fl

                shp_fl = os.path.join(shp_fldr,fl)

                #specify polygon shapefile vector

                polygonLayer = QgsVectorLayer(shp_fl, fl, "ogr")

                polygonLayer.startEditing()

                # usage - QgsZonalStatistics (QgsVectorLayer *polygonLayer, const QString &rasterFile, const QString &attributePrefix="", int rasterBand=1)

                zoneStat = QgsZonalStatistics (polygonLayer, raster, 'ipv-', 1, QgsZonalStatistics.Mean)

                zoneStat.calculateStatistics(None)

                polygonLayer.commitChanges()

                print "end: " + fl

    except:

        # Get the traceback object

        tb = sys.exc_info()[2]

        tbinfo = traceback.format_tb(tb)[0]

        # Concatenate information together concerning the error into a message string

        pymsg = "PYTHON ERRORS:\nTraceback info:\n" + tbinfo + "\nError Info:\n" + str(sys.exc_info()[1])

        # Return python error messages for use in script tool or Python Window

        print pymsg

if __name__ == '__main__':

    main()


texts: Why the Cartesian coordinates from view.ToMap() are incorrect?

texts:
Does anyone know why the coordinates output from view.toMap() are incorrect sometimes? I simplified the code to be the following. You run it and click the map anywhere and you will find that the output cartesian coordinates from view.ToMap() are sometimes different from the carteisan one shown in base map (the correct one). The Lat/Lon pairs are always correct from view.ToMap().
I actually found that X is correct sometimes if you move the map from left to right or from east to west.
If you click in Asia area and you will see that X is negative (should be positive)! Did anyone notice this before?
<!DOCTYPE html>

    <meta charset="utf-8" />
    <meta name="viewport" content="initial-scale=1, maximum-scale=1, user-scalable=no" />
    <title>Intro to MapView - Create a 2D map</title>

    <link rel="stylesheet" href="https://js.arcgis.com/4.16/esri/themes/light/main.css" />


texts:
I am using Windows 7, and have installed QGIS 2.2, 2.4, 2.8 and 2.10 to try and resolve this issue, but keep on getting the same response:
Minidump written to.....

This occurs after the programme has loaded, and I want to start adding shapefiles. I have tried shapefiles of varying sizes based on the input found here, but with no luck whatsoever. Whether they are large or small, QGIS will crash.

What could be the problem? I have used QGIS happily on other PCs with Win7.


texts:
I've been flummoxed by this query. I have a table of points, and a table of lines, and I am trying to return the point id and line id whenever the point is contained within a buffer of the line (i.e. a single point could be contained multiple times, and I would want multiple records returned).

I have been using the following query:

select 'points'.'id', 'lines'.'id' from points, lines
where Contains(Buffer('lines'.'geometry',50),'points'.'geometry')

I do get a result set returned, and it does give multiple points where it is in more than one buffer. However, the result set is limited to the first record in 'points'. I have tested the buffer selection on its own and it appears to be creating multiple polygons, so I'm unsure where the error is here, assuming this is possible to achieve. I don't think this is a specific Spatialite issue and it's probably my SQL. Any pointers much appreciated.

Thanks

Stephen


texts:
I'm working with code from the Disconnected Islands plugin in QGIS x64 and Python. Normally, the plugin ends by writing the FID of features to a CSV and their corresponding network group ID.

Instead, I would like to be able to find all the network groups with only one feature (= lines with no connections). I'm not sure how to proceed.

Preferably, I would like to be able to do what ArcGIS's Find Disconnected Features in Geometric Network tool does, but I need to use QGIS.

# get the network
layer = iface.activeLayer()
G = nx.Graph()

# construct graph
for feat in layer.getFeatures():
    line = feat.geometry().asPolyline()
    for i in range(len(line)-1):
        G.add_edges_from([((line[i][0], line[i][1]), (line[i+1][0], line[i+1][1]),
                          {'fid': feat.id()})])

# evaluate on connected components
connected_components = list(nx.connected_component_subgraphs(G))

# gather edges and components to which they belong
fid_comp = {}
for i, graph in enumerate(connected_components):
   for edge in graph.edges_iter(data=True):
       fid_comp[edge[2].get('fid', None)] = i


texts: "ERROR 000819: Tool is not assigned" While Running Spatial ETL Tool?

texts:
I have Created a spatial ETL tool to do small operations like adding an ID, updating different fields, etc.

While running it from the Edit Session of Spatial ETL it runs well and good.

http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//00210000003q000000

Now I want to run it from ArcCatalog without opening the tool, but it is throwing an error

and inside geoprocessing results it Shows

  ERROR 000819: THE TOOL IS NOT ASSIGNED

can anyone please help me with this, on how to run the tool without using the edit session.

One more question: Inside the properties of the tool, "Always run in Foreground" option is disabled, because of this I guess I have to disable background processing and run the tool in foreground. The help page " Foreground and background processing" is also not solving the problem.


texts: How to maintain many-to-many relationship between feature class and table with multiple disconnected editors?

texts:
I have an application coded with VB.NET (ArcGIS Engine) and its functionality is:

digitize a parcel and enter information for its owners, so i have a M-N relationship between them.
every parcel has it stones  so an M-N relationship ...

i have two PCs where the application is installed, and I am using them in disconnected mode.

The problem is that I want to be able to export daba and import it between the 2 PCs without losing data or the relationship.

Does anyone have any suggestions for how to do this?


texts: Angle conversion for azimuth angle due east counter clockwise to due south clockwise

texts:
Suppose a given azimuth angle is calculated due east counter clockwise (ccw), and I want to convert it to due north clockwise, so I found the following formula,

# limit the range to 0..360:
r.mapcalc "angle_cw = (450 - angle_ccw) % 360"

now if I want to keep it between -180 to +180 I used the following formula,

# limit the range to -180..180 (again, clockwise from North):
r.mapcalc "angle_cw = (630 - angle_ccw) % 360 - 180"

I want a clear explanation how these two formula actually work, so that I can write due south clockwise azimuth angle conversion formula as well while due east counter clockwise angle is given.


texts: I need to extend my digital elevation map (10m or 30m will do) into Baja from southern CA.


texts: Geoserver/Image Mosaic - A problem when adding multiple Custom Dimensions

texts:
I want to add more custom dimensions. My indexer.properties file is:

Caching=false
TimeAttribute=ingestion
ElevationAttribute=elevation
AdditionalDomainAttributes=indeksi(indeksi)
AdditionalDomainAttributes=parcela(parcela)
Schema=*the_geom:Polygon,location:String,ingestion:java.util.Date,elevation:Integer,parcela:Integer,indeksi:String
PropertyCollectors=TimestampFileNameExtractorSPI[timeregex](ingestion),IntegerFileNameExtractorSPI[elevationregex](elevation),IntegerFileNameExtractorSPI[parcelaregex](parcela),StringFileNameExtractorSPI[indeksiregex](indeksi)

Only one custom dimension (PARCELA) appears in the GUI Geoserver and should be two (PARCELA and INDEKSI)

The database (Postgresql) correctly creates all attributes and values

Any ideas what the problem is? Why is the GUI GeoServer not appear both custom dimensions?


texts: How to calculate an Illumination raster for fire ignition danger analysis?

texts:
I am doing a fire ignition danger analysis using this equation that I found in Millington (2012) http://www.jmecology.com/wp-content/uploads/2014/03/33-42-MILLINGTON.pdf

IR = 4H + 3V + 2I – E

Where:

H = human presence (i.e. within 30m of a road/trails)
V = vegetation
I = Illumination
E = Elevation

It is based in another paper by Salas and Chuvieco (1994) that I couldn't find. My problem is the illumination raster. I found another link where they give some guidelines in how to create it, but I have not found any other reference. http://ibis.geog.ubc.ca/courses/geob370/students/class06/fire/mce.html

Do you know where can I find information related to how to create this raster?


texts: GeoServer GetFeatureInfo on RGBA 70GB raster returns empty feature response

texts:
I have trouble getting feature response from GetFeatureInfo (to get raster values at a point) from GeoServer on very big raster - approx 70GB TIFF. Wierdly enough, the WMS service itself (GetMap) works fine and I am possitive, I am selecting bbox, which has data.
The raster itself is in projection EPSG:31258 while I am requesting bbox in EPSG:4326.
I tried both sending a direct request:
geoserver_url?service=WMS&VERSION=1.1.1&
REQUEST=GetFeatureInfo&
QUERY_LAYERS=layername&LAYERS=layername&
INFO_FORMAT=application/json&
X=10&Y=10&SRS=EPSG:4326&STYLES=&WIDTH=100&HEIGHT=100&BBOX=13.8480,46.6154,13.8570,46.6177

and using the Openlayers preview mode in GeoServer.
Both return a same 200 status response with empty feature set:
{"type":"FeatureCollection","totalFeatures":"unknown","features":[],"crs":null}

Preview mode returns the same but in HTML format.
The layer is set to queryable and other raster layers on same workspace yielded expected response with values of the raster in certain point.
Used GeoServer version: 2.12.0

The logs of the request do not tell me anything suspicious.
2021-01-31 19:37:06,396 INFO [geoserver.wms] -
Request: getFeatureInfo
    BaseUrl = xxxx
    Exceptions = application/vnd.ogc.se_inimage
    FeatureCount = 1
    Get = false
    GetMapRequest =
GetMap Request
 version: 1.1.1
 output format: application/json
 width height: 100,100
 bbox: SRSEnvelope[13.848 : 13.857, 46.6154 : 46.6177]
 layers: xxxx
 styles: raster
    InfoFormat = application/json
    PropertyNames = null
    QueryLayers = [org.geoserver.wms.MapLayerInfo@6b0e0587]
    RawKvp = {FORMAT=application/json, INFO_FORMAT=application/json, STYLES=, WIDTH=100, HEIGHT=100, LAYERS=xxx, EXCEPTIONS=application/vnd.ogc.se_inimage, REQUEST=GetFeatureInfo, SRS=EPSG:4326, BBOX=13.8480,46.6154,13.8570,46.6177, VERSION=1.1.1, X=10, Y=10, SERVICE=WMS, QUERY_LAYERS=xxx, TRANSPARENT=true, FEATURE_COUNT=1}
    Request = GetFeatureInfo
    RequestCharset = null
    Version = 1.1.1
    XPixel = 10
    YPixel = 10
2021-01-31 19:37:06,398 INFO [wfs.json] - about to encode JSON

Has anyone encountered a same issue? Could the raster size be the problem?


texts: arcpy.RegionGroup fails with error 999999. Workspace or data source is read only

texts:
outRegionGrp = RegionGroup(inRaster,"EIGHT", "WITHIN", "ADD_LINK", valToIgnore)

arcgisscripting.ExecuteError: ERROR 999999: Error executing function.
Workspace or data source is read only.
The table was not found. [VAT_RegionG_Con_1]
A column was specified that does not exist.
The table was not found. [VAT_RegionG_Con_1]
The operation was attempted on an empty geometry.
ERROR 010274: Error in writing to table.
ERROR 010067: Error in executing grid expression.
Failed to execute (RegionGroup).

While workspace is at c drive, this error occurs. But while workspace is at other drives, this error does not happen.

import arcpy
from arcpy import env
from arcpy.sa import *

env.outputCoordinateSystem = spatialReferance
env.scratchWorkspace = workspacePath
env.workspace = workspacePath

try:
    checkExtensionStatus = arcpy.CheckExtension("Spatial")

    if checkExtensionStatus != "Available" :
        raise LicenseError(checkExtensionStatus)

    checkOutExtensionStatus = arcpy.CheckOutExtension("Spatial")

    if checkOutExtensionStatus != "CheckedOut" :
        raise CheckOutExtensionException(checkOutExtensionStatus)

    outRegionGrpRaster = RegionGroup(inRaster, "EIGHT", "WITHIN", "ADD_LINK", inputFalseConstantValue)
finally:
    arcpy.CheckInExtension("Spatial")


texts: QGIS Line between two points based on topography

texts: Is it possible in QGIS to draw a line between two points and take the topography into account? For Example when there's a mountain or a lake between the two points, the line goes around that mountain or lake


texts:
ERROR 002760: This service definition file was created for a hosted layer and must be uploaded directly to a portal using the portal website or ArcMap.
Failed to execute (UploadServiceDefinition).

I am getting the above error when I am trying to publish a mxd to a portal as a feature service

import arcpy, os

# define local variables
wrkspc = 'C:/Ronit/Projects/DU/Hosted/'
path2MXD = r"C:\Ronit\Projects\DU\Hosted\temp.mxd"
mapDoc = arcpy.mapping.MapDocument(path2MXD)

con = 'GIS Servers/connection.ags

service = 'DemoHost'
sddraft = wrkspc + service + '.sddraft'
sd = wrkspc + service + '.sd'
summary = 'Population Density by County'
tags = 'county, counties, population, density, census'

arcpy.SignInToPortal_server("username", "password", "")

# create service definition draft
analysis = arcpy.mapping.CreateMapSDDraft(mapDoc, sddraft, service, 'MY_HOSTED_SERVICES', "", "", None, summary, tags)

# stage and upload the service if the sddraft analysis did not contain errors
if analysis['errors'] == {}:
    # Execute StageService
    arcpy.StageService_server(sddraft, sd)
    # Execute UploadServiceDefinition
    arcpy.UploadServiceDefinition_server(sd, con)
else:
    # if the sddraft analysis contained errors, display them
    print analysis['errors']

os.remove(sd)

"""


texts: ArcView 3.x Avenue Bitmaps (Tabs?) vs. ArcView 10 Python Cursors

texts:
Note: While this question has an answer, any further tips for optimizing a cursor process would be greatly appreciated. I will be monitoring for any updates.

Currently, my boss (who works in Avenue) and I (working in Python) are both attempting to solve the same problem. Rather, we have both solved it, but the speed at which our solutions operate are... disjointed, to say the least. What his script processes in 2 hours can take mine up to 6. The only real difference in syntax and implementation in logic comes from 3.x's Bitmaps and 10.x's Cursors. We both:

1) Store values from Table 1.
2) Use those values to query a row in Table 2.
3) Store values from Table 2 for insertion into Table 3 as a new row.

In both scripts, these processes are completed in two nested loops. Before I start digging into the wonderful world of code optimization, is this an expected occurrence when comparing Avenue script performance to Python? This is not the first time his scripts have greatly outperformed mine in terms of operation time, so I would like to know if there is something I should be aware of before I crucify myself for horrid scripting.

Here's my script sans extraneous bits:

import arcpy
import time
import sys
import os

def recordfindcopy(inFile,query,outFile):
    findRecord = arcpy.SearchCursor(inFile,query)
    for record in findRecord:
        copyRecord = arcpy.InsertCursor(outData) # <--- D'oh! (See answer)
        field = record.FIELD
        copy = copyRecord.newRow()
        copy.FIELD = field
        copyRecord.insertRow(copy)

StreetsFileList = [r"Path",
                r"Path"]

for sfile in StreetsFileList:
    inStreets = sfile
    inTable = r"Path"
    outData = r"Path"
    fsaEntry = arcpy.SearchCursor(inTable)
    for row in fsaEntry:
        id = row.ID
        sQuery = "ID = %s " % (str(id))
        recordfindcopy(inStreets,sQuery,outData)

EDIT: Given some of the comments so far, I wonder if there might be a better way to do this via joins, although I am dubious given the brobdingnagian (word of the day!) size of the tables. The heart of the processing is to append information from one table to any matching records in a second table and create a third table only containing the important fields. I wanted to try this using SDE, but that appears to not be an available option. Thoughts? I apologize if my questions are always so involved, but I'm trying to get to the bottom of a long-standing annoyance.

Answered: Jakub's simple suggestion alone decreased the processing time from 30 seconds per 500 records to 3 seconds per 500 records. Re-initiating the insert cursor on every insert slowed things down considerably (obviously). While this may not be the most optimization one can do for this process when put up against ArcView 3.x's speed, it is enough for my purposes at this time. Further suggestions are very welcome!


texts: uDig Spatial ToolBox NetworkAttributeBuilder - strahler order

texts:
I just upgraded my OsGeo4W, which upgraded uDig, I was hoping that would solve my issue with NetworkAttributeBuilder -stating it created stream order, but hung up of Pfafstratter..

Console output states:

Finished.
Calculate Strahler...
Finished.
Calculate Pfafstetter...

this would produce zero outputs - however, now that i have upgraded, I now have uDig 1.4.0b, and NetworkAttributeBuilder (& OmsNetworkAttributeBuilder) does not have Strahler as an output - i cannot find Strahler order in any of the Spatial Toolbox tools?  I have the same jgrasstools-0.7.7 Snapshot - i was running before I updated. OsGeo4W...

Anyone have similar issues?

In addition: I cannot get the Spatial Toolbox to run anything as I get a JDE - JDK error - the 'latest' uDig downloaded with OSGeo4W -app/uDig/ folders do not have the JDE folder as older versions did?


texts: Importing CSV data with Coordinates in one column

texts:
I'm trying to import the Lat/Long coordinates from the "Location" field of a CSV from here :https://data.colorado.gov/Water/DWR-Water-Right-Net-Amounts/acsg-f33s

Unfortunately the coordinates are in a single column and I can't seem to get QGIS to recognize them.

There is a set of columns with UTM, but I think most of my other data is in Lat/Long and I am using Natural Earth basemaps and I'm not sure how to convert the entire world data set to Colorado UTM coordinates.

My goal is to be able to load the CSV data into QGIS and place them accurately over a basemap.


texts: ArcGIS ModelBuilder not saving intermediate stage data?

texts:
In ArcGIS 10.3.1 Models I have set my workspace paramters to save intermediate data in a scratch database. When I run the model out of the Modelbuilder Editing window the intermediate stage data is saved. If I however simply double click the model in ArcCatalog only the finished results are saved and no intermediatary data is saved in my scratch database.

Can anyone explain to me why this is, and how I can change it?


texts: How to automatically extract borders from an old map?

texts:
I am new to QGIS and I am wondering whether there is a plugin which automatically recognises the borders on a map and converts the areas into polygons?

Edit: I am using QGIS 2.14.1. It is a picture of an old map with specific borders which don't exist anymore:

I was wondering whether there is a plugin that can automatically detect all of the borders so that each area is then transformed into a polygon.

So far I have only found plugins like autotrace and QGIS' advanced trace digitizing option where I would have to do it manually for each area.


texts:
Is it possible to give a value of 0 when the field is null?
You will tell me to give a default value I did it.
I'm working on a data collection project using qfield. I gave a default value to my field so that it would be always filled, yet operators delete this value and forget to give another.
So I want to populate the value of 0 when it's null.
If you have any idea tell me please.
Regards.


texts: Opening .vrt file with rasterio in Python?

texts:
I have 3 raster files (.tiff files) which are geographically 3 consecutive tiles  of a particular area. I am trying to make a .vrt file of 2 out of these 3 files using the gdalbuildvrt python command. The issue I am facing is suppose if I am making a .vrt file (named test.vrt) of the first 2 files (say file 1 & 2) and opening it with rasterio python, it's opening correctly. Subsequently, I want to create a .vrt file of the new file combinations (say file 1 & 3) and overwrite the old .vrt file (test.vrt) with this file combination and try to open it with rasterio python. The problem I am facing now is it's still displaying the the .vrt file with the old file combinations (having files 1 & 2 and not files 1 & 3). The code which I used till now is:

import rasterio
import rasterio.plot

from osgeo import gdal

vrt_options = gdal.BuildVRTOptions(resampleAlg='cubic', addAlpha=True)
test_vrt = gdal.BuildVRT(r'D:\Test TIFF\test.vrt', [r'D:\Test TIFF\Tile_58_LULC250k_2015_16_revised.tif', r'D:\Test TIFF\Tile_60_LULC250k_2015_16_revised.tif'], options=vrt_options)

test_vrt = None

raster_vrt = rasterio.open(r'D:\Test TIFF\test.vrt')

rasterio.plot.show(raster_vrt)

Also, if I am opening this new .vrt file with QGIS its displaying the correct file combinations (file 1 & 3). So what I understood is the .vrt file has changed after running the gdalbuildvrt command but there is some problem while opening it with rasterio python because it's displaying old .vrt file (file 1 & 2).

Can someone help me figure out what is the problem?


texts: Easy way to render TileMill tiles in retina mode with Leaflet?

texts:
I have a few quick little map websites that run straight off TileMill, without pre-exporting. They usually have code like this:

L.tileLayer('http://example.com:20008/tile/MyMap/{z}/{x}/{y}.png',
  {}).addTo(map);

How do you make this take advantage of a Retina type high resolution display, without just making everything (especially text!) tiny?


texts:
I'm trying to consume a GeoRSS feed from a C# project. Can I have queries against a #georss feed?

For example 'nearest item to Lon/Lat' or 'items contained in bounding box X'?


texts: Dividing tunnel into segments of 100 meters using QGIS?

texts: I have a csv file for pavement survey (potholes). It maps the pavement based on number of potholes in each segment divided into equal 100 meters. I am attaching a picture of it to view. How do I do this on QGIS? I tried different plugins but wasn't anywhere near to what I was looking for. RED indicates an area having 11-15 potholes, YELLOW indicates 5-6 potholes whereas GREEN points to 1-2 potholes.


texts: Is it possible to unit data of several polygons to a new polygon that includes the others?

texts: I have a lot of little polygons (block of houses, sezioni censuarie in Italy) that contain resident in a city. I have united the shape of the polygonsin Arcgis10, with Geoprocessing>UNION, to a grid of 1 kmq, so in every cells there are many blocks or sezioni censuarie, completes but not complits too. Every block have a population: how I can attribute the right pop, of the divided block to the grid? So I can redistribute the pop in the grid and  show this with ArScene


texts:
I  am looking for a way to visualize data from Postgres in Google Maps, but I have only figured out how to do this for Markers.
I am using PHP to retrieve the data in JSON, which I then request from Google Maps, but I can only find examples for doing this with point data.

  var geocoderFunction = function () {
        geocoder.geocode({ 'address': road[index] },
           function (results, status) {
              if (status == google.maps.GeocoderStatus.OK) {
                 new google.maps.Marker({
                    map: map,
                    position: results[0].geometry.location
                 });
           }
   }


texts: Selecting neighbours in an iteration using Select by Location

texts:
In the 'Select by location' tool, there is an option for Selected features only. When it's ticked, the resulting code for selecting touching polygons of a selected feature looks something like this:

processing.run("native:selectbylocation", {
            'INPUT': boundary_layer,
            'PREDICATE':[4],
            'INTERSECT':QgsProcessingFeatureSourceDefinition('gadm36_AFG_2_3040aec4_9c3d_472d_aad3_311351d800e8', True),
            'METHOD':0})

I need to use this code in a loop, where at the bottom of each loop, more features are selected. I don't know how to address this in the 'INTERSECT', where currently a static address is given to QgsProcessingFeatureSourceDefinition. One can use the following, which is not suitable for iteration.

'INTERSECT':QgsProcessingFeatureSourceDefinition(boundary_layer.id(), True),

In the documentation, it is mentioned:

  Parameter type: QgsProcessingParameterFeatureSource

    Accepted data types:
            - str: layer ID
            - str: layer name
            - str: layer source
            - QgsProcessingFeatureSourceDefinition
            - QgsProperty
            - QgsVectorLayer

But I don't know which one is suitable for this purpose.

Context: In a loop, I am starting with one polygon, find its neighbours using 'Select by location', touching, selected features only. Then I choose one of the neighbours based on some criteria and append its ID to a list, which already includes the starting polygon's ID. At the end of the loop, I choose both features using             layer.selectByIds(list). Going back to the beginning of the loop, now the neighbours of these two features have to be selected.


texts: Is there any way to recover data when a .qgz file will not open?

texts:
I saved my QGIS project as a .qgz file and then my computer crashed (however I can't confirm that the saving had completed successfully). After restarting the computer, however, the .qgz archive no longer works and crashes QGIS when I try and load it.

The QGIS desktop window begins with the ~"Loading project:" project tracker with the green bar in the bottom panel but stops immediately and freezes the program to the "crash report" window.

If I try and open the archive with 7zip I get an error message "Can not open file '...' as archive".

Is there any way to access or repair information from a .qgz file/archive?

alternatively ... does QGIS save a backup of the archive or just the project file on the computer somewhere similar to the .qgs~ file in the old system?

I get the QGIS crash report shown below:

Stack Trace

zip_source_function_create :
zip_source_function_create :
QgsZipUtils::unzip :
QgsProjectArchive::unzip :
QgsProject::unzip :
QgsProject::read :
QgisApp::addProject :
QgisApp::openProject :
main :
BaseThreadInitThunk :
RtlUserThreadStart :

QGIS Info
QGIS Version: 3.2.2-Bonn
QGIS code revision: 26842169e9
Compiled against Qt: 5.9.2
Running against Qt: 5.9.2
Compiled against GDAL: 2.2.4
Running against GDAL: 2.2.4

System Info
CPU Type: x86_64
Kernel Type: winnt
Kernel Version: 6.1.7601


texts: Creating 3D geographical maps for physical fields in qgis

texts:
We need to reproduce 3D maps of a certain physical property on a geographical area changing with depth (in the case of the example these are surfaces of equal subsurface temperature). The example graph from a scientific publication looks as follows:

How it is possible to reproduce this with qgis or, possibly, with an equivalent software?
The input data can be in the form of measured subsurface temperatures (or other property, for example, subsea temperature, salinity or whatever) at certain points sampled at some depths, and the surfaces are to be calculated by interpolating them, like it is done for bathymetric curves.


texts: How to open .hgts file from NASADEM data

texts:
I have the NASA SRTM1 data as .hgt files from their ftp website, but I would have these elevation data in respect to WGS84. I found by NASA's ftp site a folder named srtmOnly that contains what I want, but they have the extension .hgts not .hgt

I found just these infos in readme file

  Directory: hgt_srtmOnly

  product description: SRTM-only floating point DEM (*.hgts)

  product data type: 4-byte real

  hgt_srtmOnly contains heights relative to the WGS84 ellipsoid

  hgts  meters (with value set to -32768 for a void in hgt_srtmOnly hgts
  files)

I tried to open it with Qgis as raster but it didnt be recognized, also I have a code c++ that converts .hgt file to .csv also doesnt work with these .hgts file.

Does anyone have an idea how to get these data from these .hgts file?

Add
This is my cread .hdr file as adviced by @camilo

BYTEORDER M
LAYOUT BIL
NROWS 3601
NCOLS 3601
NBANDS 1
NBITS 32
BANDROWBYTES 14404
TOTALROWBYTES 14404
BANDGAPBYTES 0
NODATA -32768
ULXMAP 34
ULYMAP 36
XDIM 0.00027777777777778
YDIM 0.00027777777777778


texts: Expression that selects all objects that share the same value in field 1, among them selects the object with the highest value in field 2

texts:
I have a table with thousands of objects but created an example Image for this question.

I need an expression that does the following.

Select all objects that share the same values in field 1.
Among them select the one with the highest value in field 2.
From this object take the value of field 3.
Add this value to field 4 from all selected objects from step 1.


texts: osm2pgsql vs imposm rendering performance

texts:
Is there any difference in rendering performance between using a  osm2pgsql vs imposm imported PostGIS DB in Mapnik? Or asking it differently, is there any downside of using imposm if I only want to support full imports?

I'd be using OSM Bright style, which supports both, and I've seen that there is a big difference in import speed between the two, but I haven't found any information about how is the rendering part.


texts: Saving spatial reference of MapInfo TAB-file to file (like .prj)

texts: I have some MapInfo TAB-file and I'd like to save spatial reference to .prj-file. How it can be done?


texts: Google Earth Engine Error: Number of pixels requested from Image.load exceeds the maximum allowed

texts:
I have created a cumulative cost layer for the U.S. of distance to rivers based on the National Hydrologic Dataset. My source data is a 30 meter resolution image of all rivers with values equal to the Stream Order number and all non-streams equal to zero. My cost data is this same layer remapped so that all of the U.S. is equal to 1. The cumulative cost image runs fine and I can display it by:

Map.addLayer(cumulativeCost, {min: 0, max: 5e4}, 'accumulated cost');

When I try to export it, I create a geometry of a small section of the U.S., about the size of the state of Washington, and run:

Export.image.toDrive({
  image: cumulativeCost,
  description: 'CostAllRivers1',
  scale: 30,
  region: geometry,
  maxPixels: 1e13});

In the Tasks tab, I receive the error "Error: Number of pixels requested from Image.load exceeds the maximum allowed (2^31)" and in the Console tab I receive "accumulated cost: Tile error: User memory limit exceeded."

I cannot find these errors listed in the documentation, and it happens even when I make the geometry super small. I get the same error if I try to export the cumulativeCost image to an asset. Can anyone help?

I saw mention of using Export as helping with computation time-out errors, but I don't know how to implement this with my cumulative cost function.

From https://developers.google.com/earth-engine/debugging

Bad — don't do this!

var ridiculousComputation = ee.Image(1).reduceRegion({
  reducer: 'count',
  geometry: ee.Geometry.Rectangle([-180, -90, 180, 90], null, false),
  scale: 100,
  maxPixels: 1e11
});

// Error: Computation timed out.
print(ridiculousComputation);

What this error means is that Earth Engine waited about five minutes before stopping the computation. Exporting allows Earth Engine to perform the computation in an environment with longer allowable running times (but not more memory). As the return value from reduceRegion() is a dictionary, you can use the dictionary to set the properties of a feature with null geometry:

Good — use Export!

Export.table.toDrive({
  collection: ee.FeatureCollection([
    ee.Feature(null, ridiculousComputation)
  ]),
  description: 'ridiculousComputation',
  fileFormat: 'CSV'
});

Lara


texts: Forcing virtual field recalc in attribute table in QGIS

texts:
I have a layer with a virtual field, "LongLabel", an expression concatenating several of its data fields. When I edit a feature's data fields in the layer's attribute table/form, how can I force the virtual field to update for that feature? I care since "LongLabel" is the feature selector field on the left hand side of the form.

I'm unsurprised "LongLabel" doesn't update immediately, while I remain in that feature in the table. But I was hoping it would update after I switch to a different feature and come back, or when I save my edits and/or reload the table (icons at top), or - at worst - if I toggle editing on and off to force a data provider commit. No change.

The edits to the underlying data fields are correctly captured, it's just the virtual field that doesn't recalc. After the data edit, if I open a 2nd attribute table (F6 from canvas) it correctly displays the updated "LongLabel" in the 2nd table but the unupdated one in the original table.

If there's no way to do it from the QGIS UI, I'd be fine with a PyQGIS solution. An iface.activeLayer().dataProvider().forceReload() from the Python console doesn't help though. An .activeLayer().commitChanges() followed by .startEditing() toggles editing on and off, reflected in the attr table's editing icon, but still not updating the field. I can also find the attribute table widget by searching for
'QgsAttributeTableDialog/'+iface.activeLayer().id() in [x.objectName() for x in QApplication.instance().allWidgets()] but can't figure out how to force it to do anything other than .close() once I've located it.

Editing to add:
I've gone and replaced "LongLabel" with its expression definition in the layer's display name property. In this case, the LHS feature selector (that used to be "LongLabel") correctly updates when I move away from the edited feature. However, "LongLabel" as a field in the attribute form still does not update. So it's related specifically to virtual fields. Using QGis 3.10.0.


texts: How to produce shapefiles on census block groups by commuting zone?

texts:
I need census block group (CBG) shapefiles separately for each commuting zone (CZ) in the US (by the 2000 delineation, there were 709 of them).

I only found CBG shapefiles by state though, here.

By county-level FIPS, this can be linked to commuting zone delineations, e.g. from USDA ERS here.

So the logical steps would be the following:

Append CBG shapefiles into one database, with county as an attribute.
Join with CZ delineation by county FIPS.
Export separate shapefiles for CBGs by CZs in the aggregate data.

I could not get this done in QGIS. Is it simple with PostGIS? How exactly?


texts:
I have a line shp file from my friend, and want to use the QChainage to draw points every 10 km.
But the Start point and end point isn't right for me.

How can I change the Start point and End point? Search for a day and can't find the solution.


texts: Unknown coordinate system produced in 1980 by Japanese company in North Philippines

texts: I am having the same difficulty with one of the report produced in 1980 by Japanese company. The grids format also N 38,578, E 41,400. The project is 390 kilometers north of Manila. I would like to know what is the coordinate system? What grids? what datum? and how convert them to WGS1984 UTM zone 51


texts: Plot beachball (focal mechanism) for earthquake in QGIS?

texts:
I'm drawing earthquakes on the map using QGIS and I need to plot focal mechanism by beach balls.   I don't want to use svg because it need to transform new locations on the map.

Is there a plugin or tool that I can use for earthquake problems like this?


texts: Arc Hydro toolset to delineate surface drainage area

texts:
I need to determine the total amount of surface area that will drain to a given point in an urban environment. It seems that the Arc Hydro toolset can delineate the watershed that will drain to a given stream, but is there a tool that can be used to determine where precipitation runoff will go or the total surface area that will contribute runoff to a certain point?

To clarify, I don't need to calculate the amount of water that will drain to a point, just the surface area that will contribute to that amount.


texts:
Is there a way to read a TransCad Matrix file (mtx) without access to Transcad.

I am trying to convert a matrix file to a csv file.


texts: Iterate through point feature class and apply tool to specific groups of features

texts:
I am new to Python/ArcPy.

Currently I am trying to iterate through a point feature class and apply the tool "arcpy.AggregatePoints_cartography" to the group of point features for every new time step and as well if the ID (field "ID" not "OBJECTID" changes. The attribute table looks like this:

The field "Monat" includes values up to 12. And the ID-field varies between the values 1 and 2.

I tried to use a search cursor in combination with a counter "LaufID", but the Output Feature class only contains two polygons, but there should be 24 polygons. Please help?

import arcpy
arcpy.env.overwriteOutput = True
arcpy.env.workspace = r"C:/Users/.../Geodatabase.mdb"
laufID = -1

with arcpy.da.SearchCursor("Schatten", ["ID","Tag","Monat"]) as cursorP:
    for pnt in cursorP:
        ID = pnt[0]
        day = pnt[1]
        mon = pnt[2]

        if laufID == -1:
            laufID = ID

        if laufID != ID:
            arcpy.AggregatePoints_cartography("Schatten", "Schatten_Aggr", "100 meters")


texts: Creating separate list of x and y coordinates from array after using getPart?

texts:
I have an array of points by using getPart on each of the four polygons in the shapefile. I am trying to figure out a way to create a list of x and y coordinates in my for loop and name them xCoords and yCoords so I can plug them into my functions to calculate the perimeter. For some reason, I can't seem to make a list that separates the x and y points... any ideas?

Note: When I print point.X or point.Y, it gives me one coordinate. I know that I have to write some code in my for loop to get separate lists of x and y coordinates for each polygon, I'm just not sure how to.

import arcpy, os
from arcpy import env
from arcpy.sa import *
from math import *

env.workspace = r"H:\GEOG 4303\cho_lab4\lab4\data"
env.overwriteOutput = 1

theme = 'watersheds_3D.shp'
DEM = 'demlab4'

cursor = arcpy.da.SearchCursor(theme, ["SHAPE@", "Id"])
for row in cursor:
    array = arcpy.Array()
    print("Polygon Id: {0}".format(row[1]))
    for Id in row[0].getPart():
        print("Polygon coordinates: {0}".format(Id.count))
        for point in Id:
            print("X: {0}, Y: {1}".format(point.X, point.Y))
            array.add(arcpy.Point(point.X, point.Y))
        print array


texts: Expecting ArcPy to overwrite files in geodatabase file?

texts:
I'm writing a python script, I've got an mxd and its workspace is a geodatabase file with all the files I need into.
The problem I'm encountering in testing and writing code is that, despite of the command:

arcpy.env.overwriteOutput=True

all the files contained in geodatabase fail to overwrite.
There's a way to work around this?


texts: Voronoi / Thiessen polygons doesn't create a part for 1 of 4 input points

texts:
I am using QGIS to create Voronoi / Thiessen polygons for 4 rain gauges.

I am having issues with both the QGIS 'Voronoi polygons' tool and the Grass equivalent 'v.voronoi' where one of my 4 points is skipped and only 3 polygons are generated in the output.

As a test I have added an additional point and have found that both tools output a polygon with the correct number of parts (5). Is there a reason it will not work correctly on 4 or less points?

Edit - image showing the configuration of points and the output Voronoi polygons skipping the central point


texts: Creating a segment for fixing water pipes and finding disconnected areas

texts:
I'm wondering if there is a tool for this problem. I'm using QGIS 2.14.2. Basically I want to select a pipe and find closest valves there is and to turn them off. After that I need to find all pipes that are disconnected from the main grid.

I need to solve which users/pipes are left without water when closing a segment from water network to fix a pipe.

This is a very simplified picture of water network, where blue dot marks the water source, red dots are valves and yellow dots users.

So when I select the broken pipe inside the area that needs to be closed, surrounded by valves, it would take them off the grid and show all the pipes that are not connected to the water source. In reality there is usually more pipes that need to be shut down for fixing the broken pipe.

For finding the disconnected pipes I'm thinking disconnected islands tool.


texts: Reference HTML hyperlink with a field using maptips

texts:
I'm attempting to make QGIS maptips display a hyperlink as the text from another field in the attribute table. Below are 3 attempts of HTML and the different results.

When "title documents" is not NULL

When "title documents" is NULL

1

[%case when "title documents" is not NULL then '' end %]
[%case when "title documents" is not NULL then ' "title no" ' else 'Not available' end %]

2

 [%"title no"%]
[%CASE WHEN  "title documents" is null then 'Not Available' END%]

3

[%case when "title documents" is not '' then '' else '' end %]
[%case when "title documents" is not NULL then 'Click Me' else 'Not available' end %]

I'm trying to get result 3 but rather than 'Click Me', make it display "title no" which I would have thought code line 1 would achieve but obviously not.


texts: QField map disappears when centering to GPS location?

texts:
When I load a new project into QField the base map is displayed, I can digitise on and the vetor layer and it updates fine with attribute data. However the GPS position doesn't show until I click the 'canvas follows location' button, but at this point the base map diappears?
When I zoom to layer the base map reappears but the gps position disappears again!
Any one have any suggestions?


texts: Convert from WGS84 to GDA94/MGA94 on iOS

texts:
I have a mobile app (iPhone/iPad) for recording data and taking notes, including photographs, and a government client needs GPS coordinates attached to each photo.

They require coordinates in the GDA94 and/or MGA94 specification, but I can't find any concrete instructions how to convert from a standard lat/lon into those two — and sadly many of the references I did find are out of date or refer to 404 pages.

Is there a simple way to do the conversion? It doesn't need to be accurate since I'm just using a smartphone GPS chipset. But I'd like something with better than 1m accuracy.

As far as I can tell, WGS84 and GDA94 were identical in 1994 but due to continental drift are around 1.4 metres apart now (iOS GPS coordinates are in WGS84 aren't they? I couldn't find any definitive documentation for that either).


texts:
I am creating a spreadsheet of data to enter into QGIS.  My grid references are currently formatted like this: 339614, 133960 - they are the coordinates in the brackets on this site https://www.somersetheritage.org.uk/record/10571
But QGIS isn't registering them, I am struggling to find/work out what coordinate style QGIS accepts and how I can convert coordinates I have into the style required.
Please note that my spreadsheet, which a copy is saved as a .csv, is in excess of 1000 entries and is still growing.  So being able to convert them in bulk would be handy


texts: Clickable markers on PDF generated by qGis

texts: I have created a map in the qGis print-composer with point-data. I would like to display the label of each specific point if a user clicks on it. What's the best way to achieve this?


texts: How do I add a Natural Earth label layer to a TileMill project?

texts:
I'm using TileMill from MapBox to generate tiles and export Mapnik xml files. So far so good. I've downloaded the entire Natural Earth dataset, and have so far successfully imported layers for lakes, coastlines, land and rivers. Now I would like to add labels, using the "Physical Labels" dataset.

When I attempt to add a layer using the shape file 'ne_10m_geography_regions_points.shp', I get a few dozen little red dots (all of which seem to be on a coastline), but no labels. I've tried combinations of the 'ne_10m_geography_regions_polys.shp' and points, but don't get anything nearly resembling the examples on the Natural Earth download site.

I'm a bit new at this, so I'm not sure if the problem is TileMill, the Natural Earth data, or me. I suspect it's me but I don't understand why this layer would be different than the others I've successfully used.


texts: Adding Cost to Shortest Path Finder in FME

texts:
I'm in the process of creating an accessibility map that identifies where a pedestrian sidewalk network overlaps with different features along the network  (e.g. whether the network overlaps with stairs, ramps, crosswalks etc.). I've flagged segments of the pedestrian network with a "1" if they overlap with these features. The first screenshot below shows an example of the data table that I'm working with.

I'd like to find the shortest route between buildings that account for 1) length and 2) that add an additional cost for segments of the route that overlap with inaccessible barriers (e.g. stairs).

I've created a workflow in FME below where I've been able to find the shortest route, but I can't figure out how to add a cost that would indicate to not optimize the route if stair_flag==1. Does anyone know a workaround?


texts: GeoReference tiff image into the Google map?

texts:
I want to Geo-reference a tiff image contain a Map for a specific region, I need to add this region(image) and Geo-reference it on the Google map, SO I can see on the map for the region I select the image which I have inserted.

the image also need information such as the Latitude and Longitude of corners and size, and the pixel per meter measurement.

How to do the correct Geo-reference for that ?


texts: Selecting feature (point) based on specific attribute value(s) of another layer using ArcGIS Pro

texts:
How do i select a feature (point) based on a specific attribute value(s) of another layer?
Scenario:
One layer (map) polygons representing disease risk across a continent, with several columns representing different attributes (risk scale 1-4 etc.) and another layer with point features representing individual movement over time.
I would like to select point features which are overlayed on polygons with a risk scale above 2.5.


texts: GIS developer blogs ... recommendations?

texts:
I'm looking for recommendations for GIS blogs ... specifically GIS developer blogs. I'm less interested in the geography than in how apps are put together.

I'd like to keep up with ESRI ArcGIS as well as open source tools.

Which blogs are you fine people reading?


texts: How to compare aerial and satellite images for land use analysis?

texts:
I am interested to compare the land use-cover dynamics for four periods: 1950s, 1980s, 1990s and current for a certain area having a max size of 5000 ha.

For the first three periods detailed aerial photo of scale 1:20 000 are available while for the current no photo is available and the only possible option I can use is the coarser scale images such as LANDSAT images which can be downloaded for free. I have no fund to purchase detailed images.

My question is how can the information obtained from the more detailed photos be comparable to the one to be obtained using the LANDSAT images for such small area?

I would appreciate very much if anyone can suggest possible outcomes and consequences of this analysis or any suggestion for a way out?


texts: Storing GML Polygon in PostGIS 2, ERROR: unknown spatial reference system

texts:
I am receiving GML as part of a larger XML message. When I receive the XML it is unmarshalled to Java objects. I have no rights to alter the schemas or the generated Java code. Therefore, I later have to marshal the OpenGIS polygon back to XML. When I try to store it in the PostGIS database, I get the following exception:

  org.postgresql.util.PSQLException: ERROR: unknown spatial reference system

This is my query:

INSERT INTO test (polygon) VALUES (ST_GeomFromGML('
  <?xml version="1.0" encoding="UTF-8" standalone="yes"?>
  <ns1:Polygon srsName="http://metadata.ces.mil/mdr/ns/GSIP/crs/WGS84E_3D" ns1:id="P1" xmlns:ns6="http://www.isotc211.org/2005/gts" xmlns:ns5="http://release.niem.gov/niem/external/ogc/xls/1.1.0/dhs-gmo/2.1.0" xmlns:ns2="http://www.w3.org/1999/xlink" xmlns:ns1="http://www.opengis.net/gml/3.2" xmlns:ns4="http://www.isotc211.org/2005/gmd" xmlns:ns3="http://www.isotc211.org/2005/gco">
    <ns1:exterior>
      <ns1:LinearRing>
        <ns1:pos>32.23426 -110.34511 90.2</ns1:pos>
        <ns1:pos>32.23435 -110.34544 90.2</ns1:pos>
        <ns1:pos>32.23444 -110.34555 90.2</ns1:pos>
        <ns1:pos>32.23426 -110.34511 90.2</ns1:pos>
      </ns1:LinearRing>
    </ns1:exterior>
 </ns1:Polygon>', 4326))

If I try to remove the SRID from the query, I get the same exception.
The column I try to store the polygon in, looks like this:

AddGeometryColumn('test', 'polygon', 4326, 'POLYGON', 3);

Any ideas? Please let me know if you need any more information.


texts: Help me convert this US Lambert Conformal Conic Google Maps Projection to Australia

texts:
Following on from this question: Identify what projection this Australian map uses? <- this is the sort of map i want to load in

I need help changing the bits necessary in this example of Lambert Projection in Google Maps to make it work with Australian tiles instead: Lambert equal area azimuthal projection Google Map which apparently uses these formulae: Lambert Conformal Conic Projection -- from Wolfram MathWorld

Here is my attempts so far, changing what i think needs changing, but it must not be enough

its a bit confusing for me, im a web developer
i will be adding the tiles at a later stage
in my version i made it so when you click around it shows the coords, but as you do this the ones it shows are not in australia, so thus my problem
in the original demo it centers and loads tiles 001 & 101, in mine it does not
in my version i made it try to load all tiles so you can see which ones in the error console

Thanks for your time.

Heres is the original javascript it uses as seen at that site:
(you may notice it uses google maps api v2, not the latest v3 available due to its age)

var phi0 = 10*Math.PI/180;
var phi1 = 20*Math.PI/180;
var phi2 = 70*Math.PI/180;
var n = (Math.log(Math.cos(phi1)) - Math.log(Math.cos(phi2)))/
    (-Math.log(Math.tan(phi1/2 + Math.PI/4)) + Math.log(Math.tan(phi2/2 + Math.PI/4)));
var rho0 = (Math.cos(phi1)*Math.pow(Math.tan(phi1/2 + Math.PI/4),n))/
    (n*Math.pow(Math.tan(phi0/2 + Math.PI/4),n));

function ConformalConicLambertProjection() {};
ConformalConicLambertProjection.prototype = new GProjection();
ConformalConicLambertProjection.prototype.fromPixelToLatLng = function(pixel, z) {
    var x1 = (-0.0034375*(Math.pow(2,7 + z) - pixel.x))/Math.pow(2,z);
    var y1 = -0.0017254901960784311*Math.pow(2,1 - z)*
        (pixel.y - 405.68181818181824*Math.pow(2,-1 + z));
    var rho = Math.sqrt(2.143305952697287 + Math.pow(x1,2) -
        2.9280067982826044*y1 + Math.pow(y1,2));
    var t = Math.asin(x1/rho);
    var lng = 78.17779218926289*t-100;
    var lat = 28.64788975654*(-3.141592653589793 +
        4*Math.atan(2.004758218840869/Math.pow(rho,1.3644598756425434575)));
    return new GLatLng(lat, lng);
};
ConformalConicLambertProjection.prototype.fromLatLngToPixel = function(latLng, z) {
    var t = Math.PI/180*(100 + latLng.lng())*n;
    var rho = (Math.cos(phi1)*Math.pow(Math.tan(phi1/2. + Math.PI/4.),n))/
        (n*Math.pow(Math.tan((Math.PI/180*latLng.lat())/2. + Math.PI/4.),n));
    var x1 = rho*Math.sin(t);
    var y1 = rho0 - rho*Math.cos(t);
    x = Math.round(Math.pow(2,z)*(128 + 290.9090909090909*x1));
    y = Math.round(Math.pow(2,-1 + z)*(405.68181818181824 - 579.5454545454546*y1));
    return new GPoint(x,y);
};
ConformalConicLambertProjection.prototype.tileCheckRange = function(tileIndex, zoom, bs) {
    if(0 <= tileIndex.x && tileIndex.x < Math.pow(2,zoom) &&
       0 <= tileIndex.y && tileIndex.y < Math.pow(2,zoom-1)) {
         return true;}
    else { return false; };
};

function load() {
  if (GBrowserIsCompatible()) {

    var lambertTileLayer = new GTileLayer(
        new GCopyrightCollection(), 1, 4);
    lambertTileLayer.getTileUrl = function(tile, zoom) {
        return "Tiles/ConformalConicLambertTile_" + tile.x + "_" + tile.y + "_" + zoom + ".png";
    };
    var lambertMap = new GMapType([lambertTileLayer],
        new ConformalConicLambertProjection(), "ConformalConicLambert",
        {tileSize:256});

    var map = new GMap2(document.getElementById("map"), {mapTypes:[lambertMap]});

    // Add the controls
    map.addControl(new GSmallZoomControl());
    map.setCenter(new GLatLng(37.9251,-100), 1);

  } else {
    document.getElementById('map').style.backgroundColor = '#DDDDDD';
    document.getElementById('map').innerHTML = 'Sorry, your browser does not appear to be compatible with Google maps.';
  }
}


texts: Get Attributes from Multiple overlapping WMS layers

texts:
I am working with Openlayers to consume WMS from GeoServer. I have loaded multiple layers and get attributes for one layer using getGetFeatureInfoUrl. But this approach gets attributes for the specific layer on which getGetFeatureInfoUrl method is called.

I have overlapping layers , hence I need to get attributes from all overlapping layer at a specific point.

Below is how I am doing it :

var wmsSource = new ol.source.TileWMS({
    url: url,
    params: { 'LAYERS': 'localhost:village_maps', 'TILED': true },
    serverType: 'geoserver',
    // Countries have transparency, so do not fade tiles:
    transition: 0,
    transparent: true,
    opacity: 0.5,
    //format: new ol.format.GeoJSON()
});
map.on('singleclick', function (evt) {

      var viewResolution = (view.getResolution());
      var url = wmsSource.getGetFeatureInfoUrl(
        evt.coordinate, viewResolution, 'EPSG:3857',
        { 'INFO_FORMAT': 'application/json'}); //text/html

      if (url) {
          var parser = new ol.format.GeoJSON();
          var lookup = {};
          $.ajax({
              url: url,
              dataType: 'json',
              method: 'GET',
              success: function (response) {
                  var result = parser.readFeatures(response);
                  //if ((result.length > 0) && (wmsLayer.getVisible() == true)) {
                  if (result.length > 0)
                  {
                      var coord = evt.coordinate;
                      var content = '
';

                      for (var i = 0, ii = result.length; i < ii; ++i)
                      {
                          // Gets attributes from wmsSource only, But I have multiple other sources as well !
                          var text = result[i];
                          var props = text.getProperties();
                          var txt = text.getId();
                          content = content + 'ID:' + txt + '
';
                          content = content + 'District:' + props.district + '
';

                      }

                      content = content + '
';
                      content_element.innerHTML = content;
                      overlay.setPosition(coord);
                  }
              }
          });
      }
  });

Is it like for every layer I need to url again and then render the json ?  or is there a better way to accomplish this.


texts: Violation of foreign key constraint while copy-pasting features

texts:
I am building my first urban drainage model using GISwater, which connects QGIS and EPA-SWMM. I am using QGIS 2.18.15.

I have problems with copying and pasting elements from original conduits .shp layer to EPA SWMM DATA->Hydraulics->Arc->Conduits.

I receive next message:

Provider errors:
  PostGIS error while adding features: ERROR:  insert or update on table "inp_conduit" violates foreign key constraint "inp_conduit_arc_id_fkey"
DETAIL:  Key (arc_id)=(406474) is not present in table "arc".
CONTEXT:  SQL statement "INSERT INTO inp_conduit VALUES(NEW.arc_id,NEW.barrels,NEW.culvert,NEW.kentry,NEW.kexit,NEW.kavg,NEW.flap,NEW.q0,NEW.qmax, NEW.seepage)"
PL/pgSQL function ravbar.update_v_inp_edit_conduit() line 85 at SQL statement

Does anyone have a solution to this problem?


texts: Why does my phone's GPS suffer from "urban canyon" problems in Beijing but not New York?

texts:
In a separate question on Travel.SE I asked why my phone's GPS didn't seem to work in Beijing, China. Turns out that it was simply an "urban canyon" problem, and that the GPS did in fact work when I moved to higher and/or more open areas. What I'm wondering, is why is my phone capable of receiving GPS satellite signal in New York City, a much more dense city than Beijing especially in terms of skyscraper density downtown, but not so in Beijing? It seems if anything that it should be the other way around, or that I shouldn't be able to receive signal in either city.

Please note: I am not asking about location services on the phone, which can be influenced by things such as nearby WiFi networks and certain service carriers. I am specifically asking about GPS satellite reception coverage.

Clarification:

Outside of downtown Beijing my phone is able to receive GPS signal just fine. This question is not an issue of the GPS receiver on my phone altering the signal or any type of Chinese law. I know that the GPS receiver works perfectly. I'm simply confused as to why downtown Beijing has particular problems compared to a denser city such as downtown NYC.


texts:
I try to use QGIS to display some data stored in PostGIS, some of the schemas don't show except public; and it's a testing database, so I have only one super user: postgres. I use PostgreSQL9.3.

Anyone has some suggestion?


texts: R or ArcGis: average per polygon of the shortest distance to a river from each point in that polygon

texts:
All the shapefiles necessary can be found here. These are shapefiles of lakes, and rivers. Downloading one shapefile as an example should be enough to get me on the right track.

The nuts file can be downloaded here.

What I want is to calculate the distance from each possible point in one nuts3-polygon to the nearest river. If all the distances are calculated from each point, I want to determine per nuts3-polygon what the average, minimum and maximum distance is to the nearest river, together with the standard deviation (if possible).

My paint skills are probably not the best ever, but the idea is something like this picture below (for each of the dots). Then I would calculate the average of all the brown lines that originates from a point in the polygon.

I can work a bit with R and ArcGis. Online I do not find a solution that matches what I am looking for.


texts: I want to use QGIS with Qt 5.But i don't know the exact relationship between Qt and QGIS.Will there some problems if i use qgis with qt5 or it will not depand on it...


texts: Convert polygon with filled region inside hole into multiple polygons

texts:
Suppose you have a polygon with an inner ring representing a hole, and inside that hole is another inner ring (an "island"):

The above shape can be generated using:

SELECT
  ST_MakePolygon(
    ST_ExteriorRing(ST_Buffer(center, 10, 10)),
    ARRAY[
      ST_ExteriorRing(ST_Buffer(center, 6, 6)),
      ST_ExteriorRing(ST_Buffer(center, 3, 3))
    ]
  ) AS geom
FROM ST_MakePoint(0, 0) AS center;

While this can be represented as a single polygon, I have a requirement to convert this into multiple polygons. I am working with a library that supports polygons with holes, but does not support "islands" inside those holes - thus, I need to be able to extract those islands into separate polygons. And in case those islands also have holes that have even more "inner" islands - I need to be able to be able to extract those inner islands as well, such that no polygon in the result set contains "islands". Multiple holes (inner rings) are allowed - just no islands inside any of those holes (no inner ring wholly contained inside another inner ring).

How can I achieve this using PostGIS?

Also just to expand my knowledge... is there a more precise terminology I can use instead of "island"? I'm aware of "inner ring", but an inner ring can represent either a hole, or an "island" inside a hole. Just wondering if there is a more specific term for when an inner ring represents an island vs a hole.


texts: Error placing a static image with OpenLayers 3?

texts:
Working with the last version of OpenLayers 3, I'm not able to load correctly a static image over the map. After some test (with similar configuration) I found that the image looks fine on OL2, while on OL3 it has a displacement of several kilometers, as you can see on the following image:

I've done a plunker with that example. Could it be a bug or I'm missing something? http://plnkr.co/lzAKNGdkXkz6ZCM9H6ml

/* OL3 */
  var graphicOL3 = new ol.layer.Image({
    source: new ol.source.ImageStatic({
      url: 'https://raw.githubusercontent.com/VictorVelarde/ImagenEstuario/master/Composite.png',
      imageSize: [864, 1024],
      imageExtent: [-791350.472638389, 5376837.581117962, -778150.4726383891, 5398392.773800777]
    })
  });
  var osm = new ol.layer.Tile({
    source: new ol.source.OSM()
  });
  mapaOL3 = new ol.Map({
    layers: [osm],
    target: 'mapaOL3',
    view: new ol.View({
      center: [-784751.048, 5389384.957],
      zoom: 13
    })
  });
  mapaOL3.addLayer(graphicOL3);


texts: QGIS: disabling point deplacement when zooming out

texts:
I have a project in which I have multiple sets of data at the same coordinates, so I'm using the QGIS point-displacement renderer to display these points.  When zoomed in on the points this works great, and gives me just what I want:

But as I zoom out and the points become more tightly clustered, the algorithm generates ridiculous results:

I would like to tune this behavior so that, for example, QGIS would disable point displacement after a certain zoom level.  I looked at the rule-based renderer, but it doesn't appear to be possible to combine rule-based and point-displacement styles.

Can you suggest a way of resolving this behavior?


texts: Joining Multiple Features to One Feature

texts:
I have two layers - an index layer (polygon) and a town boundary layer (polygon).  The desired output is a FC with each town represented as only one feature and that FC's table having a field with values including all index layer features that intersect each town.

Is there one tool in ArcGIS that can accomplish this or am I looking at a multi-step process?


texts:
I need to calculate a lot of distances between several cities in Germany which are identified by a postal code, e.g. 10115 for Berlin. I cannot use the Google Maps Distance Matrix API because it would be way too expensive for the amount of requests (even with caching).

Does anybody know a self-hosted/offline API (maybe based on OSM?) where I can request a distance matrix/list of distances between cities identified by postal codes? Unfortunately, I couldn't find anything by Googling.


texts: Create shapefile with images as attributes using JavaScript

texts:
I want to make a feature in a react.js web application where I can export/download shapefiles from the dashboard. I have implemented shp-write which takes a GeoJSON input. It is working fine with textual data. I would like to add 3 attributes which should contain an image in the shapefile. I have saved the required Image in my local storage. There is no JavaScript NPM module to do so

My implementation

var shpwrite = require('shp-write');
var options = {
    folder: 'myshapes',
    types: {
        point: 'mypoints',
        polygon: 'mypolygons',
        line: 'mylines'
    }
}
// a GeoJSON input for features
shpwrite.download({
    type: 'FeatureCollection',
    features: [
        {
            type: 'Feature',
            geometry: {
                type: 'Point',
                coordinates: [0, 0]
            },
            properties: {
                name: 'Foo',
                image: 'C:\Users\DELL1\Downloads\Survey.png'
            }
        },
        {
            type: 'Feature',
            geometry: {
                type: 'Point',
                coordinates: [0, 10]
            },
            properties: {
                name: 'Bar',
                image: 'C:\Users\DELL1\Downloads\Survey.png'

            }
        }
    ]
}, options);

The Current output has image path as a text type field. It doesn't contain image.

My Expected output should contain the image itself.

What approach can I use in JavaScript to achieve this result.


texts: efficient workaround to treat missing values in gridded data in `RasterBrick` object

texts:
I have dealt with E-OBS gridded dataset and how to cropped raster grid for a particular country (I used 0.25-degree regular grid daily mean temperature observation). However, I have read about E-OBS grid dataset' publication and they told about maximum 20 percent of missing values might be occurred. After I cropped the raster grid for my interested regions and want to render plain text data with its metadata (geo-coordinate pair and daily level temperature observation), original missing values become 0.

It is important to me because it is difficult for me to differentiate which one is temperature observation (it could be 0 Celsius degree) and which one is missing observation. I want to render missing value with -999 or ....

To avoid this problem, I want to preserve missing values in a cropped raster grid, and treat those missing values when I intend to calculate the yearly average temperature for each coordinate of grid.

How can I deal with missing values in raster grid? Any idea?

Update:

I tried @Rodrigo' solution down below:

> tg1980 <- raster::brick("data/tg_0.25deg_reg_1980-1994_v17.0.nc")
> tg1980 <- reclassify(tg1980, cbind(NA, -999))
Error: cannot allocate vector of size 1.9 Gb

But now I got a memory problem in R.  Any quick solution for that?