texts:
I have a simple code here some map with zoom+pan effect. It work fine on Opera, Chrome even IE. But Firefox do not show it.-

If you open "Inspector" panel, you'll see has same values (translation, scale, path, d**) on Firefox than the other, but don't render properly.

Obviously is Firefox's bug cause on other browser works nice and D3 generate the code correctly.-

And here same map without zoom+pan effects, Firefox show it nice. Some google do not gave me any reference, console operation give the right value.

Any idea what's going wrong?

Here an snapshot. Can you see path.county's size?

(*) code take it from Mike's site: http://bl.ocks.org/mbostock/5914438
He use:

var projection = d3.geo.mercator()
    .scale((1 << 12) / 2 / Math.PI)
    .translate([width / 2, height / 2]);

and

var zoom = d3.behavior.zoom()
    .scale(projection.scale() * 2 * Math.PI)
    .scaleExtent([1 << 11, 1 << 14])
    .translate([width - center[0], height - center[1]])
    .on("zoom", zoomed);

(**)

<svg width="600px" height="500px" transform="translate(5579.288888888889,-2536.2020481280824)scale(32768)"><path class="county" d="M239.07661387245844,54.......


texts: Trilateration using 3 latitude/longitude points, and 3 distances?

texts:
I want to find out an unknown target location (latitude and longitude co-ordinates). There are 3 known points (latitude and longitude co-ordinate pairs) and for each point a distance in kilometers to the target location. How can I calculate the co-ordinates of the target location?

For example, say I have the following data points

37.418436,-121.963477   0.265710701754km
37.417243,-121.961889   0.234592423446km
37.418692,-121.960194   0.0548954278262km

What I'd like is what is the math for a function that takes that as input and returns 37.417959,-121.961954 as output.

I understand how to calculate the distance between two points, from http://www.movable-type.co.uk/scripts/latlong.html I understand the general principle that with three circles like these you get exactly one point of overlap. What I'm hazy on is the math needed to calculate that point with this input.


texts: Why does QGIS 2.4 crash at every shutdown?

texts:
I am using QGIS on two Win7, 64 bit machines for years. After upgrading to QGIS 2.4 it crashes on both machines at every shutdown with the message "Minidump...". The crashes occur independant of what I have done before or whatever project I am working on. Until now there seems to be no loss of data, but it is annoying. (I do not have installed the LecoS-plugin on either machine.)

Any ideas, solutions or similar problems?


texts: Stop OpenLayers featureselected event listener blocking map drag event

texts:
When trying to drag the map, if the user clicks on a vector feature the featureselected event is triggered and the map dragging does not happen. Please try to drag the map starting from any feature in the example below.

http://openlayers.org/dev/examples/select-feature-multilayer.html

Is it possible to disable the featureselected event so the map can be dragged?

Solution

i found this FAQ http://trac.osgeo.org/openlayers/wiki/SelectFeatureControlMapDragIssues

It works.


texts:
I'm trying to build a mosaic, for my XYZ Tile Server, from a VRT file that map about 9000 tiffs of satellite images that i previously warped to EPSG:4326 and compressed. Unfortunately the mosaic created appears to have only transparent tiles. These are the steps:
// 1. Warp all the tiffs to EPSG:4326
gdalwarp -t_srs EPSG:4326 -r average -dstalpha %mypath%!infile! %mypath%!outfile!

//2. Build the VRT with all those tiffs listed in a txt file
gdalbuildvrt -overwrite C:\tiffsList.vrt --optfile C:\tiffsList-warped-compressed.txt

//3. Build the mosaic from the VRT
gdal2tiles.py -z 1-12 -w all C:\tiffsList.vrt C:\xyz-tiles

What am I doing wrong?


texts: Extract Point Cloud bounding box with PDAL in Python?

texts:
Data:

I have a point cloud file, let call it "pc_file" (it's not a .las but a .ply file ; https://en.wikipedia.org/wiki/PLY_(file_format) ).

Goal to achieve:

I am searching a fast way to extract the bounding box in the xy plane.

To do that, from the working directory, I can:

A)

Run pdal info --boundary pc_file > boundary.json
Parse the json file through a python script where I need the bbox:

import os
from shapely import wkt

working_dir = '/absolute/path/to/my/WorkingFolder/'
boundaryFile = os.path.join(working_dir, 'boundary.json')

with open(boundaryFile) as f:
    data = json.load(f)

wkt_boundary1 = data[u'boundary'][u'boundary'].encode('ascii','replace')
bbox1 = wkt.loads(wkt_boundary1).bounds

This is clean, but the main drawback is that I need two different commands; one bash line followed by a run of the python script (of course I can stack them both in a bash file but I will end up with a temporary 'boundary' file written on the disk, which I do not want if possible). That was why I was searching to do that all at once without any disk writings:

B)
 1. Run an other python script where I directly call the pdal module ( https://pdal.io/python.html ):

    import os
    import numpy as np
    from io import StringIO
    import pandas as pd
    import pdal
    from shapely import wkt

    working_dir = '/absolute/path/to/my/WorkingFolder/'
    PCFile = os.path.join(working_dir, 'pc_file')

    json_str = u"""
    {
      "pipeline": [
        """+'"'+PCFile+'"'+""",
        {
            "type": "filters.hexbin",
            "threshold":10
        }
      ]
    }"""

    print("json string: {}".format(json_str))
    pipeline = pdal.Pipeline(json_str)
    pipeline.validate()
    pipeline.loglevel = 8 #really noisy
    count = pipeline.execute()
    arrays = pipeline.arrays
    metadata = pipeline.metadata
    log = pipeline.log

    json_obj = pd.read_json(StringIO(metadata))['metadata']['filters.hexbin']
    wkt_boundary  = json_obj[1][u'boundary'].encode('ascii','replace')
    bbox = wkt.loads(wkt_boundary).bounds

But this second solution takes really longer as it computes a close envelope-like boundary ( https://pdal.io/stages/filters.hexbin.html ):

The question:

I did not found a way to call a "pdal info" like command directly in python from the pdal module (which would be, in my humble opinion, a clean way to extract point cloud metadata from python).
The pdal module seems to only deal with "pipeline" structures, thus calling "filters".
The pdal "info" command does not exist as a filter ( https://pdal.io/stages/filters.html ).

Is there any or some filter to extract the bounding box directly?

(Don't hesitate to tell me if the question would have better chances on SO directly.)


texts: Using Interactive Overlay to Export Geographic Coordinates

texts:
I can barely spell JavaScript, but I have envisioned an application working with or embedded in Google Earth that employs the following work flow:

This application is intended on working at scales of 1–10 km.
Apply a user-created overlay image in the Earth Viewer window, aligned to the local geographic grid, that shows detailed elevation data of the underlying terrain using chromatic scales.
Insert a user-created transparent overlay containing a polygon or other geometric feature on top of the first overlay.
Position, x-y deform, and/or rotate the geometric overlay so that its figure conforms to an underlying terrain feature revealed by the elevation overlay. Any of these alterations can be easily performed on Google Earth overlays.
Export to a database geographic data relating to the geometric overlay, such as its center coordinates, the rotation angle commanded, and/or x-y deformations performed on the geometric overlay for further computations and analyses.

Potential computations using such data include determining a feature's orientation to True North, its perimeter measure, area, aspect ratio (length ÷ width), eccentricity/asymmetry, and other descriptive information.

Could one create such a tool using the Google Earth Engine environment available? And if so, is there a ready resource for getting started on such an application? I haven't been able to make much headway with what is available at the GEE project.


texts: Find Maximum value per month in rasterstack/brick and corresponding timestamp R

texts:
When dealing with raster images, I found the raster and rts packages in R to be very flexible for analysing my raster data unlike ArcGIS. In my learning efforts, I have had to do some stacking/bricking of individual images to form a continuous time series.

Now, I would like to extract some values of interest from the stacks just as one will do with a time series object in R.

My raster data (A) is rainfall every 30mins. So layers in my rasterstack are separated by30minstime step. For comparison with another raster data (B) set which is at 6-hourly time step, I summedA` to 6-hourly as follows:

dates1=seq(as.POSIXct("2015-04-01 00:00:00"), as.POSIXct("2015-11-30 23:59:59"), by="6 hours",tz="GMT")

index=rep(seq(1,976,by=1),each=6) #11712/12=976 nlayers/# of intervals for accumulation

ras <- setZ(A,dates1)
sixhrly <- stackApply(ras,indices=index, fun=sum)

Output from a pixel of sixhrly looks like:

         `time`          `value`

    2015-11-28 00:00:00 5.452091e-02
    2015-11-28 06:00:00 4.030154e-01
    2015-11-28 12:00:00 3.107421e-01
    2015-11-28 18:00:00 8.851760e-02
    2015-11-29 00:00:00 9.171816e-02
    2015-11-29 06:00:00 8.803389e-02
    2015-11-29 12:00:00 3.286928e-02
    2015-11-29 18:00:00 0.000000e+00
    2015-11-30 00:00:00 5.635915e-02
    2015-11-30 06:00:00 1.299293e-03
    2015-11-30 12:00:00 3.354174e-02
    2015-11-30 18:00:00 6.680167e-03

1) using stackApply(which I find to have the most speed on my system), how can I find the maximum value per pixel month? The output is a rasterstack, with layer1=April, layer2=May etc...

2) How can I obtain the timestamps for each max value? This is important if I want to know what time the max rainfall occurred and if A and B.

A more generic way for doing this is most welcome as well since data can be on hourly and minutes time step over many years and one has to extract some values of interest.


texts: Changing land use land cover classes in ArcGIS Desktop?

texts:
I have a tiff file containing raster data of land use land cover (LULC) classes.

Each pixel has a pixel value in the form of RGB values of the corresponding LULC class, need to change the pixel values of certain regions on the map and make a new map with revised LULC classification.

How can I do that in ArcGIS Desktop?


texts: Changing extent of raster in loop using R

texts:
I have raster objects of a certain extent (0, 360, -90, 90) and I want to change that extent for all of them to a different extent (-180, 180, -90, 90).
The thing is that I have my objects (raster layers) in a list and I want to change the extent with one step in a loop.
My approach (example)
st_ext <- extent(-180, 180, -90, 90)

  for (i in 1:NROW(r)){
    r[[i]] <- setExtent(r[[i]], st_ext)
  }

When I do it like this, the extent does not change but when I use the exact same command where I put in "i" manually, the extent changes like it should.
r[[1]] <- setExtent(r[[1]], st_ext)

Why does this happen and how I can make it work with the loop?
I also tried
extent(r[[i]]) <- st_ext

and
xmin(r[[i]]) <- -180
xmax(r[[i]]) <- 180

But when using all of those those in a loop, the extent remains unchanged.


texts: Has anyone successfully installed HAZUS MR4?

texts:
I am attempting to install HAZUS-MH MR4 (v1.4) on a Windows XP SP3 computer with ArcGIS 9.3 SP1 installed.  The installation seems to be running fine up until it attempts to start services.  Then it displays this error message:

I realize it's a bit of a long shot, but has anyone else installed HAZUS, recieved this error and found a successful solution for it?  I've pretty much exhausted the solutions on the FEMA site, and am hoping someone may have a pointer in the right direction.


texts:
I have a tiff file and I'd like to convert it to mbtiles using GDAL 2.2.2, compiled with sqlite3.
The first command seems to work :

gdal_translate -of mbtiles 000-2017-0915-6528-L93-0M20-E100.tif test.mbtiles
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.

Now when I try to compute overviews with gdaladdo (like in doc):

gdaladdo -r average test.mbtiles 2 4 8

ERROR 4:test.mbtiles' not recognized as a supported file format.`

However MBTILES is in the list of supported format in gdaladdo --formats.

Does anyone have an idea ?


texts: Getting map by id or class using OpenLayers?

texts:
How can I get hold of a OpenLayers3 map object associated with a specific id or class (unique)?

I searched and found this thread, which says something like:

var map = $('.my_map_class').data('openlayers'); //assuming there is just one map on the page

But I got an typeerror saying that map is not defined. I have tried replacing $ with jQuery, and it didn't work either.

Any pointers?


texts:
I'm trying to carry out Tasseled Cap Transformations in GRASS GIS. I have a Landsat 8 image with bands 2-7 that are needed for the calculation. I have imported the bands into GRASS and attempted to carry out the calculation but when I do, it outputs an image size of 1x1. I have checked all of the input bands and they are all 8202 x 8113 (rows x columns).

The screenshot below is the output of r.info on the output file. I think it is doing the calculation on the whole map rather than on a pixel by pixel basis, but I've no idea why. Does anyone know how to fix this?

I'm using Fedora 27 and GRASS GIS 7.2.2. I've checked online but no one else seems to be having the same problem.

r.in.gdal input=LC08_L1TP_208023_20170408_20170414_01_T1_B2_2157.TIF output=B2

r.in.gdal input=LC08_L1TP_208023_20170408_20170414_01_T1_B3_2157.TIF output=B3

r.in.gdal input=LC08_L1TP_208023_20170408_20170414_01_T1_B4_2157.TIF output=B4

r.in.gdal input=LC08_L1TP_208023_20170408_20170414_01_T1_B5_2157.TIF output=B5

r.in.gdal input=LC08_L1TP_208023_20170408_20170414_01_T1_B6_2157.TIF output=B6

r.in.gdal input=LC08_L1TP_208023_20170408_20170414_01_T1_B7_2157.TIF output=B7

i.tasscap sensor=landsat8_oli input=B2,B3,B4,B5,B6,B7 output=tct_test


texts:
I have a script which scans a directory and outputs basic raster data information such as the file name, format, number of bands, and etc. I need a way to make it so if the directory does not contain raster data (i.e., anything other than raster data), a message is displayed stating that the directory doesn't have the correct data type.

I know ArcPy has a Describe() function that I could use to determine the type of data in a folder, but am not sure how to implement it. This is what I have so far:

rasterList = arcpy.ListRasters("*", "ALL")
filesType = arcpy.DataType('RasterDataset') # Can use `DatasetType` as well.
                                            # I've tested this function to describe
                                            # raster data and ArcPy prints out
                                            # 'RasterDataset', that is why I have it
                                            # there in the brackets.
for name in rasterList:
    if rasterList == filesType:
        print ("\nFilename:"), name
    else:
        print ("This directory does not contain any raster data.")

Any suggestions?


texts:
So, I generated a map on tilemill and exported the mbtiles file. I've setup Tilestache, and all required dependencies. I tried the OSM preview provided by default which worked fine.

When I tried to load a MBtiles file through the Tilestache config file(.cfg), the terminal doesn't show me any errors (infact it shows me all the HTTP calls perfectly fine), but the preview output at 127.0.0.1:portno/{layername}/preview.html seems to be a blank map!

When I did explicit HTTP calls to specific png images, the images are empty. In other words, I see empty outputs from it. I also tried different MBtiles files and I'm still facing the same output.

What could be the problem here??

My Tilestache config file:

{
  "cache":
  {
    "name": "Test"
  },
  "layers":
  {
    "roads":
    {
        "provider": {"name": "mbtiles", "tileset": "ward.mbtiles"},
        "preview": {"zoom": 15, "ext": "png"}

    }
  }
}


texts: How to automatically set an adequate zoom level for tiled data from stamen in R

texts:
I have a project where I need to automatically download data from OSM/Stamen in R. So I also need to be able to automatically set an adequate zoom level for the tiles to be downloaded. It should show sufficient detail given the area of interest. My idea was e.g. to use the extent of the bounding box to set a good level i.e. to define some breaks for the 20 or so levels from OSM which have described "qualitative" breaks on wikipedia, see here. So I wondered if there is already something done in this way or if I have to write a function from scratch.

Here is a minimal reproducible example with a manual choice of an adequate zoom level

library("ggmap")

# manually choose zoom level
zoomlevel<-12
# see all zoom levels here: https://wiki.openstreetmap.org/wiki/Zoom_levels

map <- get_stamenmap(bbox = c( -75.18343,5.382292,-74.95132,5.592376), zoom = zoomlevel, maptype = "terrain")

mymap<-
  ggmap(map) +
  theme_minimal()

mymap


texts: How to identify in googlemaps with below data


texts:
Is there max length limitation on WKT name?

see the bold text below:

GEOGCS["GCS_Wake_Eniwetok_1960",
    DATUM["Wake_Eniwetok_1960",
        SPHEROID["Hough_1960",6378270,297]],
    PRIMEM["Greenwich",0],
    UNIT["Degree",0.017453292519943295],
    AUTHORITY["EPSG","37229"]]

Can you please point me to link or official documentation for WKT specs limitations or length?


texts: Iterate calculate field in multiple tables in ModelBuilder

texts:
I would like to iterate calculate field in multiple tables using ArcGIS Model Builder.
The figure below explains details of the calculations. My existing model calculates value only for the first table.

Any direction would be valuable.


texts: QGIS Buffer Not Recognizing All Features in Layer

texts:
I'm completing online course GST 101 for QGIS (my version is 2.18.3) and I'm stuck on the Buffer tool.

All of my layers are set to CRS "NAD83(HARN)/New Mexico Central (ftUS) (EPSG:2903)" as instructed in the lab.

The goal is to use the query builder to filter the municipal limits in Bernalillo county down to only Albuquerque. Then after filtering create a buffer layer of 5280 ft around it.

The problem is I get an empty layer from the buffer tool with no entries in the attribute table. So I tried Buffer on the entire layer (with 0 ft buffer to pinpoint the problem areas) and this is the result.

There appears to be a number of features in this layer that aren't recognized because I would have expected the entire layer to be selected.

I've also tried the query builder to filter Albuquerque and save to a new layer. Then try to buffer the new layer, still results in an empty attribute table.

I've tried the shp file included in the lab as well as downloading the file directly from here:
City of Albuquerque GIS Data (Municipal Limits / Various municipal limits in Bernalillo county / SHP File).

I'm still fairly new at GIS but I'm following a youtube video and pdf instruction but I'm clearly on a newer version of QGIS.

Any suggestions on what the problems may be?


texts: Identify polygon overlaps within single shapefile

texts:
Using either QGIS (1.9) or ArcInfo (9.3.1) I'm trying to identify unique overlaps of prehistoric building footprints from an archaeological site recorded within a single polygon shapefile. While the Overlay>Interset and Geoprocessing>Intersect tools in ArcInfo and QGIS, respectively, work to some degree both return multiple results for each unique overlap. An appropriate parallel for my problem I suppose can be illustrated by a simple 3 unit Venn diagramme in which there are 3 unique intersections or overlaps, yet when running an Intersect routine (particularly that in ArcInfo) I often get up to 12 unique polygons returned - 2 for each shared 'leaf-shaped' intersection, and a further 6 for the centre overlap (each intersection being counted twice).

I'm trying to produce a shapefile illustrating those unique intersections. Ideally I would love for the attribute table of this produced shp to contain details of each intersection i.e. Feature_X intersects with Feature_Y, but I fear I may be asking for the moon on a stick.

Many thanks in advance for any and all advice,
Marc

Edit: Following the excellent suggestion of @PolyGeo, I have attached a rough sketch below to hopefully better illustrate my current problems and what I'm hoping to achieve. Forgive the hand-drawn sketch, I'm away from my main PC.


texts: Clustering geometry points, return a point's label if the point does not belong to a cluster?

texts:
I am clustering geometry points using the k-mean algorithm. The result is a list of clusters. Some clusters contain multiple points. Others only contain a single point. On the client side, we would like to display markers on a map. Clusters containing multiple points will be displayed as a marker containing an integer equivalent to the number of points in that cluster. Clusters containing a single point will be displayed as a marker containing the point's name.

With SQL and postGIS functions, I would like to fetch clusters and the name of any point that is the only point in a cluster.
I am nearly there. As shown in the attached image, the following SQL query returns a list of clusters, their kmean index, the number of points in the cluster (called the count) and their centroid point. This is all great, but for all clusters that contain only a single point, I'd also like to include the name of the point.

SELECT kmean, count(*), ST_Centroid(ST_Collect(location_point)) AS centroid,
       CASE
           WHEN count(*) > 1 THEN 'cluster'
           ELSE 'voter_name should go here'
       END "member"
FROM (
    SELECT st_clusterkmeans(location_point, 3) OVER () as kmean,
           voter_name, location_point
    FROM voter_locations
     )
     as ksub
group by kmean
order by kmean

This is the voter_locations table where the geometry points are stored.

If, in the query, I replace the text 'voter_name should go here' with the voter_name variable, then I get an error

ERROR:  column "ksub.voter_name" must appear in the GROUP BY clause or be used in an aggregate function
LINE 3:            ELSE voter_name

However, if I do add the voter_name variable to the group by clause, then the response includes one row for every record, like in the following photo...

 This isn't good because we want one row for every cluster...

How can I tweak my query to achieve the desired results?


texts: Using ArcGIS Desktop field calculator (python) to parse date

texts:
I'm trying to parse a datetime field using python (but have never used python for it before).  I have a field that contains a valid date, but not a valid time... and other fields that have valid hour and am/pm information.  I need to grab the date (mm/dd/yy) from the datetime field and then combine that with the valid hour and am/pm fields into a new string field.

Below is the vbscript way that I would do this (and it works fine, but I need to build it using python):

    (DatePart ( "m", [STARTDATE])) &"/"&
    (DatePart ("d", [STARTDATE] )) &"/"&
    (DatePart ("yyyy", [STARTDATE] )) &" until "&
    (DatePart ( "m", [ENDDATE])) &"/"&
    (DatePart ("d", [ENDDATE] )) &"/"&
    (DatePart ("yyyy", [ENDDATE] )) &" (from "&
    [StartHour] &" "&
    [StartAMPM] &" until "&
    [EndHour] &" "&
    [EndAMPM]&" )"

Here's what I've written so far using python (but doesn't work):

    def dateFull(oid)
    startDateObj = datetime.datetime.strptime( !STARTDATE! , "%m/%d/%Y")
    endDateObj = datetime.datetime.strptime( !ENDDATE! , "%m/%d/%Y")
    return startDateObj + " to " +  endDateObj + " from " + !StartHour! + " " + !StartAMPM! + " to " + !EndHour! + " " + !EndAMPM!

Can someone please help me figure out what I'm doing wrong here?


texts: GeoServer SLD Filter rule before PointStacker transformation

texts:
I'm trying to apply a filter rule in an SLD, before applying the GeoServer PointStacker transformation (which groups points that are too close to distinguish their icons). This transformation is nicely described and demonstrated here and here etc.

Without my filter, the transformation is working fine, and without the transformation, my filter is working fine. My filter must go in the SLD instead of in CQL because it allows certain points to be conditionally displayed at certain zoom levels, depending on an attribute field. I can't seem to find where and how to put in the filter condition without breaking the XML validation (or even a related question on SE).

Below is what I have tried already. I've put the filter in after the transformation, so as to show what I'm trying to achieve. (Of course this does not work for two reasons: 1\ the point-stacker transformation does not return my original fields used in my filter, so the filter no longer works there. 2\ it stacks points which should have been filtered out already by the time the transformation is done.)

I've also tried moving the <ogc:filter> tag and contents to just before the <Transformation>, as well as just inside it - at the start or the end. All these attempts gave XML validation errors.

EDIT: removed fluff from the code here and formatted. The code is exactly as per the first link in par. 1 above, but with a filter added.

<?xml version="1.0" encoding="UTF-8"?>
<StyledLayerDescriptor version="1.0.0" xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd" xmlns="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <NamedLayer>
    <UserStyle>
      <FeatureTypeStyle>

        <Transformation>    <!-- exactly as per in the references above -->
        </Transformation>

        <Rule>
          <ogc:Filter>      <!-- I want this filter to be applied before the Transformation  -->
          </ogc:Filter>

          <PointSymbolizer>
          </PointSymbolizer>
          <TextSymbolizer>
          </TextSymbolizer>

        </Rule>

      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>

GeoServer v2.7.2 with OpenLayers 3.9.0. The WPS extension is installed (required for the transformation).


texts: Bing's geocoder with GeoPy not returning a zipcode

texts:
I'm just getting started with GeoPy and have been running some tests
to determine the quality of the different available geocoders. I have
noticed that Bing's geocoder is more accurate than Google's (for my
purposes, at least), but instead of returning a zip code, Bing returns
<built-in function zip> rather than the zipcode.  Any thoughts as to
why/ways to correct this?

Here's my test code:

from geopy import geocoders
g = geocoders.Google()
b = geocoders.Bing('<DEV KEY HERE>')
services = [g, b]

def search_list(addylist):
    for address in addylist:
        for x in services:
            for place, (lat, lng) in x.geocode(address,
exactly_one=False):
                print "%s: %s: %.9f, %.9f \n" % (x, place, lat, lng)

addy1 = '513 N 36th Ave , Seattle'
addy2 = '121 S First Street , La Conner'
addy3 = '854 NE 58th St, Seattle'
addylist = [addy1, addy2, addy3]

search_list(addylist)

My output is:

<geopy.geocoders.google.Google object at 0x18e75d0>: 513 36th Ave,
Seattle, WA 98122, USA: 47.606811900, -122.287247300

<geopy.geocoders.google.Google object at 0x18e75d0>: 513 N 36th St,
Seattle, WA 98103, USA: 47.651385400, -122.351573800

<geopy.geocoders.bing.Bing object at 0x18e76d0>: 513 N 36TH St,
Seattle, WA <built-in function zip>, United States: 47.651340485,
-122.351692200

<geopy.geocoders.google.Google object at 0x18e75d0>: 121 1st St, La
Conner, WA 98257, USA: 48.391576600, -122.496045800

<geopy.geocoders.google.Google object at 0x18e75d0>: 121 N 1st St, La
Conner, WA 98257, USA: 48.394719800, -122.494525700

<geopy.geocoders.bing.Bing object at 0x18e76d0>: 121 S 1ST St, La
Conner, WA <built-in function zip>, United States: 48.391513824,
-122.495391846

<geopy.geocoders.google.Google object at 0x18e75d0>: 854 NE 58th St,
Seattle, WA 98105, USA: 47.670977000, -122.317815000

<geopy.geocoders.bing.Bing object at 0x18e76d0>: 854 NE 58TH St,
Seattle, WA <built-in function zip>, United States: 47.670978546,
-122.317817688

Bing seems to be more consistent, only returning a single correct result rather than two results.  The only issue is that Bing doesn't seem to want to return a zipcode, for some reason, rather it returns a string containing the address and <built-in function zip>.


texts: Symbolizing vector features using RGB values stored in three fields using ArcPy?

texts:
I need to symbolize a vector layer in Arcmap based on RGB (or HSB - can do either) values found in fields in the attribute table. I'm working with up to 10 different datasets, each with over 100 features. As a result, I can't manually adjust each feature's fill attribute, nor can I use a template .lyr file and then apply its symbology.

Is there a method via Arcpy that I can use to perform this?

Looking at the documentation it appears that symbology isn't exposed to Arcpy but I'm not familiar enough to know for sure.

Otherwise, are there any third-party scripts or tools that can accomplish this?


texts: GeoServer GetLegendGraphic request ServiceException

texts:
I have a layer in GeoServer which is stored in a PostGIS database. The geometrytype of the layer is GEOMETRY and now the table only contains LINESTRING. But I what to add points and polygons to the table, this is the reason I choose GEOMETRY as the geometrytype.
My layer is working fine in the Layer Preview...

But when I do a GetLegendGraphic request I get the following ServiceException:
java.lang.IllegalArgumentException: class com.vividsolutions.jts.geom.Geometry is not supported by this method
class com.vividsolutions.jts.geom.Geometry is not supported by this method

My SLD is quite simple:

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0" xmlns="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc"
  xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
  xsi:schemaLocation="http://www.opengis.net/sld http://schemas.opengis.net/sld/1.0.0/StyledLayerDescriptor.xsd">
  <NamedLayer>
    <Name></Name>
    <UserStyle>
      <Abstract></Abstract>
      <FeatureTypeStyle>
  <Rule>
  <Title>Test</Title>
  <ogc:Filter>
   <ogc:PropertyIsEqualTo>
    <ogc:PropertyName>feature_type</ogc:PropertyName>
    <ogc:Literal>Test</ogc:Literal>
   </ogc:PropertyIsEqualTo>
  </ogc:Filter>
          <LineSymbolizer>
      <Stroke>
     <CssParameter name="stroke">#FFA618</CssParameter>
     <CssParameter name="stroke-width">2</CssParameter>
     <CssParameter name="stroke-dasharray">5 2</CssParameter>
      </Stroke>
    </LineSymbolizer>
  </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>

My GetLegendGraphic request is:
http://localhost:8080/geoserver/wms?REQUEST=GetLegendGraphic&VERSION=1.0.0&FORMAT=image/png&WIDTH=16&HEIGHT=16&LAYER=namespace:layer

I have tried the request on a GeoServer 2.0.1 and a GeoServer 2.0.2 and in both cases I get ServiceException as above mentioned.

What can I do to get the request working?


texts: ArcGIS Label Manager 10.2 takes 20-30 seconds to load

texts:
Why does label manager take 20-30 seconds to load? I am on 10.2 with a SSD and 24gb of RAM.

Does anyone else have issue loading the label manager in ArcMap? Is there any way to speed this process up?

It almost seems like the longer ArcMap has been installed, the longer this process takes.

I just loaded a shapefile with 8 simple polyline features and it took 20+ seconds to load the label manager.

I should also note that this has been hapening for me since 10.1. I also did a fresh install of my OS and 10.2 a couple months ago


texts:
I want to convert the coordinates (639117.49, 6865775.92) from EPSG:4326 to EPSG: 3949, by using the following code:

from pyproj import Proj, transform

P4326 = Proj(init='epsg:4326')
P3949 = Proj(init='epsg:3949')

x,y = transform(P4326,P3949,  639117.49, 6865775.92)

The values should be: 1639190.94757331, 8187869.9312972 but I get inf, inf.
What is the mistake?
Otherwise, can I do it directly via PostGIS without using Python?


texts:
I'm rather new to QGIS and I have the following issue:

I have a raster image with 255 different categories all having their own colour. However, I want to cut out some of this data as it is not relevant to me. When going to properties > style > paletted, I'm basically searching for a way to delete some of the categories I do not need. Does anyone know how I can do this?


texts: Crash after performing random forest using enmap box and dzetsaka

texts: i want to perform random forest classification using tools like dzetsaka and enmap box on qgis 3.2 version but after classification completes upto 50 percent qgis crashes what can be the issue


texts:
I just upgraded to QGIS 3.0 (OS X as well on Windows 10) and can not find out how to do some georeferencing. I am told by the software, that GDAL georeferencing is still installed and that I could not remove it, and that it should be in the category "raster". But I can not find it there, also not in other parts of the menu.

new profile created, but the same issue...:


texts: GDAL Python API - statistics valid percent information

texts:
I would like to know if the STATISTICS_VALID_PERCENT (percentage of valid (not nodata) pixels) information is accessible through the GDAL Python API?

It seems to be missing from the GetStatistics Function.

As a work around, I have been scraping this information from the PAMDataset (XML file) generated through the GetStatistics(False, True) call.

Is there a better way?


texts: Clip polygon layer with another polygon layer - KeyError: 1006

texts:
I have two different polygons layer : the first one represent buildings and the second one represent cities.

My goal is to extract the buildings which are in the cities.

I tried to use the Clip function but it returns :

Uncaught error while executing algorithm
        Traceback (most recent call last):
          File "/usr/share/qgis/python/plugins/processing/core/GeoAlgorithm.py", line 203, in execute
            self.processAlgorithm(progress)
          File "/usr/share/qgis/python/plugins/processing/algs/qgis/Clip.py", line 71, in processAlgorithm
            source_layer.crs())
          File "/usr/share/qgis/python/plugins/processing/core/outputs.py", line 326, in getVectorWriter
            crs, options)
          File "/usr/share/qgis/python/plugins/processing/tools/vector.py", line 618, in __init__
            uri = GEOM_TYPE_MAP[geometryType] + "?uuid=" + unicode(uuid.uuid4())
        KeyError: 1006

I have been able to use this function to extract points (after having created a point layer with the centroïds of the buildings) and then keep only buildings which have the same IDs as the centroids extracted but I had to do it with PostGIS and it is quite long to do.

Is there any way to do it with a simplier method ?


texts: Drawing lines showing connections between geocoordiante points

texts: I have an ArcMap document with 100 locations plotted using their latitude and longitude values. I want to show connections between the 100 points in the graph. Essentially 100*100 = 10,000 lines. Is there a way to construct a graph showing these connections using ArcMap?


texts: Visualizing queue length with colored Polylines - Folium

texts:
I have a bunch of PolyLines that i want to color code based on the length.
I managed to change the color of the entire line, based on the length of a line segment (between 2 of 6 points), but I am having trouble figuring out how to color parts of the line instead of the entire line.

I am to visualize queue lengths, so I want to change the color of each line/queue based on the length, so if shorter than 5 m = green, 5-10 m = yellow, >10 m = red.

What I've done so far is:

def color_length(x):
    if x < 0:
        raise ValueError
    elif 0 < x < 5:
        return "green"
    elif 5 < x < 10:
        return "yellow"
    else:
        return "red"

And then I put it into:

folium.PolyLine(points, 'color_length(x)', weight=2.5, opacity=1).add_to(mymap)

But as mentioned, when doing this, it returns only one color for entire line, and not for the segment, which the "def" refers to.


texts: Issues authenticating with Copernicus Open Access API

texts:
I'm having issues with the Copernicus Open Access API when it comes to authentication.

If I put e.g. https://scihub.copernicus.eu/dhus/odata/v1/Products('2b17b57d-fff4-4645-b539-91f305c27c69')" in my browser I get a HTTP Basic login form and when I enter my data I get a valid response (i.e. I'm authenticated successfully).

The documentation states that I can also use curl. But when I execute curl -u {username}:{password} "https://scihub.copernicus.eu/dhus/odata/v1/Products('2b17b57d-fff4-4645-b539-91f305c27c69')" I get a 401 Not authorized error. The same problem persists with the Java library I'm using to query results.

Now if I explicitly set the authorization header in curl like curl -H "Authorization: Basic {base64encodedUsernameAndPassword}" "https://scihub.copernicus.eu/dhus/odata/v1/Products('2b17b57d-fff4-4645-b539-91f305c27c69')" everything works but this approach seems like a hack and I cannot apply it to the Java library without forking.

Is there something I'm missing here? Why it doesn't accept the -u parameters?


texts:
I am using ArcGIS 10 (includes spatial and geostatistical analyst).

I have a raster (GRID) with a value range of -167 to 18.

I would like to query the raster and create a new one with only the values between -117 and -69.

I cannot figure out how to use the raster calculator to do this task.

Any suggestions?


texts:
I'm trying to calculate WDVI index but for some reason the computation failes and I don't find my mistake. The function works but all the pixels are 0.

var addWDVI = function(image) {
  var red=image.select('B4');
  var RedEdge=image.select('B8A');
  var WDVI = RedEdge.subtract(red.multiply(RedEdge.divide(red)))
  .rename('WDVI')
  .copyProperties(image,['system:time_start']);
  return image.addBands(WDVI);

};

This is the link to my code-
https://code.earthengine.google.com/43e76e0665ceb3dc20dcbc8a8cae5e25

My end goal- to add this index to my code


texts: Georeferencing to x y coordinates with one lat long origin?

texts:
We are running parallel transect lines for a survey, which will roughly cover a rectangular area. When making the transects, we will periodically capture spatial locations of features and events we encounter.

Due to technical issues, we can only record the latitude/longitude coordinates of the very beginning of the first transect (i.e. one corner of the rectangle). The locations of all our logging of other features/events will be recorded relative to the origin, with x y coordinates (running along the east-west north-south axes, respectively).

I would like to georeference all our data points from the survey into a map, where I can click on any data point and get its coordinates in latitude and longitude.

How should I do this?

I prefer a solution that could be done in QGIS 1.7+, but ArcGIS 10 is OK, too.


texts:
I have a polygon layer and a table without geometries. I join the table to the polygons to use one of the table attributes to label the polygons. I don't use the table for anything else.
Is there a way to completely hide the table from the user? I don't want it to show in layer tree:


texts: Previewing pages of field group section in report layout

texts:
I used the report layout in QGIS 3.0 to re-create two atlases I made using 2.18.  I'm using these two tutorials from North Road and Anita Graser.

Atlases seem to work the same way in the regular layout in 3.0 as they did in the print composer of 2.18. But the report layout works a little differently, and I'm still trying to figure it out.

I added two field group sections to my report, and set up each section as the atlas was set up in my 2.18 project. When I exported the report, it had all the pages I expected, just as if I had created atlases in two different layouts.

Now I want to preview each page of the report and make sure it looks okay before exporting the report again. In 2.18 and in the layout of 3.0, there's an Atlas toolbar that let you preview the atlas and flip through the pages. I don't see any way to use the Atlas toolbar in a report, though.

How do I preview the pages of a field group section in a report in QGIS 3.0?


texts: GDAL 2.0.1/QGIS 2.14.0 on OSX won't open MODIS HDF?

texts:
I'm on OSX 10.9.4, using GDAL 2.0.1 from the terminal. I want to mosaic MCD12Q1 HDFs together with GDAL but my build can't even read them.
The file: http://e4ftl01.cr.usgs.gov/MOTA/MCD12Q1.051/2001.01.01/MCD12Q1.A2001001.h00v08.051.2014287161513.hdf
Trying to gdalinfo that file:
gdalinfo failed - unable to open 'MCD12Q1.A2001001.h20v13.051.2014287163801.hdf'.

gdalinfo --formats shows that my build supports:
      HDF5 -raster- (ros): Hierarchical Data Format Release 5
      HDF5Image -raster- (ro): HDF5 Dataset

So I assume these MODIS HDFs are using the HDF4 standard, which isn't listed in my formats list, and that's why GDAL has no idea what to do with them.
My QGIS Essen 2.14.0 says it's running against GDAL 1.11.4, which can't read them either; I add the HDF and it asks me to pick a raster layer in the file, but when I select one nothing appears and the error log says "Unable to load."
I saw this question about the same issue on Windows, but apparently HDF4 support is an optional flag when you install GDAL? Anyone else have this issue with MODIS HDFs on OSX? Is there a way to get HDF4 support into my GDAL builds? I installed both GDAL and QGIS using KyngChaos's installers.
EDIT: QGIS would not load these HDF4 files because their directory path had quote marks " in it. Moving them to a directory without quotes fixed this.


texts: Can't replace attribute with field calculator

texts:
I have an attribute table where one of the columns (laeskaerm) has been edited/created by checkboxes on 2 different machines.

Result: it contains both 1, 0, true, false and NULL

I'd like to replace 1 with true, and 0 with false. I've tried:

replace("laeskaerm",'1','true')

and

replace("laeskaerm",'1','true')

But I can't seem to make it work (QGIS 2.2.0).

Do anyone have any idea whats wrong?


texts: How to add data to relationship table in ArcGIS 10.1 database design?

texts:
I have created 2 tables HORSE and VENDOR and a relationship table VENDOR_HORSE
With the foreign key in Horse being the Primary key in VENDOR (VendorID)
and the foreign key in VENDOR being the primary key in HORSE (HorseID).
It is a many to many relationship and ArcGIS automatically creates the link table -.

yet I had to create the Relationship table which is simple and has extra fields (besides the horseID and vendorID) as Notes (text), Service (coded domain),

and Date (most important to solve many to many matter and have a relational database.
I cannot add data to the VENDOR_HORSE table.
I click on start an edit session, just as I do to add rows to the HORSE and VENDOR tables, with no problems, yet I cannot add rows to the Relationship table.
ESRI help gave instructions for adding rows yet they do not work.
 That is,

I selected a row in HORSE,
I selected a row in VENDOR,
then double clicked on OBJECTID to have the table fields appear in
the Attribute window
and NOTHING. What am I missing?
Help!
Thanks
Bette


texts: Getting log error about default locale settings when using GRASS with QGIS?

texts:
I downloaded python 3.6, installed it, dowloaded QGIS 3.0.2, installed GDAL, then QGIS, then opened QGIS and tried to use GRASS processing tools. Initially I got the following error:

env: python3: No such file or directory
env: python3: No such file or directory
Execution completed in 0.50 seconds
Results:
{'lake': <QgsProcessingOutputLayerDefinition {'sink':/Users/yalngarpiril/Desktop/MAP/water/ocean/ocean.tif, 'createOptions': {'fileEncoding': 'System'}}>}

I then followed the instructions to add a custom variable on kyngchaos.com. I quit and restarted QGIS and opened the same GRASS processing tool and got a different error:

WARNING: Default locale settings are missing. GRASS running with C locale.
Starting GRASS GIS...
ERROR: Path '/<UNKNOWN>/<UNKNOWN>' doesn't exist
Default locale not found, using UTF-8
Exiting...
WARNING: Default locale settings are missing. GRASS running with C locale.
Starting GRASS GIS...
ERROR: Path '/<UNKNOWN>/<UNKNOWN>' doesn't exist
Default locale not found, using UTF-8
Exiting...
Execution completed in 0.88 seconds
Results:
{'lake': <QgsProcessingOutputLayerDefinition {'sink':/Users/yalngarpiril/Desktop/MAP/water/ocean/ocean.tif, 'createOptions': {'fileEncoding': 'System'}}>}

I'm on a 2010 MacBook Pro running macOS Sierra 10.13.4.

What causes my current error?


texts: Changing page layout from portrait to landscape using ArcPy with ArcGIS Pro?

texts:
Changing page layout from portrait to landscape using ArcPy with ArcMap? is answered by saying that changing a layout size in ArcMap requires you to use .Net or comtypes rather than ArcPy (arcpy.mapping).

With ArcGIS Pro now available, can it be done there using ArcPy (arcpy.mp)?


texts:
I work with this perimeter\area fields:

i created it with "Add geometry columns" vector tool. How can i create a dynamic perimeter\ area fields, so when i create a new feature or change existed feature- those fields will be change also.

I try to seek the solution by create a new area field and try to find the correct expression in the field calculation, but with no success.

I can use again the "Add geometry columns" vector tool, so the attribute table will be update- but this solution is very clumsy.


texts: Change Measure Tool units through ArcObjects

texts:
I'm activating the measure tool through code by using the ControlsMapMeasureToolClass.
I want programmmatically to change the default length units of tool.
How do  do that?


texts: ArcGIS Dissolve failing with small polygons

texts:
I have an arcpy script which clips (Clip_analysis) a shapefile of polygons to a grid square, converts any multipart polygons created to singlepart polygons (MultipartToSinglepart_management), runs dissolve (Dissolve_management) and then gets the area of the dissolved polygon to get the proportion of the grid square that is covered by polygons.

In the arcpy, the clipping is done using a 400m x 400m rectangular arcpy.Polygon, created from a point array.  When replicated in arcGIS desktop, I used a polygon shapefile with one feature (the same rectangle as used in the arcpy).

This script has been run many times without issue, but now I get an error in a particular (400m x 400m) grid square:

Error 999999: Error executing function.
Invalid Topology [Y coord limit exceeded.]
Failed to execute (Dissolve)

I can replicate the error using shapefiles/geodatabases and ArcGIS tools (10.1).  The grid square is nowhere near the extents of the coordinate system (the grid square is in London and the Coordinate system is British National Grid (EPSG: 27700)), and the grid squares North, East, South and West of that square all work fine.

The grid square only contains two small triangles from the same original polygon, and neither extend further than 75cm from the base of the grid square.  If I manually edit the polygon to extend over 1m from the base of the grid square, dissolve works fine.

I've tried repairing the geometry, but that didn't help.  I have tried various XY Resolutions and Tolerances (including not set, 1cm and 1m) with no success

What could be causing this error?  Is there a work-around?

When I run using the tool dialog, I get the same error, so it looks like a Dissolve tool or data problem, rather than a python problem.


texts: Increase maximum zoom level for a WebTiledLayer in ArcGIS

texts: I am using the WebTiledLayer from the Javascript ArcGIS API to render a custom WMTS source. The WMTS I am using has zoom levels going from 18 to 21, however I cannot zoom in further than what appears to be zoom level 19. Does anyone know how to increase the maximum zoom level for a WebTiledLayer?


texts: Error while creating spatial database connection

texts:
I'm using ArcGIS 10.0, and I've already created my database using "ArcSDE for PostgreSQL Post installation" ..

My database is successfully created in PostgreSQL and now that I'm trying to create a spatial database connection from ArcMap, so I went to :

ArcCatalog ==> Database connections ==> Add spatial Database connection..

I've filled all the textboxes with the needed information and when I test the connexion before I submit, I get this error :

  Failed to connect to the specified server.

  Entry for SDE instance not found in services file.

Although I think I've provided the right information as below :

  Server : localhost

  Service : sde:postgres:localhost

  Database : my_project


texts: Reprojecting elevation dataset for "tilted Earth"

texts:
I would like to create maps of a "tilted Earth", such as the following example where Australia is at the North Pole:

http://i1048.photobucket.com/albums/s379/laskaris_mods/tilted_earth.jpg

I can reproject maps in image formats like .jpeg like this with the software I am using, Geocart 3. However, what I need to do is reproject elevation datasets in netcdf format (the ETOPO .grd data) in just this way, so that I have exact elevation data for my "tilted Earth".

Is there a GIS software that can do this?


texts: Using more than one group for single layer in GeoServer?

texts:
I have a huge number of layers in my geoserver that I have painstakingly imported and they all work correctly when viewed with GeoExplorer etc.

I want to know if there is any way to create set maps which show groups of these layers. I know how to group the layers together but don't know if i can use the same layer in multiple groups.

If this is entirely impossible, what is the best way to create set maps- using a new workspace/store for the layers that require duplicating?

The maps will eventually be linked into a server that hopefully will have different interactive maps created from the above.


texts: Export properties of a feature collection Google Earth Engine

texts:
I would like to export only the properties of a filtered feature collection, I made a try which is probably wrong because I'm not selecting the properties I want (area, FinalDate,Id,InitialDate ) and it generates a very big csv file. So the question is how to select only the properties I need?
var dataset = ee.FeatureCollection('JRC/GWIS/GlobFire/v2/FinalPerimeters');

var mask = ee.FeatureCollection('RESOLVE/ECOREGIONS/2017')
                              .filter(ee.Filter.eq('ECO_ID', 569));

// Filter the table geographically
var filtered = dataset.filterBounds(mask);
print('First ', filtered.first());

var visParams = {
  palette: ['f5ff64', 'b5ffb4', 'beeaff', 'ffc0e8', '8e8dff', 'adadad'],
  min: 0.0,
  max: 600000000.0,
  opacity: 0.8,
};

// Visualization
//
var image = ee.Image().float().paint(filtered, 'area');
Map.centerObject(mask, 5);

Map.addLayer(image, visParams, 'GlobFire Amazonas');
Map.addLayer(filtered, null, 'for Inspector', false);

// Export the FeatureCollection to a csv file.
Export.table.toDrive({
  collection: filtered,
  description:'GlobFireChacoSeco',
  fileFormat: 'csv'
});


texts: Creating new slot in RasterLayer using R?

texts:
I am trying to use a function from a package in R that requires a slot in my raster/stack with time information. Ideally, this simple code should do the trick:

lu_2012 <- raster("f2012lcs_new.rst")
lu_2012@t <- 2012

But it doesn't, since this is not the right way to create slots:

lu_2012@t <- 2012
Error in (function (cl, name, valueClass)  :
  ‘t’ is not a slot in class “RasterLayer”

Is there a way to create a new slot in a RasterLayer/stack like that?


texts:
I'm a beginner SQL and GIS user and this my first time using stackexchange, so apologies if I'm asking the question properly. I'm trying to open some MSSQL tables in QGIS. I established a connection without any problems, but I am having trouble adding layers to QGIS. Some tables will be added properly, while others will not. I keep getting the error message The Layer dbName = 'x' host = x, estimatedmetada = true, srid = 0, table = "dbo"."practicetable". "practicetable" sql=is not a valid layerand can not be added to this map. I think the problem has something to with the lack of metadata. When I look at properties, I see a whole bunch of stuff in metadata for tables that work, while it only says error in metadata in tables that don't load.

Here's how the table I'm most interested in was created.

create table practicetable
(name varchar(10),
longitude float(5),
latitude float(5),
srid varchar(50))
insert into practicetable
values('home' , 39.78, -86.14, 'EPSG:4326 - WGS 84')

When I cut and paste the output into a csv, I am able to add it as a delimited text layer without problems.


texts:
I have managed to create a 3D-landscape from DEMs and are trying to get my building layer in 2.5D in this landscape. My buildings are retrieved as polygons with height defined in the attributes. When I choose height, the buildings get what seems to be an unlimited height in the map. When this is manually changed to a value of "0.0001-0.0003" it suddenly starts to resemble some realistic proportions. But I want these buildings to represent an actual scaled representation and not an subjective representation of the height.

The building vector layers height attributes are represented in meters and have been retrieved from OpenStreetMap @ QGIS 3.12.1. I can also say that the same problem occurs when I am creating new vector layers for buildings.

Is there a way to change some sort of a "reference scale"?

Any tips on how this might be fixed? I'm still learning, so even "the wrong solution" will most likely give me more knowledge in the field of cartography.


texts: Are CartoDB Map views by month or by year?

texts:
We have 10000 map views with the free version of CartoDB.

Are they reset every month or is it the maximum views for the all utilisation ?


texts: What is the most efficient way to search a geodatabase for NULL-like records?

texts:
The problem: I have a geodatabase with several datasets and many more feature classes within.  The fields within the feature classes have been populated through joins with shapefiles and manual edits.  Often times string fields will become populated with whitespace (i.e. '', ' ', '   ', etc) or the string "Null", and numeric fields will become populated with a zero (0).  I would like to find these records and replace them with a true NULL value.  I have the following code, which uses an UpdateCursor, but it still very slow and doesn't catch all of the NULL-like records.  Does anyone know of other ways to accomplish this task?

GDB = arcpy.GetParameterAsText(0) #input geodatabase
arcpy.env.workspace = GDB
datasetList = arcpy.ListDatasets() #list datasets

for dataset in datasetList:
        arcpy.env.workspace = os.path.join(GDB, dataset)
        fcList = arcpy.ListFeatureClasses()
        for fc in fcList:
                arcpy.AddMessage("Processing %s..." % fc)
                #count features
                arcpy.MakeTableView_management(fc, "tempTableView")
                count = int(arcpy.GetCount_management("tempTableView").getOutput(0))
                if count > 0:
                        fieldList = arcpy.ListFields(fc)
                        for field in fieldList:
                                arcpy.AddMessage("...%s" % field.name)
                                rows = arcpy.UpdateCursor(fc)
                                for row in rows:
                                        count = 0
                                        if row.isNull(field.name):
                                                continue # if already null
                                        elif field.type == "Text":
                                                value = row.getValue(field.name)
                                                if value.lstrip(' ') == '' or value.lower() == '<null>': # looks for whitespace or '<null>'
                                                        row.setNull(field.name)
                                                        count += 1
                                        elif field.type == "ShortInteger" or field.type == "LongInteger" or field.type == "Float" or field.type == "Double":
                                                value = row.getValue(field.name)
                                                if value == 0:
                                                        row.setNull(field.name)
                                                        count += 1
                                        if count > 0: # update row if records have changed
                                                rows.updateRow(row)
                                del rows
                else:
                        arcpy.AddMessage("...NO RECORDS FOUND.")


texts: How to change icon symbols in QGIS? Can icons have JavaScript animations?

texts:
There are a lot of SVG files located in the "QGIS Dufour" Program Files folder on my PC. I was thinking that I could find the symbol that I want to change in the GUI in this listing in InkScape so I can actually have my own custom data points. However, when I transfer the edited .svg file back to the folder it's contained in I get an error: "Error 0x80070570"

I don't know what these data are, but this is the data associated with the symbol I want to replace:

<!DOCTYPE qgis_style (View Source for full doctype...)>
- <qgis_stle version="1">
 - <symbols>
  - <symbol alpha="1" type="marker" name="star3">
.
.
.

It goes on. I'm not sure what this is.

I just want a biohazard symbol and a school symbol.

Just as an aside, can these images have JavaScript animations?


texts: Clipping with results going to existing feature class using arcpy.Clip_analysis()?

texts:
I am using arcpy.Clip_analysis() in arcpy immediate window to clip feature. This function creates new object as a result of clipping. That take more processing time. Can I do the clipping and save the clipped object to an existing feature class?

My sample code:

arcpy.Clip_Analysis("Road", "clipRectangle", "D:\output.gdb\road2")

Now, some of the suggestions below is to use "in_memory" to clip faster by writing to the RAM instead of hard disk. I tried the following:

arcpy.Clip_Analysis("Road", "clipRectangle", "in_memory\road2")

That works fine, but now I will have again to write them in the hard disk to have them saved. Wont that be the same as if I clip directly to hard disk as first code?

Is that the only way or there are some other approaches that can do faster?


texts: Add CAD feature class to Arcgis using arcpy

texts:
I'm trying to add a CAD layer (point, polyline, annotation...etc) into ArcMap by python code, i'm using Arcgis 10.1, IDLE(Python 2.7),
I faced many codes but they all gives errors, so i'm blocked here !

I tried this :

arcpy.AddCADFields_conversion("C:/Saadia/SAFTOP/Parcelle et rivrains/simulation arcgis/python-script python/Progr_work/plan.dwg","ADD_ENTITY_PROPERTIES")

And also this :

arcgisscripting gb = arcgisscripting.create()
gb.ImportCAD_conversion("'C:/plan.dwg'", "C:/workspace")

But no one works !

I tried this :

import arcpy
arcpy.env.workspace = "D:/Users/saadiya/Documents/ArcGIS/Default.gdb"
arcpy.env.overwriteOutput = True

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
addPoint = arcpy.mapping.Layer(r"C:\plan.dwg\Point") # reference to point layer
addPolyline = arcpy.mapping.Layer(r"C:\plan.dwg\Polyline") # reference to Polyline layer
arcpy.mapping.AddLayer(df, addPoint, "BOTTOM")
arcpy.mapping.AddLayer(df, addPolyline, "BOTTOM")

But i got this error :

  RuntimeError: Objet\xa0: CreateObject can not open the map

I changed it to this :

mxd = arcpy.mapping.MapDocument(r"C:\Saadia\BD\Projet_Marrakech.mxd")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
addPoint = arcpy.mapping.Layer(r"C:\plan.dwg\Point") # reference to point layer
addPolyline = arcpy.mapping.Layer(r"C:\plan.dwg\Polyline") # reference to Polyline layer
arcpy.mapping.AddLayer(df, addPoint, "BOTTOM")
arcpy.mapping.AddLayer(df, addPolyline, "BOTTOM")

it fixed this error :

df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
IndexError: list index out of range


texts:
In QGIS I am working on a system to collect sewer data. At the moment I have like 25 manholes investigated. The information that I need are Material, shape and diameter of the sewer pipes.

At the moment I have three value maps:

Material:

Concrete
PVC
Spirosol
Stoneware
Asbestos Cement

Shape:

Round (Concrete)
Egg shaped (Concrete)
Square (Concrete)
Round (PVC)
Etc….

Diameter:

250mm (Round concrete)
300mm (Round concrete)
750/750mm (Square Concrete)
Etc.

What I would like is that when I select concrete, only the shape options for concrete (Round/Egg shaped/Square) become available. When I then select for instance square, only the diameter for concrete square becomes available.

What I am looking for is some sort of If than statement in order to make it easier to register data.

Can somebody tell me how to do this?


texts: Convert ArcMap label expressions from VBScript to Javascript or Python

texts:
I have a large number of layers (hundreds) with label expressions written in VBScript.
I am using ArcMap 9.3.1, and I'm now starting to export the layers as MPKs for use in our ArcGIS Runtime (WPF) application.

I get an error that says that the layers can not be exported to MPK with Runtime support because it has label expressions written in VBScript.

It will take to long, and too much manpower to manualy convert the label expressions from VBscript to Python or Javascript. I'm looking for a tool to do this automatically.

Another possible solution can be to enable the export of layers with VBScript expressions to MPK format. This solution would be the best case scenario!

Any help will be very much appriciated!


texts: How do I combine raster data in different coordinate systems and different cell size into one data set with fine cell size with ArcMap?

texts:
input:

CORINE land use data (cell size 100m x 100m, LAEA coordinates, EU countries) and GLC ("Global Land Cover") data (1km x 1km, WGS84, Europe, land use classes changed to match CORINE data with reclassify tool)

output:

one land use (discrete values) raster data set with 100m x 100m cells in LAEA coordinates.

details:

For EU countries I want the data to remain exactly as it is, with LAEA coordinates and 100m cell size. For non-EU European countries I would like to have the land use data from GLC but in an expanded 100m x 100m raster that is a seamless continuation of the CORINE raster. (I know that the coarser information does not normally justify just splitting one cell up into 100 and giving all of them the same value as the original larger cell. However, for my purposes, this approach makes sense and I really want to do it this way.)

How can I do this with ArcMap? What tools do I need to combine in what order?
Project GLC data first, then merge?

When I tried to project the GLC data, I always had to have unspecified raster size in order for it to work (10^-3? To account for differently oriented traster cells, I guess). That is not what I want at all. For my purposes, I really need a raster with equally sized "squares". Outside of EU, the information does not need to be absolutely accurate, if some of it gets distorted a bit, that's okay, I just really need a raster with squares. Also, I need to make sure that CORINE data stays exactly the way it is in EU countries.


texts:
For a web app I want to source roadmap and satellite type maps from a web API which will be used with geometrical data from a database to generate reports. The reports would be viewed on screen or printed out.  Print would require 300 dpi resolution to avoid blurring. Width/height approximately 17cm would be ideal.

Sources of map data may be Google Maps, MapBox, OpenStreetMap or another. OpenStreetMap has fair usage policies and their servers aren't designed for commercial use.

I would like to find out what Server Side Rendering options exist.


texts:
Cartography question.

I do have administrative boundries for counties in polygon FC. What tool do I use to prioritise display of features. In example, when lines overlap- I want to show only 1 dashed line instead 2 overlapping lines with shifted dashes.

Similary, when I do have 2 separate polygon feature classes, I want to prioritise certain class over another in way, that higher rank is displayed and lower is not. Of course in both cases rules would affect only overlapping segments of the polygon or polyline, if it would be more convinient.

[EDIT]
On the picture overlapping lines in area 1 should look like line in the area 2.


texts: Is GPS more accurate on specific hours of the day?

texts:
I have heard that the GPS system is more accurate, and less accurate on some specific hours on the day.

If this is true, how can I get the information on when GPS is more or less accurate at my location? And how big is the difference?

I guess this depends on where in the world you are, I'm living very north, so it's maybe has a bigger effect here. I also think it has to do with how many GPS satellites you get signals from at the moment.


texts:
I have problems setting up a clustered Geoserver. So far got four instances (version 2.2.4) running on a ubuntu linux box with a shared data directory behind HAProxy.
I use one of the geoservers to configure the cluster (master) and reload configuration in the other instances (slaves). In practice I often get errors when reloading configuration on the slaves. Only restarting each instance in the cluster distributed the configuration to all instances. It seems that i am doing something wrong.
Has anyone made some experiences in clustering geoserver?

EDIT:

Thank you Brad and sorry for getting back so late. Here I provide you with some details:
the cluster was set up according to http://www.tomcatexpert.com/blog/2010/07/20/basic-tomcat-clustering-grails-applications
After installation a directory /var/lib/geoserver_data was added and in the startup scripts I set CATALINA_OPTS DGEOSERVER_DATA_DIR=/var/lib/geoserver_data -DGWC_DISKQUOTA_DISABLED=true -DGWC_METASTORE_DISABLED=true

So all instances share the same data directory.

So far so good. Until additional data is added:
For adding data only one of the geoserver instances is used. Adding - lets say a postgis layer - works fine. I can see the layer in the preview and use it as expected.
But when i go to the admin interface in of the other instances, hitting 'reload configuration' an error message greets me, saying:

  java.lang.IllegalArgumentException: Settings already exist for
  workspace 'abc'

Removing workspace 'abc' leads to

  "Settings already exist for workspace 'def'"

and so on.


texts: labeling point symbols with leader lines

texts: Is it possible to label the "edge" of a point symbol in ArcMap when using leader lines? I would like to offset the leader so the line does not overlap the point symbology. In the past I have converted my points to polygons to leverage anchors, but want to know if there is another cartographic method to addressing this issue. If have included a picture displaying the overlap I want to remove.


texts:
I have a GeoExt form and want to do file upload here. But it never works with GeoEXT. I have this:

buttons: [{
text: 'Загрузить',
handler: function(){
if(fp.getForm().isValid() && Ext.getCmp('opinions').getValue()=="blabla"){
form_action=1;
fp.getForm().submit({
         url: '<portlet: actionURL>',
         waitMsg: 'WAITING...',
 success: function(fp, o) {
          DO SOMETHING
          }
         }
        }
      }]

This is correct?

What arguments use in success: function(fp, o)? First argument is FormPanel but I don't understand what is second. Another question how to do asynchronous file upload in this case?


texts: DefineProjection_management - "Object: Error in executing tool'

texts:
I included
coordinateSystem = arcpy.SpatialReference(4283) and it works when I point the workspace path to a GHB file. However I get an error below when including multiple GHB folders into the workspace z:/temp.

name 'inData' is not defined.

My code is below. Can somebody suggest a fix?

import arcpy
from arcpy import env
import os

arcpy.env.workspace = "Z:/temp"
fcs = arcpy.ListFeatureClasses('*', "ALL")

for file in fcs:
   print arcpy.env.workspace, file

try:
    for file in fcs:
      inData = file
     coordinateSystem = arcpy.SpatialReference(4283)

 # updates/overwrites coordinate system information and map projection stored with a dataset
# accepts dataset and coordinate system information
      arcpy.DefineProjection_management(inData, coordinateSystem)

except arcpy.ExecuteError:
   print arcpy.GetMessages(2)
   arcpy.AddError(arcpy.GetMessages(2))
except Exception as e:
   print e.args[0]
   arcpy.AddError(e.args[0])

The ghb folders and shape files are stored in their folders as per the screenshots.


texts: Managing changing historical location name with relational database/ PostGIS

texts:
I work with income data linked to 19th century U.S. town names, and I have to join the income data with town points (year 2014) shapefiles. Many of the towns changed name for several times, I manually search for the current location of historical town names, and want to record such info in PostGIS.

I thought about using synonyms for a historical town names and its current town name in PostgreSQL, but it seems that PostgreSQL does not have CREATE SYNONYM statement....

Another way to do this might be changing the historical town names into its current town names in my income table, and then join with the 2014 town points shapefile, but that is my last resort since I want to keep the original town names.

Is there any other way to work around this issue?


texts:
In ArcGIS, when I want to Label my layer, there are a number of label styles I can choose from.

I would like to have that Feature in QGIS as well. For example, when I label my layer in a certain font, color, format and buffer, I would like to save my self-made label and just add it to another layer. Would that be possible?


texts: Why is the Turf Javascript library's 'transformTranslate' function distorting geometries?

texts:
I'm working on a webapp using MapBoxGL and TurfJS and I've encountered an issue using Turf's transformTranslate, which is defined here.

Now, the issue is that a geometry that is being transformed is getting distorted.

For example, this is the original geometry:

And this is the translated geometry:

And the relevant section of code is this:

  console.log(typology)

  var typologyCentroid = centroid(typology)
  var mapCenter = this.$store.state.map.getCenter()
  var mapCenterToTurf = lngLatToTurfPoint(mapCenter)

  var typologyDistToMapCenter = rhumbDistance(typologyCentroid, mapCenterToTurf)
  var typologyBearingToMapCenter = rhumbBearing(typologyCentroid, mapCenterToTurf)

  var translatedPoly = transformTranslate(typology, typologyDistToMapCenter, typologyBearingToMapCenter)

  console.log(translatedPoly)

The actual goal is to get a feature and shift it to be in the centre of the current map display. I can't understand why the translate is twisting and distorting the feature though.

Edit: it turns out that I can achieve what I want by just using a simple shift (adding deltaX and deltaY to each coordinate). I'm surprised that the Turf function doesn't work as expected, though... I'm leaving the question up.


texts: How to read postgres sql Table (No geometry) in PYQGIS?

texts:
In my database there is a table without geometry. how to read that from from pyqgis.
self.uri = QgsDataSourceUri()
self.uri.setConnection(self.Hostname, self.Portid, self.Dbname, self.Username, self.Password)
db = QSqlDatabase.addDatabase("QPSQL")
db.setHostName(self.uri.host())
db.setDatabaseName(self.uri.database())
db.setPort(int(self.uri.port()))
db.setUserName(self.uri.username())
db.setPassword(self.uri.password())
self.dbresult = db.open()
self.uri.setDataSource(self.Schema, "project_ids", "geom", "")
lay = QgsVectorLayer(self.uri.uri(True), "project_ids", "postgres")

print("Layer fields Names",lay.fields().names())

is shows not showing any fields How to read the No geometry Tables in pyqgis?


texts: I recently upgraded to QGIS v2.0.1 and noticed that the EasyPrint plugin is not supported for the latest version.  Will this plugin be updated to work in this version, or are there any alternatives for sheet mapping?


texts: Opening old QGIS model in current QGIS versions

texts:
There is an older answer by @underdark to this post: Finding centrelines from polygons in QGIS?
She links a QGIS model that apparently was created in QGIS 2.12. It has the file extension .model
When I try to open the model in QGIS 3.16, it appears that only models in the new .model3 file extension can be opened. Other models do not appear in the opening dialog and thus can't be selected. I tried to simply change the file extension, but that did not work. When trying to open it, I get an error: model could not be loaded.
Is there a way to convert the model to run it in newer QGIS versions?


texts: Symbology/points not showing on ArcGIS Online web map but data is present

texts:
I have a layer that has been symbolized in an ArcDesktop project.  The project was saved and Shared As --> Publish Service up to our ArcGIS Online account.  All layers and their symbology have uploaded flawlessly up to this point.  However, this particular layer does not.  The data is there and if you know where one of the points is located and click on it, a pop up window with the correct data appears.  The symbology shows up in the Content bar on the left of the web map interface just fine.

My first attempt was to upload the project with all layers.  The first 4 layers showed symbology correctly but the remaining 11 layers that all have the same symbology (which is different from the first 4 layers) are the ones not showing.  (link to web map: http://arcg.is/1S5iWtg)

Several attempts after that were to upload in various ways that would keep the symbology.  Including, uploading smaller groups, uploading individual layers at a time, smaller number of natural breaks (from 7 to 4), using more simplistic symbols, none of which worked.  While in ArcGIS Online, if I try to change the symbols, as soon as I click on the Change Style icon, all the points then appear as small black dots (which was not the specified symbology).

Here is the example of the web map with just one of the layers uploaded that are not showing up.  (link to web map: http://arcg.is/1S5ltDE)


texts: Exporting Map into PDF Distort Grey Colour

texts: Sometimes when I export large maps into a PDF format from ArcMap 10.3, the PDF doesn't retain the colour. I can print maps off ArcMap 10.3 and the maps are perfectly fine. What I mean by that statement is what I see on the screen and the printed versions are the exact same. However, when I export the map into a PDF, the colour isn't the same as what is seen on ArcMap. I've realized that PDF distorts only the layers which uses grey colour; for example, parcel lots and road line network. I am working with only vector layers. I've spent hours playing with the export setting and still haven't been able to rectify the problem. Please assist.


texts: Arcmap Python add-in wizard run code on start-up/open document

texts:
I have the following code that I want to run each time a Arcmap is open:

Windows 7
ArcMap 10.3.1

I have followed the instructions here: http://resources.arcgis.com/en/help/main/10.1/index.html#//014p00000018000000

It's suppose to run on an OpenDocument event. However the code does not run. If I copy paste into the python window in Arcmap it works fine, any ideas?

import arcpy
import pythonaddins
import getpass

class OnOpen(object):
    """Implementation for RunOnOpen_addin.OnStart (Extension)"""
    def __init__(self):
        # For performance considerations, please remove all unused methods in this class.
        self.enabled = True
        OnOpen.openDocument(self)
    def openDocument(self):
        mxd = arcpy.mapping.MapDocument("CURRENT")
        text_elements = arcpy.mapping.ListLayoutElements(mxd,"TEXT_ELEMENT")
        username = getpass.getuser()
        initial = username[0]
        lastname = username[1:]
        name = initial.title() + ". " + lastname.title()

        for element in text_elements:
            if element.name == 'Producer':
                element.text = "Producer: " + name

OnOpen()


texts: Using Arcpy to copy features in a loop to a differing workspace

texts:
I've tried to create a script to loop through a folder and then clip the features. This works. My problem is that my new files are created in the same workspace (in my case Data) but I would like to save them in a different workspace (a geodatabase).

# -*- coding: utf-8 -*-

import arcpy
import os
from arcpy import env

env.overwriteOutput = True

env.workspace = r"C:\Users\Victor\Carto\Data"

clip = r"C:\Users\Victor\Carto\XXX.shp"

list = arcpy.ListFeatureClasses()
print list

for l_list in list:
 filename = os.path.splitext(l_list)[0]
 clip_filename = filename + "_GL"
 arcpy.Clip_analysis(l_list, clip, clip_filename, 0)


texts: Lambert_Conformal_Conic: converting from shapefile to geojson fails

texts:
I have a shapefile that also contains a proj file that looks like so:

PROJCS
[
    "MGI_Austria_Lambert",
    GEOGCS
    [
        "GCS_MGI",
        DATUM
        [
            "D_MGI",
            SPHEROID["Bessel_1841",6377397.155,299.1528128]
        ],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]
    ],
    PROJECTION["Lambert_Conformal_Conic"],
    PARAMETER["False_Easting",400000.0],
    PARAMETER["False_Northing",400000.0],
    PARAMETER["Central_Meridian",13.33333333333333],
    PARAMETER["Standard_Parallel_1",46.0],
    PARAMETER["Standard_Parallel_2",49.0],
    PARAMETER["Latitude_Of_Origin",47.5],
    UNIT["Meter",1.0]
]

I want to convert this shapefile into a geojson like so:

ogr2ogr -f GeoJSON \
    -s_srs STATISTIK_AUSTRIA_STADTREGIONEN_3MIO.prj \
    -t_srs EPSG:4326 \
     STATISTIK_AUSTRIA_STADTREGIONEN_3MIO.geojson \
     STATISTIK_AUSTRIA_STADTREGIONEN_3MIO.shp

This fails with the following error message:

ERROR 6: No translation for Lambert_Conformal_Conic to PROJ.4 format is known.

Changing the proj file to contain either PROJECTION["Lambert_Conformal_Conic_2SP"] or PROJECTION["Lambert_Conformal_Conic_1SP"] instead of Lambert_Conformal_Conic results in valid geojson files that seem to have the correct coordinates BUT the resulting coordinates differ in their 2nd or 3rd decimal places, e.g.,

LON               , LAT
16.526053978650904, 48.32876392574797

vs.

16.527168091786308, 48.32902755664574

My questions are:

Which one of the two (if any) is the correct result?
Is there a intuitive explanation for the difference between Lambert_Conformal_Conic, ~_1SP, and ~_2SP that also a non-cartographer is able to understand.
There is an inaccuracy in the value of the GEOGCS UNIT value in the proj displayed above compared to the value given for the respective projection at spatialreference.org. Is this the reason for the differences?
The publisher of the dataset states that there is a more accurate (but non-free) version (1:500000 vs. 1:3000000) of the data. May it be that the inaccuracies due to the data source is bigger than the inaccuracies due to inaccurate re-projection?


texts: R 'overlay' two rasters to create one single band raster and preserve cell values

texts:
I have two single band rasters r1 and r2.

and

Put together in a GIS they look like this:

Raster extend and cell size of r1 and r2 match perfectly. The white areas are nodata cells. There is no overlapping of cells that have a cell value.

I would like to 'overlay' the two rasters to create one single band raster!

I already tried:

r_overlay <- overlay(r1, r2, fun=sum)

and

r_mean <- mean(r1, r2)

and

r_brick <- brick(r1, r2)

but none of the approaches seems to work. The goal is to create a single band raster containing the values of both rasters r1 and r2 and nodata values where there are nodata values now.


texts: QGIS Two Raster Data Sets with black line between the two sets where they join?

texts: Where the two Digital Raster Quad sets (topo maps 7.5 Quad) overlap there is a significant thick black line (not a single line but wide and is wider the further south you go on the map).  I am looking to remove this line or join the two so the map is not interrupted.  The line also goes through the middle of the area I am working.  Directions on how to solve this would be greatly appreciated!


texts: Dynamically generated map overlay with map server - what's the best option?

texts:
I'm facing a problem and would like to ask your help finding the best online mapping solution!

Project details:

The User can pick 1 from several hundred territory based functions to map its values. A good example: Population.
Every function represents data on a group of territorial units. These units can be very large (i.e countries or continents) and can be very small (settlements). Every unit is represented as a (multi)polygon.
The user can filter on what territorial units he wants to display. In other words he can filter the areas to be displayed. Based on the values of the selected function and the selected areas, the map to be created will be a 30% transparent Choropleth map overlay.  This adds a dynamic behaviour, so static maps are not an option.
The User is allowed to select huge numbers of territories. I.e: all settlements of a country. This introduces several thousands of polygons to display.
I have high resolution ESRI shapefile data on each unit. I can convert them to any GIS format.

The current implementation: uses Google Maps with Fusion Tables layer added. The Fusion Table is generated dynamically with regard to the following procedure:

A pre uploaded base fusion table is used, which holds all geometry data for all possible areas. An areaId primary key is part of the table
A dynamically generated fusion table is uploaded runtime which holds the values and area ids for the selected function.
When this dynamic table is created, we can get the fusion table to display by joining the dynamic table with the base table .
Colors are calculated runtime and applied as conditional formatting in the view layer

The problem: When you generate a new fusion table google needs some time to cache the generated map tiles on server side. So far i didn't find a solution to get feedback on these map tiles are already generated or not when displaying the map to the User. When displaying the map for the first time(and most probably the last time..since they are generated runtime) the user won't see the Choropleth map("Data may still be loading" message on map tiles), he has to refresh the page several times to get the map displayed correctly. And altogether this solution is very slow. Takes 10-25sec to upload data to Fusion Tables service and another 2-20sec to get the tiles cached on google's servers.

What i have tried:

KML Layers on Google Maps: wouldn't take KMLs larger than 10M. The caching issue is still present
Creating polygons on client side: is an overkill for the browser for large shapes like continents and for the cases when 3000-4000 smaller shapes like settlements are displayed.

Looking for: a mapping server solution that i can install and use on the project's server to dynamically generate these map tiles. The solution would qualify if the whole process wouldn't take more than 10sec. The raster layer behind the Choropleth map should be preferably Google Maps but not necessarily.

Any suggestions are most appreciated and thanks for your help in advance!


texts:
I need to do some sort of spatial clustering where each spatial coordinate has a value attribute. So, for example, each point has a salary attached to it.

So each cluster would have a total salary if I summed each point's salary. But I want to have each cluster to be between a certain total salary. For example, New York can spatially be one cluster but based on the total salary it should be three.

I need to do spatial clustering with a constraint on that total salary so it falls between two values.

Right now all I have is some Python for hierarchical spatial clustering but I don't know how to take into account the attribute!


texts: Disapearing features after specific zoom degree in ArcGIS 10.2

texts:
I have faced a Problem during my Analysis. I have some Polygon Features and after up to specific degree they would be disapeared.

I have checked all the Setting and no zoom degree in defined.

What could be the reason?


texts: I am new to QGIS.  I wanted to negate (behaving like NOT gate in a digital circuit) the MASK I have obtained earlier. That means I want to convert all 0s to  1 and all 1s to 0s. Is there any way to do it using raster calculation?


texts: What are the spectral indices that use addition of Landsat 8 bands for finding vegetation?

texts: NDVI is found out using the division of subtraction and addition of bands. Is there a simpler index for finding out vegetation(or other physical characteristic)? For example, DVI=NearInfraRed-Red is a simple linear operation for finding vegetation. There might be other such simple operations on bands for finding different physical characteristics of a region.


texts: Raster projected in Google Earth Engine showing weird values within R

texts:
I classified a raster in GEE and I am trying to work with the same within R. I have 7 land cover classes and when I view the output in R, I get something that does not make sense. Please see the values:

> lc2018
class      : RasterLayer
dimensions : 19261, 18167, 349914587  (nrow, ncol, ncell)
resolution : 10, 10  (x, y)
extent     : 633830, 815500, 1090410, 1283020  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=43 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
source     : E:/Land Cover Data/Reprojected Image_WG.tif
names      : Reprojected_Image_WG
values     : -2147483648, 2147483647  (min, max)

How do I ensure that the values tab goes from 1 to 7? When I plot the raster, it looks fine to me.


texts:
Do you know a way to map all the river gauges / stations in UK?
I found some data from the Environment Agency but it lacks of the longitude and latitude.


texts: Is Heatmap Color Ramp uniform across different rasters for ArcGIS Desktop?

texts:
I have about 16 rasters spanning the same extent and containing the same variable.

I use the display symbology in ArcGIS to generate a stretched Color Ramp which is the heatmap i want.

As the Colorramp colors adjust to maximum and minimum values within a given raster , does it makes sense to compare two different heatmaps with same colorramp but slightly different statistical distribution ?

Is there a way to generate a heatmap that is independent of statistical distribution and purely a function a variable value ?


texts: What kind of unit of measurement for scale does SQL Server 2012 Spatial use?


texts: A LiDAR user is telling me that his LiDAR data acquired from my organization has only integer LiDAR in certain areas. He is requesting non-integer LiDAR for those areas. First, can anyone help explain integer vs non integer LiDAR to me?


texts: Auto Panning to Cursor Position using ArcGIS API for JavaScript

texts:
Does anyone know if it's possible using the ArcGIS JS API to make it so the map auto-pans based on the user's cursor position if they get within a close distance to the edge of the map viewing area?

This would allow the user to continue editing and not have to do some other action (such arrow keys, or a pan icon) to get the map to pan to an area they want to continue to edit in.  I've been looking for some examples of this and I can't seem to find any.


texts: QGIS overview map of atlas - setting a color

texts:
I have an overview map of atlas in atlas print layout. How can i set an automatic change of color for atlas map in overview map?
E.G. I have 600 atlas maps. In oveview map every atlas map is a polygon with a black line and and no filling. I want for every atlas map change color of line of that 1 different polygon (atlas extent) in overview map.


texts:
I know this is pretty much a duplicate of Can you hide attribute columns in QGIS? but I have about 100 columns I want to hide so right clicking each column - Hide is not ideal.

Isn't there a way to hide many columns, then unhide the ones I want to see, like this in ArcMap:


texts:
Background

I'm trying to make contour/isoline/isochrone maps to show what areas are within walking distance to schools for my pet project to help with the Safe Routes to School program. I followed the pgRouting workshop and underdark's posts (all linked to from my project page) to import OpenStreetMap into a PostGIS database, then used the driving_distance() function to create a point layer that had a cost attribute, then ran the contour plugin of Quantum GIS to create contour lines, which I then styled in QGIS to create this image (pardon the spurious lines, not sure where they came from):

What I'd like to do

I'd like to replace the QGIS contour plugin and styling feature with GRASS, as I'd like to create this image dynamically for a web service, and I don't think QGIS is meant for that (at least performance-wise, I think it may be possible with Python). So I'm trying to figure out how to do this same thing with GRASS, which appears to be used quite a bit in web services.

From my reading, it appears I'll want to use something like v.surf.rst, though that seems like it would generate more of a heat map with a gradient, rather than the type of image shown above with solid colors for each band (in this case each band is 1/4 mile).


texts: Calculating time difference between consecutive rows does not work within graphical modeler

texts:
I have a dataset with time-stamped GPS points. I would like to calculate the time difference in seconds between these rows.  I successfullymanaged to do so by running the piece of code below in the field calculator outside the graphical modeler. However, when I want to calculate the time difference with the field calculator in a graphical model, it does not work.. It only generates NULL values, except for row with id 0. I used the exact same code. Any ideas what the problem is?

CASE WHEN id=0 THEN 0
ELSE
second(age(attribute($currentfeature, 'datetime'),
       attribute(get_feature_by_id('[layername]',$id-1), 'datetime')))
END


texts:
I have a point layer with geographic boundaries of an area and another point layer with road network. I want to calculate the distance of each point of road network layer to the points of the geographic boundary and create a raster with the results.

Any ideas??

I use ArcMap 10.1


texts: Modis Reprojection Tool "Input header file, Error in module"

texts:
I am trying to run Modis Reprojection tool for conversion from HFD to Geotiff for Modis terra Level 2 data 3KM swath i.e "MOD04_3K". Every time I try to input a file I get the following error:
Error in module: ReadHdrFile"

Any suggestions what can be done in order to solve this problem.


texts: QGIS 3.0.3 python error "/bin/sh: gdal_polygonize.py: command not found"

texts:
I have recently downloaded QGIS on my Mac to preform an analysis for my dissertation.

I am trying to convert my raster data to polygon data using the gdal_polygonize.py tool provided by QGIS.

Whenever I try to use any of the gdal software I get the error message:
"/bin/sh: gdal_polygonize.py: command not found"

Can anyone help me solve this problem?
I really need to use this feature for my analysis.


texts: "Bring All to Front" is missing in QGIS 3

texts:
I am using QGIS 3.4 and 3.10 and a plugin's window ends up hiding behind the canvas. In QGIS 2.10 I had the option "Bring All to Front" to bring all windows to front.

Is there an alternative for it in QGIS 3?

This is on Mac OSX High Sierra 10.13.6.


texts: Can a GeoJSON GeometryCollection contain another Collection?

texts: From a reading of the spec it seems that a GeometryCollection is a Geometry and can contain geometry objects? But I want to be sure that this is really allowed before I file a bug report.


texts: How to connect Qfield to PostgreSQL by online?

texts:
I created a database on PostgreSQL, connect it to QGIS through PostGIS.
Next, I package for the QGIS project and send it to my mobile. But it gets an error, my phone cannot read the database, it calls the error "PostGIS???". My postpose try to collect data by Qfield and add it online to PostGIS. Is there any document talk about this topic, step by step? I am not an IT and I am a beginner in this field.


texts: Filling null areas in raster to match extent of second raster using ArcGIS Pro

texts:
I have two global rasters and would like to expand one so that it matches the extent of the second by filling in NoData values with values taken from or averaged from nearby pixels. I have been using this expression in Raster Calculator in ArcGIS Pro:
Con(IsNull("fill_raster.tif")&"extent_raster.tif"==1, Int(FocalStatistics("fill_raster.tif", NbrCircle(15, "CELL"), "MEAN")),  "fill_raster.tif")
This has been working to fill most of the small areas of NoData, but there are still larger or more isolated regions where there are no cells in the "fill_raster" with values, so using the circular neighborhood at a specified size does not reach them. In the image below, I am trying to fill the green areas with data from the grayscale raster, but some of these green areas had no cells with data within the circular neighborhood and remain NoData in the grayscale raster even after applying the focal statistics expression.
Is there a way to specify that the values should come from the nearest cells with data instead of any cells within the specified neighborhood size?


texts: Change/filter markers from custom check boxes

texts:
So I'm working on a map using leaflet and I want to use custom check boxes, not the default controls. I'm returning the correct json information when I check the boxes. I'm just not sure how to swap out the json reference and trigger marker refresh. Any thoughts would be welcome. The project can be seen here: https://codepen.io/sharperwebdev/pen/PQvMqR?editors=1010

Basically what I'm asking is how do I swap out the original geoJSON for a reduced version on an external (outwith the map) click event.


texts: Can I publish a QGIS processing plug-in written in GeoPandas?

texts:
Or do I need to rewrite my script using "pure Python" and PyQGIS? My script is basically a recreation of ArcMap's Tabulate Area tool (vector layers only for now). The goal is to make it official in the QGIS plugin repository
Code snippet
data = []
for a in range(gdf1.shape[0]):
    muni = gdf1.loc[gdf1['NAME_2'] == gdf1['NAME_2'][a]]
    intersection = gpd.overlay(muni, gdf2, how='intersection')
    area = intersection['geometry'].map(lambda p: p.area)
    type_in = intersection['DESCRIPT'].values.tolist()
    type_all = gdf2['DESCRIPT'].values.tolist()

    ...

df3 = pd.DataFrame(data, columns = cols) # creates df
df3.insert(0, 'zones', rows) # creates df's rows
df3.to_csv('munis_landCover_tabulated.csv')


texts: Getting No module name exifread error from QGIS Photo2Shape plugin?

texts:
I had Photo2Shape plugin installed on QGIS 2.4 (Windows 64 bit), following installation of 2.6.1 and an upgrade of the plugin itself Photo2Shape is now broken.

Error message on install is "No module named exifread".

I've done an uninstall of all versions of QGIS, deleted all User folders, and done a clean install of 2.6.1 thinking this would reset all my plugins but they're still there? And Photo2shape still not finding Exifread.


texts: Changed Projected Coordinate System back to original but, Shape_Area field is still off by about 55 SF

texts:
I recieved some GIS parcel data a little over a year ago. This data is the basis for the majority of my work. About a week ago I noticed that the Projected Coordinate System had somehow been converted from NAD_1983_StatePlane_Massachusetts_Mainland_FIPS_2011_Feet into WGS_1984_Web_Mercator_Auxiliary_Sphere. I am not sure when this conversion happened as I have not needed to measure any parcels for quite some time. I noticed the change when I went to compare some of our assessing records against the GIS data and they were VERY different from one another. I have projected the data back into NAD_1983_StatePlane_Massachusetts_Mainland_FIPS_2011_Feet and the Shape_Area is now very close to the original data (I kept a backup copy of the original data on a disc) but is still short from the original by about 55 SF for the parcel I checked. Not sure if I missed something but, if anybody can help me to get the Shape_Area of these parcels to match the original data I would be very grateful.

Note - a few parcels have been changed from the original data so I cannot simply reload it. Also the to project my data back to its original Projected Coordinate System i went into ArcToolbox -> Data Management tools -> Projections & Transformations -> Project. Lastly, I'm using ArcGIS 10.2.2 and the file I'm working with a File Geodatabase Feature Class.

I wasn't sure whether to start a new thread for this or, attach add it to one of the many threads about coordinate systems that are already on here. I felt as though this is different from most of the questions on coordinate systems that I have seen on here, so I made the decision to open a new thread.


texts: Reading National Elevation Dataset (ArcGrid/GridFloat/IMG) with Python only tools?

texts:
I've found high precision elevation (1/3 - 1/9 arcsecond resolution) data from the National Elevation Dataset provided by the USGS. It comes int IMG, Arcgrid, and GridFloat formats. I know ArcGrid corresponds to a paid software package, but I'm trying to stick to using freely available tools.

I have GPS data that I'm trying to correlate with ground level.

Are there any python libraries which will let me transform gps data into altitude using a datafile in IMG, ArcGrid, or GridFloat formats?


texts:
I'm fairly new to GIS, working on a project where I overlay counties on google maps, and overlay cities on top. I'm using fusion tables, and running into Z-index issues. To get around this I'm hoping to take the county tiger gis shape files and "adjust" the lines to only draw areas where there are no "places" from tiger gis. It seems like I'm looking for a reverse clipping, any hints on how to accomplish this with ogr2ogr?

A copy of the script I'm using to get to this point can be found at: https://github.com/TinyHouse/shape_to_sql/blob/master/shape_to_sql