texts: JSONs made from shapefiles are not plotting (using folium or Leaflet alone)

texts:
My goal is to plot two shapefiles as layers on an interactive map.
The files I'm using live at https://github.com/AnnaMag/maps.

I transformed shapefiles to .geojson using

ogr2ogr -f ...

e.g. ogr2ogr -f geoJSON Wards2011.geojson Wards2011.shp

Neither of the files is being plotted using folium or Leaflet alone.

Code example for 1 shapefile:

#center in SA:
lllon = 15
lllat = -35
urlon = 33
urlat = -22

state_geo = r'data/Wards2011.geojson'
map = folium.Map(location=[(lllat+urlat)/2, (lllon+urlon)/2], zoom_start=6,
tiles='cartodbpositron')
map.geo_json(geo_path=state_geo)

map.create_map(path='test.html')

Replacing the geo_path with:

state_geo = r'data/sal_population.json'

is not working either.

Is there anything wrong with the json files or the approach?

UPDATE

As suggested I reduced the complexity:

in Wards shapefile I left only 4 necessary features (incl. geometry, of course)
Simplified the shapefiles

ogr2ogr out.shp in.shp -simplify 0.000001
transformed to geojson

ogr2ogr -t_srs EPSG:4326 -f GeoJSON Wards2011.geojson out.shp
transformed to topojson

topojson --spherical  --properties -s 0.001 -q 1E6 -o Wards2011.topojson Wards2011.geojson

Neither of the files is working. What should be fixed to make it work?


texts: Trying to make calculation with gdal_calc.py?

texts:
I have a tif image that I'm trying to run gdal_calc.py.
In the source file the values are between 1,118-9,704 and I want to divide all the values by 10,000 to get results between 0-1.

This is the command I run:

gdal_calc.py --type=Float32 -A C:\z\input.tif --outfile=C:\z\output.tif --calc="A/10000"

The problem is that the output is not working, it seems that it does not know how to get decimal numbers.

I display the output in ArcMap software

Anyone have an idea?


texts:
I have data at the ZIP3 level and I need to join it to a state file. I also need the data within the zip3 file to be averaged, but for some reason when I use the join tools, i'm not getting the results that I need.

Is there a way to aggregate data from a zip3 level to a state level without swtiching over to raster data (which is what the spatial analyst tool requires).

Thanks!


texts: Dissolve polygons for each column in ArcGIS

texts: I am wondering if there is a way to repeat the dissolve process for each column in order to have just one shapefile as outcome. In the image an example of input data and desired outcome.


texts: Convert fields of points into separate polygons based on their outlines

texts:
I have several fields of points that look like so:

They are actually thousands of points that appear as polygons because they're so close to one another. They're research points 720 feet from one another in all directions.

So clearly on this map we can see 7 would-be polygons. How can I get qgis to create these polygons without me having to draw them by hand?


texts:
First of all, I am new to map making, and I am mostly relying on tutorials to learn. As part of a project, I am making an airway map using QGIS and I do it by uploading the points with lat/long thru a csv file, and then connect them thru Points2One plugin. After that, I used pgrouting to create topology and calculate shortest path with dijkstra. However, when input the start and end points for calculation, the pgr_dijkstra won't return anything unless I used the "source" and "target" points. Below is the route network:

If I use point 5 & 6, it will yield result since in the table, there are the source and target. but if I use 29 and 23, it won't give result.

Is pgr_djikstra supposed to work like that? Anything that I am doing wrong?


texts: Setting ArcMap layout size based on web unit (pixel) standard?

texts:
I need to export my map from ArcMap as SVG and display it on a web page in an exact size of 500x500 pixels, like this:

Can you please let me know how I can set up the size of the layout to be exactly 500px in height and width?

I already found out how to change the size of the layout from File -> Page And Print Set up but in Map Page size there only 4 options available: Points, Inches, Centimeters and Millimeters.

Is Points equal to px in web? If not how can I do this?


texts:
I dont even know if QGIS is the right program for what I need..

What I need:
I need a program with maps(Denmark to be specific), that I could feed with a database of users through a CSV file.. e.g. give it 500 names, adresses etc. and have them dotted on the map or shown in different graphical ways, like perhaps a heat map.

What I did:
tried mappoint from microsoft, quickly realised this wasn't the right program, while it was able to accept my CSV file, the way it displayed my users was of no use to me, its a glorified 'route planner' and not suited for handling multiple selections at the same time.

Where I am at now:
I downloaded and installed QGIS thinking this was it, now opening it, I have no idea what to do, and it doesn't look like it holds any standard maps. so maybe this program has nothing to do with the solution I need?
if not, could you direct me to one?

Note: I'd really like to avoid uploading these CSV file and not expose my users information to the internet, which is why I found solution like google maps not to be the way to go.


texts: correct name of the files in the layers panel

texts:
I wrote a script for clipping all the raster files in a directory using a polygon shapefile. The script works fine, but the clipped layers are then uploaded in QGIS with the generic name Clipped (mask).

How can I do to preserve the correct name of the original files.

import os
import processing
from qgis.core import *

raster_filepath = "R:\\S1\\Season_2014_2015\\1_splitted\\VV"

mask_filepath = "R:\\RRE-    GIS\\Vector\\Border_scenes\\Rothamsted\\Rothamsted_box.shp"
mask_layer = QgsVectorLayer(mask_filepath, 'mask', 'ogr')

outputDir = "R:\\S1\\Season_2014_2015\\1_splitted\\VV"

# create the new folders clip
dirCreateClip = outputDir + "\\" + "clip"

if not os.path.exists(dirCreateClip):
    os.makedirs(dirCreateClip)

for i in os.listdir(raster_filepath):

    layercount = os.path.join(raster_filepath, i)
    raster_layer = QgsRasterLayer(layercount, 'raster')

    fileOut = dirCreateClip + "\\" + i
    fileOutClip = fileOut[:-4] + "_clip.tif"

    processing.runandload("gdalogr:cliprasterbymasklayer", \
    raster_layer, \
    mask_layer, \
    None, \
    False, \
    True, \
    True, \
    5, \
    4, \
    75, \
    6, \
    1, \
    False, \
    0, \
    False, \
    "", \
    fileOutClip)


texts: Converting shapes in a personal geodatabase to GeoJSON

texts:
Sorry if the question is too broad or not exact, but I'm pretty new in this area and not really sure about the namings and so on.

I have got an MS Access Personal Geo Database which is created either by ArcGIS Desktop or ArcSDE (I'm not sure as there is no clear statement about it on the source). This DB contains a table with a column named shape which is supposed to be Polylines with M value according to the manual.
Now I want to convert these polylines to something more usable specially if possible GeoJSON.

Unfortunately I can't decode the binary format in this column, some of them start with these 8 bytes: 17 00 00 00 which is according to this PDF the code for a PolyLineM but some other records (most of them) start with these 8 bytes: 32 00 00 60 which seems to be no valid value for any shape.

Also trying to open the binary data using tools that support ESRI's shapefile format results in errors which I think is obvious because it is not a full shapefile.

Can anybody point me to a Java library or another open source solution that can convert these records at least from command line? or at least tell me what format should I look for if I want to convert these records?


texts: Only Draw Geometric Objects Above 8th Parallel North

texts:
I have many geometric objects (in both SRID:4326 - the_geom - and SRID:3857 - the_web_geom - in PostGIS extension tables of PostgreSQL), including country boundaries, for the country of Ghana and I want to make a hard cut-off of the 8th parallel north latitude - i.e. I only want to draw the_web_geom (SRID:3857) object if it falls above this line.

Is there a way in either Mapnik, by adjusting the extent property of the layer, or by amending the PostgreSQL query to only draw the geometric features that fall above this line of latitude? I believe one solution, for the point objects, would be ...WHERE ST_Y(the_geom) > 8.0, but not sure if this is the best approach or how this would work with polygon objects.


texts: Where is the edit widget for photos in QGIS3?

texts:
I’d like to use qgis2web to display a photo when you click on a polygon in the map using qgis2web, but I can’t seem to find the correct settings.

I have a table with a column named photos. Each row has a field with a html link to an image. I’ve read that in version 2.8 you would simply choose “photo” in the field edit widget, but there is not a photo option in QGIS 3. The replacement method appears to be the “Attachment” option in the Attributes Form in the Layer Properties.

I have tried every combination I can think of and still can’t get it to work. The closest I have come is an html link that when clicked on, would open another window to display the image.

Using QGIS version 3.0.1-Girona and the latest plugin update for qgis2web

Can anyone let me in on the secret to making this work?


texts: How to use a specific AREA (km2) to create a buffer-like feature around a polygon (Area to Buffer)

texts:
How to create a buffer (or buffer-like feature) based on a specific area (in km2)?

I have a specific area (e.g. 14.5 km2) and would like to lay this area as a buffer (or something similar to a buffer) around a pre-existing, irregularly-shaped polygon.

This buffer should be of equal width around the polygon, so that I can extract the distance between the outer margins of the polygon and the newly created buffer.

Is there a way to do this? I'm using ArcGIS Pro 2.4.0.


texts: Labeling mulitple interests stacking problem

texts:
I am using ArcView 10.1 and am using this code to no avail in the Labeling Properties  [Company_Ne] &VBNEWLINE& [Exp_Date_O]

I know there is a way but am not familiar with the scripting for it.  I am trying to label each parcel but there are multiple tracts that have different expiration and company net acres.

One tract may have only one Lessee,just one exp date and company net acre.

another tract have multiple (ranges 2-9+) shared interests, so they will have different exp dates and different company net acres.  Is there a way to include all those in there without them stacking up on top of each other, but layered underneath one another?

There are 2 columns I am using for the labeling.  I am trying to make them look like the top Left (centered left side).


texts: How can I stop my map algebra being limited to the size of my smallest raster?

texts:
I have five input rasters and I'm using map algebra to find areas where all of them have a value over 75, areas where all of them are over 25, and areas where one or more are under 25. I can do this with:

Con((Raster("RISK-CH.tif") > 74.9) & (Raster("RISK-MI.tif") > 74.9) & (Raster("RISK-RE.tif") > 74.9) & (Raster("RISK-SE.tif") > 74.9) & (Raster("RISK-SO.tif") > 74.9), 100, Con((Raster("RISK-CH.tif") > 24.9) & (Raster("RISK-MI.tif") > 24.9) & (Raster("RISK-RE.tif") > 24.9) & (Raster("RISK-SE.tif") > 24.9) & (Raster("RISK-SO.tif") > 24.9), 50, 0))

The problem is that all of my rasters are slightly different shapes and sizes, and some have small holes, but I dont want to completely discount areas with no data. The output is the size and shape of the smallest input. I won't allow them to make the highest category, but if all the available datasets are above 25 I want the cell to still make my middle category.

I've tried using

  ((Raster("RISK-CH.tif") > 74.9) |
  IsNull(Raster("RISK-CH")))

but that gives me the exact same output, I've tried

((Raster("RISK-CH.tif") > 74.9) | (Raster("RISK-CH.tif") IS nodata))

and

((Raster("RISK-CH.tif") > 74.9) | (Raster("RISK-CH.tif") = NODATA))

but both of these throw errors.

Is the issue to do with my boolean logic, or can a raster calculator output not be larger than the smallest inputs?


texts:
I  have a Java Application. I am using QGIS to solve one such module of my Java Application. I have written code in python console within QGIS.Is there any way to integrate both of these?
I am thinking to create a plugin and then call this plugin from My Java application . The output File is .csv (and not shape Files). Is it possible to do ?


texts:
Is there a plugin for supporting opening KML files with colors, as currently the KML format is supported, but the colors and styles that comes with it are dropped?

I've created a KML style file:
http://www.files.com/shared/599d872e33689/sierraLeone.kml.zip

Uploading it to KML viewer, it's displaying properly, but in QGIS it fails. The kml viewer site:
http://ivanrublev.me/kml/


texts: Converting from UTM to WGS84 using GeoTools?

texts:
I have created a function that reads a shapefile file and I get its geographic coordinates well. This works fine. My problem is that I need those coordinates in wgs84 format (latitude, longitude), but I cannot get it to work.

I used a 'Query' object, but its setCoordinateSystem method does not do anything.

This is my code:

public Geometry getGeometryDistrict(String refCatDistrict, String shapeFileDistrict){

                Geometry union = null;

                File shapeFile = new File (shapeFileDistrict);

                try {

                   Query query = new Query();
                   query.setCoordinateSystem(DefaultGeographicCRS.WGS84);

                   Filter filter = CQL.toFilter("obj_co_id = '"+refCatDistrict+"'");
                   query.setFilter(filter);

                   FileDataStore storeDistrict = FileDataStoreFinder.getDataStore(shapeFile);

                   SimpleFeatureSource featureSource = storeDistrict.getFeatureSource();

                   SimpleFeatureCollection collection = featureSource.getFeatures(query);
                   SimpleFeature fDistrict = null;

                   ArrayList<Geometry> geometries = new ArrayList<>();

                   SimpleFeatureIterator itrDistrict = collection.features();

                   while(itrDistrict.hasNext()){

                      fDistrict = itrDistrict.next();

                      Geometry geomDistrict = (Geometry)fDistrict.getDefaultGeometry();
                      geometries.add(geomDistrict);

                   }

                   GeometryFactory factory = new GeometryFactory();

                   GeometryCollection geometryCollection = (GeometryCollection) factory.buildGeometry( geometries );

                   union = geometryCollection.union();

                   itrDistrict.close();
                   storeDistrict.dispose();

               }catch(IOException | NoSuchElementException | CQLException e){
                   Throwable cause = e.getCause();
                   System.out.println(cause);
               }

                return union;
            }

What should I change?


texts: Is there a Free Aggregate Points tool out there for ArcGIS 10?

texts:
I need to do some Getis-Ord analysis but first I need to aggregate my points. I know of XToolsPro, but the aggregate points function is locked and you must purchase a license to use it.   I am wondering if there is a free tool out there, or a workaround in ArcMap that I am missing to aggregate points.

Some context:  What I have are over 8,000 points, many of which have different attributes, but may share the same (approximate) geospatial location.  For this analysis, I just need to aggregate any points overlapping one another based upon a field containing a numeric value to be summed.  For example:

Point A overlaps Point B and Point C.
Point A field value = 1.4
Point B = 2.4, and
Point C = 5.2.

The Aggregate of points A,B, and C would return point D with a value of 9.0

I hope this makes sense.  Maybe I could get away with using something really basic but just haven't pieced it together yet.  Any help you can offer would be appreciated!


texts:
The server is very slow to produce an image when the client requests an image at a scale of 900 layers but only one layer has the requested bbox. I think that GeoServer loads information about all the layers that fall under the scale rule in the style, and then starts filtering them by BBOX. I think that if you set a bounding box in the style rule, and I think that then the GeoServer will not load information about all layers that do not fall into the BBOX of the request from the client.

  I tried to use cql filter, but I can't change the request from the
  client and add &CQL_FILTER=test to it&

http://192.168.0.163:8080/geoserver/gwc/service/wms?STYLES=&VERSION=1.3.0&SRS=EPSG%3A3857&TRANSPARENT=TRUE&LAYERS=cite%3AALLGroup&EXCEPTIONS=INIMAGE&SERVICE=WMS&REQUEST=GetMap&FORMAT=image%2Fpng&CRS=EPSG%3A900913&BBOX=6770791.9617694,7697820.240683,6771097.7098825,7698125.9887961&WIDTH=256&HEIGHT=256

I'm trying to set a bounding box in style. Can I set a rule in the layer style that specifies the BBOX I need?

  Each layer is a pyramid of images created using gdal from tiff. 100
  000 * 100 000 pixels.

<?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"
 xmlns:gml="http://www.opengis.net/gml">

  <!-- a Named Layer is the basic building block of an SLD document -->
  <NamedLayer>
    <Name>default_raster</Name>
    <UserStyle>
    <!-- Styles can have names, titles and abstracts -->
      <Title>Default Raster</Title>
      <Abstract>A sample style that draws a raster, good for displaying imagery</Abstract>
      <!-- FeatureTypeStyles describe how to render different features -->
      <!-- A FeatureTypeStyle for rendering rasters -->
      <FeatureTypeStyle>
        <Rule>
            <MinScaleDenominator>00001</MinScaleDenominator>
            <MaxScaleDenominator>10000</MaxScaleDenominator>
            <RasterSymbolizer>
                <Opacity>1.0</Opacity>
            </RasterSymbolizer>
        </Rule>

        <gml:BBOX>
        <PropertyName>GEOMETRY</PropertyName>
        <gml:Box srsName="urn:x-ogc:def:crs:EPSG:4326">
          <gml:coord>
             <gml:X>-10</gml:X> <gml:Y>0</gml:Y>
          </gml:coord>
          <gml:coord>
             <gml:X>10</gml:X> <gml:Y>10</gml:Y>
          </gml:coord>
        </gml:Box>
      </gml:BBOX>

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

getting an error:

  line 28: cvc-complex-type.2.4.a: Invalid content was found starting
  with element '{"http://www.opengis.net/gml/3.2":BBOX}'. One of
  '{"http://www.opengis.net/sld":Rule,
  "http://www.opengis.net/sld":VendorOption}' is expected.


texts: Capturing data in clip from smaller extent

texts:
What I am trying to perform is a clip what captures the shoreline on a town in CT. I have the town boundary and I can use the buffer to extrude the polygon but I have wondering if there was an easier way.

as you can see its the negative space in the boundary of the town polygon that I am hoping the fill.


texts: Merging OpenLayers3 overlay with a get feature on click command that closes

texts:
I have this JSFiddle that is working fine on its own to select cities in the UK.

I have then taken this OpenLayers3 example and put it into a JSFiddle.

I am looking to combine the two so I have a close button and some other text when I click on each city.


texts: Combine Shapefile with address information

texts:
I downloaded a UTM32 shapefile (.cpg, .dbf, .prj, .shp, .shx file, where other available georeferences are GK3 and TM32) that shows the exclusive division of Germany in its municipalities. Above, I have a large dataset of subjects with their postal addresses (street, house number, zip code, city).

I would like to assign each subject to exactly one municipality using the postal address. I want the corresponding municipality's id (saved as attribute in the shapefile) to be attached to each subject in my database and export it, say, in a csv file.

However, I am not sure how to combine the shapefile with the addresses. I can imagine two approaches:

Using Nominatim or other services I could geocode the postal addresses, i.e. convert them into lat/long coordinates. Let's assume for a second that I already did.

How can I add coordinate information to the shapefiles in order to achieve the assignment of subjects?

There are also shapefiles for German Streets (Open Street Map etc.). Let's assume I have access to such data.

Is it in principle possible to place the municipality and the street layer on top of each other, locate each address on the street layer and identify the corresponding municipality on this basis?

Which one is practical?

I am using QGIS 2.6 and Mac OS X Yosemite.


texts: Creating empty geodatabase with all feature datasets according to another geodatabase structure using ArcPy?

texts:
I have a file geodatabase which includes a lot of feature datasets and each feature dataset has lots of Features. I am going to do clipping through all the Features using ArcPy Code. However, the clipped Features are created in the same feature datsets. I would like to know how I can chnage this script to have a new FGDB with the same feature datasets but just include the clipped Features.

import arcpy, os
database = r'C:\Users\...\WEB_clip_4half.gdb'
clip_feature = r'C:\...\Kapsel_.shp

try:

    arcpy.env.workspace = database
    datasets = arcpy.ListDatasets(feature_type='Feature')

    for dataset in datasets:
        arcpy.env.workspace = os.path.join(database,dataset)
        features = arcpy.ListFeatureClasses()
        for feature in features:
            if '_clipped' in feature:
                arcpy.Delete_management(in_data=feature)
        else:
             arcpy.Clip_analysis(in_features=feature,clip_features=clip_feature, out_feature_class=feature+'_clipped')
except Exception as e:
  print("Found an error - {0}".format(e))


texts: How do I access versioned spatial data in SDE with SQL queries?

texts:
We're running SDE and have several feature sets with versioning turned on. When I query those layers via ArcGIS Engine or ArcGIS Server, all expected features are returned.

We have numerous non-GIS applications that use SQL queries to access information from the spatial data. However, we don't get all the features that I'm expecting to be returned. How do I access the features that have been added/deleted (and are, therefore, "hidden" from a straight SQL query)?


texts:
could someone help me with this method of arcobjects? I don't know what the parameter "hint" means. Thanks.

ISurface.Locate Method
Returns the intersection of the query ray and the displayed feature.
[Visual Basic 6.0]
Function Locate(
    ByVal pRay As IRay, _
    ByVal hint As Long _
) As IPoint

Parameters
pRay [in]
  pRay is a parameter of type IRay
hint [in]
  hint is a parameter of type long
ppLocation [out, retval]
  ppLocation is a parameter of type IPoint


texts: Can I manually set the size of GeoServer's tiles?

texts:
I'm using GeoServer to generate a street name overlay that we can use with our various basemaps.  It works great, but the names repeat too often to be truly beautiful.  I've tried bumping up the repeat option in the SLD, but high values are being ignored, and the label is generated once per tile.

The simplest workaround I can think of is to simply bump up the size of the tiles for the layer, but I can't find an option that'll let me set that.  Does this option exist in GeoServer?

Thanks in advance!


texts: Geotools Point in Polygon incorrect result vs QGIS

texts:
I need to find the point in a polygon to get the corresponding feature and the result does not match the display in QGIS. QGIS shows the information correctly but with geotools I can't get it right.

Code:

public Map<String, String> getLocalityRecognizer(double lat, double lng) throws MalformedURLException, IOException, NoSuchAuthorityCodeException, FactoryException {

    Map<String, String> nombreCordon = new HashMap<String, String>();

    FeatureIterator<SimpleFeature> shapeSource = getShapeFeatures();

    while (shapeSource.hasNext()) {
        SimpleFeature feature = shapeSource.next();

        //This fragment return the result
        if (feature.getBounds().contains(lat, lng)) {
            nombreCordon.put("localidad", Normalizer.normalize(feature.getAttribute("nombre").toString().toLowerCase(), Form.NFC));
            nombreCordon.put("cordon", feature.getAttribute("cordon").toString());

            return nombreCordon;
        }
    }

    shapeSource.close();

    return null;
}

private FeatureIterator<SimpleFeature> getShapeFeatures() throws MalformedURLException, IOException, NoSuchAuthorityCodeException, FactoryException {
    File shpFile = new File("C://Users//Sdev//Documents//shp//voys_meli.shp");

    // Read
    DataStore inputDataStore = DataStoreFinder
            .getDataStore(Collections.singletonMap("url", shpFile.toURI().toURL()));

    FeatureSource<SimpleFeatureType, SimpleFeature> source = inputDataStore.getFeatureSource("voys_meli");

    FeatureCollection<SimpleFeatureType, SimpleFeature> inputFeatureCollection =
            new ReprojectingFeatureCollection(source.getFeatures(), CRS.decode("EPSG:4326"));

    return inputFeatureCollection.features();
}

With the same point, the code returns "MERLO", but the point intersect "ITUZAINGO".

CRS is EPSG:4326

SOLUTION

    while (shapeSource.hasNext()) {
        SimpleFeature feature = shapeSource.next();

        MultiPolygon geo = (MultiPolygon) feature.getDefaultGeometry();

        GeometryFactory fac = new GeometryFactory();

        Point p = fac.createPoint(new Coordinate(lat, lng));

        if (geo.contains(p)) {
            ...
        }
    }


texts:
I have a modelbuilder model setup that, in short, takes all selected features from ArcMap and exports it to a dwg.  I need to add two things to this model which I'm not sure can be done in modelbuilder.  The ArcMap document has several feature classes which don't necessarily overlap so first I need the model to run only on the feature classes that HAVE selected features and ignore the rest of the feature classes.  Second, I need the model to NOT run if more than 10000 features are selected from any of the feature classes in ArcMap.  This can probably be easily done in Python but I have no experience in Python hence I'd like to confirm if it can be done in modelbuilder.  If not, I'll have to explore Python.

Edit: I need a statement where if the number of features are greater than 0 and less than 10000 then it should be true and proceed.  I think the "Get Count" and "Calculate Value" tools can be used to do this in modelbuilder, just not sure how.


texts: Clipping faster (i.e. in bulk/batch) in ArcGIS for Desktop?

texts:
I have 100 different input vector layers and one clip feature vector layer.

I have to clip them 100 times but its a large amount of work.

Is there any Python script (ArcPy) in ArcGIS for Desktop to accelerate my work and avoid the one by one clipping?


texts: Using MultiValue in Python script tool parameters?

texts:
My script works fine with a single feature class selected for clip_fc = arcpy.GetParameter(0) .

However when I set clip_fc to MultiValue in the script properties and then select several feature classes, I get an error "Parameters are not valid.".

Can anyone show me where I am going wrong?

import arcpy
import os.path
from arcpy import env

arcpy.env.workspace = r"C:\GIS Home\project_1\clipshps\Combined_mf.gdb"
arcpy.env.overwriteOutput = True
clip_fc = arcpy.GetParameter(0)
clip_name = str(clip_fc)
output_dir = r"C:\GIS Home\project_1\clipshps\outputs_clip"

for fds in arcpy.ListDatasets():
    for fc in arcpy.ListFeatureClasses('','',fds):
        out = os.path.join(output_dir, fc + "_" + clip_name[40:] + "_clip.shp")
        arcpy.Clip_analysis(fc, clip_fc, out, "")


texts: Add GIS functionalities to existing java based tool/application

texts:
I have a tool/application and want to add GIS functionalities to it. I will be mainly working with the vector data.
The required functionalities are:

Adding/importing the GPS point data and generation of polygons from the points.
Classification in maps based on various attributes.
Labelling of the features.
Generation of maps.
Spatial queries (attribute and location based).
Network analysis (finding shortest route between two points).
Add google/bing map as a base layer.

I came across that QGIS is able to all the functionalities mentioned above. Can anyone suggest me if it is possible to integrate the QGIS with java based tool/application or are there any java based GIS tools (user-friendly) able to perform these requirements.


texts: Filter out duplicate Sentinel-2 Images form Earth Engine Image Collection by date

texts:
I'm stacking all the images from Sentinel-2A from the L1C image collection, but when viewing the images beforehand, I see that there are duplicates images, same imaging date but different processing dates.

Sample:

features: List (182 elements)
0: Image COPERNICUS/S2/20150820T082006_20150820T082942_T36RYV (1 band)
1: Image COPERNICUS/S2/20150820T082006_20160929T004310_T36RYV (1 band)
2: Image COPERNICUS/S2/20150830T082006_20150830T082754_T36RYV (1 band)
3: Image COPERNICUS/S2/20150830T082006_20161006T204405_T36RYV (1 band)
4: Image COPERNICUS/S2/20150830T082754_20160405T115445_T36RYV (1 band)
5: Image COPERNICUS/S2/20150909T081736_20150909T082324_T36RYV (1 band)
6: Image COPERNICUS/S2/20150909T081736_20161014T144949_T36RYV (1 band)
7: Image COPERNICUS/S2/20150919T081736_20150919T082756_T36RYV (1 band)
8: Image COPERNICUS/S2/20150919T081736_20161022T211502_T36RYV (1 band)
9: Image COPERNICUS/S2/20150929T081736_20150929T082542_T36RYV (1 band)

As you can see, most of the images appear 2 times, some even 3 times.
So I want to filter out these images, and retain unique dates with the most recent processing date.

I'm not sure id there are any differences in the values, but to be on the safe side I don't want to reduce images together, but filter them out.


texts: How to retrieve time-series aggregated data from rasters based on a polygon in PostGIS?

texts:
I created a database in PostGIS with some MODIS data stored in the schema.table modis10v04.gpp. A simple query to show the structure of my data is given as follows:

SELECT rid, filename, acquisition_range FROM modis10v04.gpp LIMIT 5;

My output is something like this:

As you can see, there are repeated info in the first 4 rows, this is because of my rasters are divided in tiles. To clarify, my dates range from 2003-01-01 to 2016-12-31.

Then, I have my geometry data in the schema.table public.shape10v04 containing polygons with the column "id".

What I want is the average of my data for each polygon (within my shape10v04) for each date (from initial date to final date, e.g. 2005-10-10 to 2007-10-10).

I came across the post referred here:
Calculate mean value of polygon from raster in PostGIS, but it stills blurry in terms of implementing it for a time-series of rasters.


texts: Is there anyway to stack or merge .hdf files in R Studio?

texts: I would like to be able to merge, mosaic, or stack several .hdf files. Does anyone know the package or library i should be looking for the tool to do this?


texts: Effective way for creating Open Source Web mapping applications

texts:
Sorry for asking such a basic question, but these days there is lot of advancement on open source technologies.

My requirements for my application are below

Basic Map application with zoom, pan, etc
Dynamic Attribute Query based on the Layer selected on TOC
Input Data is from mostly Shape files / PostGIS

I know about GeoServer + GeoEXT or MapServer + Fusion or any other alternative if so please suggest

But need advice from real programmers for extending the application in further and also wanted to know best open source web mapping practice.

Thanks in advance


texts: Adding mouse move to my map using OpenLayers

texts:
I'm a new JavaScript coder.
I'm adding a mouse move to my map using OpenLayers but it's not working to display coordinate in my mouse.

1 - The first step I wrote this code in my js file:

map.on('pointermove', function(event) {
    var coord3857 = event.coordinate;
    var coord4326 = ol.proj.transform(coord3857, 'EPSG:3857', 'EPSG:4326');
    $('#mouse3857').text(ol.coordinate.toStringXY(coord3857, 2));
    $('#mouse4326').text(ol.coordinate.toStringXY(coord4326, 5));
});

2- Then I added this in the HTML file :

Coordonnées géographiques de la souris dans les différents systèmes de projection:

    EPSG:3857 > 0 / 0
    EPSG:4326 > 0 / 0

It's not working, and I don't know why.

I have two messages :

  ReferenceError: $ is not defined ($the dolar is replacing document.getElementById('')

If I change the $ to document.getElementById('mouse3857'), I have a message of:

  TypeError: document.getElementById(...) is null.

So do I need to do something else in HTML file or calling jQuery file?


texts: Similar to Free sources of mapping data for France, does any one know a source for free downloadable flood and precipitation maps for Belgium?


texts: GeographyFromText and ST_SetSRID+ST_MakePoint different results?

texts:
Why am I getting different results?

SELECT ST_DISTANCE(
  ST_GeographyFromText('SRID=4326;POINT(13.401111 52.383343)'), --berlin
  ST_GeographyFromText('SRID=4326;POINT(9.997559 53.533778)') --hamburg
);

= 262077.401614256

SELECT ST_DISTANCE(
  ST_SetSRID(ST_MakePoint(13.401111, 52.383343), 4326),
  ST_SetSRID(ST_MakePoint(9.997559, 53.533778), 4326)
);

= 3.592724162238034

The first result looks a lot closer to reality ;-)

Edit: had to cast the points in the second query to geography.

SELECT ST_DISTANCE(
    ST_SetSRID(ST_MakePoint(13.401111, 52.383343), 4326)::geography,
    ST_SetSRID(ST_MakePoint(9.997559, 53.533778), 4326)::geography
);

Still don't really understand why I have to do this. Isn't SRID=4326 enough for ST_DISTANCE to know what type of point I am passing?


texts:
I have a time that indicates a ship was detected at a certain point. I want to go through a list of the tracks of ships in the area(id'ed by their mmsi) and find the closest time right before the detection time and after the detection for each ship. This will help me ID which ship in the area was responsible for that detection at that time. The problem I have though is that I don't know how to find the closest time points.

import arcpy, os, datetime
arcpy.env.overwriteOutput = True

Source_File = "L:\\gathr\\indonesia\\Sara\\SAISData\\AAAS\\SAIS2.gdb\\AIS_Jan2016_CorrectTIme"
field = "mmsi"
mmsi_List = [row[0] for row in arcpy.da.SearchCursor(Source_File,"mmsi")]
mmsi_List1 = [x.encode('UTF8') for x in mmmsi_List]
UniqueMmsi = set(mmsi_List1)
UniqueMmsi_List = list(UniqueMmsi)
print UniqueMmsi_List
for mmsi_query in UniqueMmsi_List:
    out_layer = str(mmsi_query)
    in_features = "L:\\gathr\\indonesia\\Sara\\SAISData\\AAAS\\SAIS2.gdb\\AIS_Jan2016_CorrectTIme"
    #need 2 where clauses, one for the closest time before the detection time
    #and the other for the closest time after the detection time
    where_clause1 =???????????
    where_str1 = str(where_clause1)
    arcpy.MakeFeatureLayer_management(in_features, out_layer)
    arcpy.SelectLayerByAttribute_management(out_layer, "NEW_SELECTION", where1_str)


texts: Projection discrepancy with GDAL (1.9) and proj4

texts:
I boiled the problem down to that Python example:

import sys
try:
from osgeo import gdal
from osgeo import osr
except:
import gdal
import osr

if __name__ == '__main__':
version_num = int(gdal.VersionInfo('VERSION_NUM'))
print( "GDAL Version: " + str(version_num))

wkt = """
PROJCS["Popular Visualisation CRS / Mercator",
    GEOGCS["Popular Visualisation CRS",
        DATUM["Popular_Visualisation_Datum",
            SPHEROID["Popular Visualisation Sphere",6378137,0,
                AUTHORITY["EPSG","7059"]],
                TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6055"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4055"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
    AUTHORITY["EPSG","9001"]],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
    AUTHORITY["EPSG","3785"]]"""

X = 2662981.191207001451403
Y = -1935374.062384711578488

srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)

srsLatLong = srs.CloneGeogCS()
ct = osr.CoordinateTransformation(srs, srsLatLong)
(int, lat, height) = ct.TransformPoint(X, Y)

# Report results
print('X: %s\t\tY: %s' % (X, Y))
print('Lat: %s\t\tLon: %s' % (gdal.DecToDMS(lat, 'Lat', 2), gdal.DecToDMS(int, 'Long', 2)))

my output is:

GDAL Version: 1900
X: 2662981.19121                Y: -1935374.06238
Lat:  17d 1' 1.01"S             Lon:  23d55'19.08"E

However, when I run this from proj4, I get a different answer:

invproj +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs
2662981.191207001451403 -1935374.062384711578488
23d55'19.081"E  17d7'29.797"S

The latitude is off by 6'28" or so... ??????
What could be going wrong?


texts: Differences in costs of open-source hosted GIS software?

texts:
I recently saw a listing of the available GIS products at thegismarketplace.com, which led to a couple of questions.

What is the difference between a paid version of GeoServer or MapServer and the free ones?
I also noticed that the site is charging for QGIS, which I don't understand either since QGIS is free.

There is a related discussion here: comparison of open sources gis servers


texts: Unicode Decode Error while geocoding MMqgis

texts:
Running Qgis, trying to geocode a CSV file using MMQGIS

Qgis gets to the 66th element (there are 119) and returns the error:

Traceback (most recent call last):
  File "C:/Users/BPolinsky.MDC/.qgis2/python/plugins\mmqgis\mmqgis_dialogs.py", line 1198, in run
    shapefilename, notfoundfile, fields, service, 1)
  File "C:/Users/BPolinsky.MDC/.qgis2/python/plugins\mmqgis\mmqgis_library.py", line 1901, in mmqgis_geocode_web_service
    attributes.append(unicode(row[z], 'utf-8').strip())
UnicodeDecodeError: 'utf8' codec can't decode byte 0xc9 in position 10: invalid continuation byte

Python version:
2.7.4 (default, Apr  6 2013, 19:54:46) [MSC v.1500 32 bit (Intel)]

the question is, is the error something that I can change by either using a different plugin, playing around the the content / map settings, or would I need to go into the python of Qgis? If so, how would I go into its python? Why am I getting this error?

The other questions I've seen had to do with python mostly, not python filtered through qgis.


texts: Configuring Travel Modes in ArcGIS Network Analyst?

texts:
I am a new user of ArcGis and I am trying to do the exercises provided on the ESRI website for Network Analyst (https://desktop.arcgis.com/en/desktop/latest/guide-books/extensions/network-analyst/exercise-2-creating-a-multimodal-network-dataset.htm).

When doing the second exercise, after configuring the PedestrianTime evaluator I should configure the Travel Modes, but I don't get any window with this, and instead it jumps to the part from where I can configure directions.

How do I get the window needed for Travel Modes?


texts:
I am a new user of GRASS GIS.
How can I mosaic different satellite imageries using GRASS GUI?
I was trying to use i.image.mosaic, but could not succeed. I have been looking for an online tutorial but could not find one.


texts: Loading / georeferencing OS rasters in QGIS

texts: I have got some OS rasters but there isn't georeferencing file (TFW) with them. When I try to load them, it will obviously load in the wrong place. They are loading fine in a different GIS software (CadCorp) because the tile name is recognised (example:SJ86SW). Is that possible in QGIS? If not, while georeferencing the file can I snap to the corners of the raster to get the absolute position of these points?


texts:
I want to plot a diagram like this:
http://matplotlib.org/examples/shapes_and_collections/scatter_demo.html

Unfortunately, I have no idea what these diagrams are called and therefore I do not know where to search for an explanation of how they are created.

Does anyone know how they are called and whether they can be created in matplotlib?


texts:
I am trying to generate a map in Qgis which makes it possible to see the temperature of a region during a heatwave. This was a great succes. However, I want to add shades as well, since shades lower the temperature.

Now I have the following map:

.

On this map we see a part of a city, with red being warmer and blue beging colder. Tree shades for instance make some parts of a street cooler. As you can see some area's have become very pixel-ish. I would like to smooth this map a bit so that for instance one orange pixel within a group of blue pixels, takes over the value of its neighbour and turn into blue as well.

Does anybody has an idea for this?


texts: How to convert multipoint attribute to polyline in QGIS?

texts:
I have a multipoint layer in QGIS with 860 points.

Each points has some attributes value from a join, and I need to use one of those attributes (a column) to draw polylines. I need one polylines for each different value of that column (so 180 polylines)


texts: How to retrieve the processed file name, file path and save the result in a csv file using arcpy?

texts:
We have a list of feature-classes which are in multiple gdb's in multiple folders. This code is opening the csv and doing "GetCount_management" for the listed feature_classes alone and printing the output. The code is printing correctly. Say, if there are 25 feature classes in the csv, it's listing all the counts of all the 25 files. But while saving it in csv, only the 25th count i.e the last file's count alone is getting saved. I don't know whats the issue.

import csv
import arcpy
from os import path
csvpath = r"D:\RTXQAQC\Open2.csv"
output = r'D:\Test\Test_2.csv'
csvdir = path.split(csvpath)[0]
with open(csvpath, "r") as csvfile:
    fcs = [r[0] for r in csv.reader(csvfile)]
for fc in fcs:
    result = arcpy.GetCount_management(fc)
    count = int(result.getOutput(0))
    print(count)
with open(output, 'wb') as f:
    csvwriter = csv.writer(f,delimiter=',')
    csvwriter.writerow([fc,result])
    f.close()


texts: Error running SEXTANTE in QGIS on Scientific Linux

texts:
I've installed SEXTANTE into qgis on my Scientific Linux machine and I am getting the following when attempting to run Warp;

INFO|Mon Feb 24 2014 15:56:20|GDAL execution console output|/bin/sh: gdalwarp: command not found

Has anyone else come across this?


texts:
I would like the following to select (highlight) the "project" layer that intersects the grid, and then pull a record from the selected "Grid". I do not need to save the layer.

Note: Running the script from within ArcMap in the Python Console window.

import arcpy

# Make a layer from the feature class
arcpy.MakeFeatureLayer_management("project","project_lyr")

# Then add a selection to the layer based on location to features in another feature class
arcpy.SelectLayerByLocation_management("project_lyr","INTERSECT","Grid")

This code runs and no errors are given, yet it does not select(highlight) any features. Any Idea why? Also runs really slow.

Basically I'm trying to do this:

This runs really fast almost instant, but I want to do the same using python. I just need a value from the selected feature's table, I don't need a copy of the layer.

Part of the solution may be to use the in-memory workspace which is

  Writing data to the in-memory workspace is often significantly faster than writing to other formats...

Just not sure how to implement in-memory workspace


texts: How to group XY data based on certain categories in ArcGIS

texts:
I'm a newbie with ArcGIS.

I have a .csv file with the boundary latitude and longitude for various polygon( or islands) and each of these polygon fall within a few categories of region. I need to generate the following file types  (for polygon) for further use - .shp, .shx, .dbf, and .prj. But before i do so i need to group these polygons based on their respective regions and this is something i have been struggling with in ArcGIS 10.4.


texts: Selecting Untitled polygons in Google Earth

texts: I have an area where I'm trying to mark some polygons in it and under temporary places folder they were adding as shown below. Now how do I select some of the polygons and move them to different folder.


texts: Converting RFT legend template into a LYR-file for ArcGIS Desktop 10.x

texts:
I have access to an ArcGIS Server 10.5 Mosaic-service which has a colour-ramp symbology created from an RTF-template. I need to show the legend in ArcGIS but the mosaic is symbolised in the table of contents as a normal RGB raster.

How can I recreate the colour-ramp in ArcGIS from the RFT-file?


texts: How do I develop some of WFS and WMS methods?

texts:
I want to develop a GIS application. My application has some options as following:

Users with different access to layers
Layers with different source like postgis and  shp
OpenLayers-3.6.0 in client side
OpenLayers layers get its data via WFS and WMS protocols.

I want to use Spring java framework as application structure and Hibernate for Non Spatial Data transactions and geotools library for Spatial Data.

The best solution for me is using GEOSERVER, but I'm not allowed to use it.

I don't need all method of WFS and WMS, I need methods as following:

getFeatuers, transaction and getCapabilities methods in WFS Protocol.
getMap, WMSGetFeatureInfo and 'getCapabilities methods in WMS protocol.

I have some questions:

Can I separate WFS and WMS module of geoserver and add these to my project and use these? If yes, Any reference?
How can I develop the methods of WFS and WMS Protocols(I don't want to have a WFS and WMS server, I only want to implement some of its methods)?


texts:
I have written an application that does a lot of geoprocessing using arcpy. Once started, the application itself runs at an acceptable speed, but it is very slow to boot. It can take 30 seconds or so to start up the GUI, so I'd like to find ways to make it start faster.

As far as I can tell, the main source of delay is when I set the dataframe as follows:

self.df = arcpy.mapping.ListDataFrames(self.mxd, "")[0]

Does anyone have any tips / tricks to make either the whole application or the ListDataFrames run quicker when the application is started?

I use ArcGIS 10.1, python 2.7 and tkinter for the GUI. Feel free to ask if I've forgotten to mention anything.

EDIT:

As requested, here is the (very!) stripped down version of the start my code:

import arcpy
import sys
from Tkinter import *
import ttk
import tkMessageBox
import tkFileDialog
import webbrowser
import time
import csv
import os
import textwrap

arcpy.env.workspace = r".\..\Append\temp.gdb"

Title = "Automated Map Generator"

class Application(Frame):
    def __init__(self, master=None):
        Frame.__init__(self, master, background = "white")
        self.grid()

        self.createWidgets(master)

    def createWidgets(self, master):

        # setup map document
        self.mxd = arcpy.mapping.MapDocument(r".\autolim_mapping.mxd")
        arcpy.env.overwriteOutput = True
        self.df = arcpy.mapping.ListDataFrames(self.mxd, "")[0]

#[Lots of missing code here]

root = Tk()
root.tkraise()
app = Application(master=root)
app.master.title(Title)
app.mainloop()

How do I know what it's the ListDataFrames that takes so long? I put a whole heap of print statements through the startup code and noted how long they took to come out... I bet there are better ways to do this, but I'm pretty new to python.

EDIT 2:
@dassouki - thanks! Here are the output times from after 'import time':

Start Class Application (Frame):  1.53200006485
Start def Create Widgets:         1.59400010109
Start  mapping.MapDocument:       1.6099998951
Start env.overwrite:              1.64100003242
Start mapping.ListDataFrames:     1.64100003242
Finish mapping.ListDataFrames:    12.2660000324


texts: how to display the legend when the map contains external icons in R

texts: I haven't found anything about the legend when the map contains icons external R (in png), is't possible to use the legend function or another function in R.


texts:
I am trying to create a Python array using data from an ArcGIS attribute table. Here's the code I've written so far:

import arcpy
import numpy
input = "c:/data/from/arcgis"
arr = arcpy.da.TabletoNumPyArray(input,('OBJECTID','SoilEvap5'))
a = arr["OBJECTID"]
b = arr["SoilEvap5"]
np.columnstack((a,b))
print arr

When I try running this in the ArcGIS Python window, it returns:

Runtime error
Traceback (most recent call last):
  File "<string>", line 7, in <module>
AttributeError: 'module' object has no attribute 'TabletoNumPyArray'

Why is this error showing up and how can I fix it?


texts: Looping based on a spreasheet using python (in ArcGIS 9.3)

texts:
Here's what I am trying to do: I want to move feature classes from one .gdb (OLD) to another (NEW), but only the feature classes that are present in an excel spreadsheet.

read from excel spreasheet (old name column)
compare to old .gdb feature classes to find match
move matched feature class to new .gdb
rename feature class in new .gdb based on new name in excel spreadsheet (new name column)
continue doing the above until the end of the spread sheet


texts:
I recently updated my Mac to Sierra OS. In the past, I had use QGIS without any problem with my previous OS version (Lion). Since I installed Sierra, I can open shapes, visualise maps and export images, and also save my project. However, if I want to re-open the project, I see nothing on the display, even though I see the project shapes to the left. So I reinstalled QGIS from king chaos, I now have 2.16 version. Still the same. I can work with QGIS, save project, but when I re-open the project, the display is empty.

Does anybody know what is happening and how to fix the problem?


texts:
I have a Landsat raster time-series (2000-2016) of 102 rasters for the same location. Since the original Landsat scenes don't overlap exactly, I defined a polygon to use for cropping all the rasters to the same extent, since I need to stack them.

I work in R and used this function to batch process my time-series:

processLandsatBatch(x=listzip, pattern=glob2rx('*.zip'), outdir=output_gz, srdir=srdir, delete=TRUE, vi='ndmi', mask='fmask', keep=0, e= Riau_sp, overwrite=TRUE)

Riau_sp is my extent polygon.

However in the output I see that the rasters still don't have exactly the same extent, and I get ERROR when I try to stack them. Here is an example of two of my rasters extent after processing:

extent: 101.6222, 102.8827, -0.6339126, 0.6637299  (xmin, xmax, ymin, ymax)

extent: 101.6222, 102.8827, -0.6339154, 0.6637271  (xmin, xmax, ymin, ymax)

The y-values are just a bit off, even though both were cropped with the same extent-polygon.

Any idea what has happened, and how I may fix it?


texts: Why is Make XY Event Layer executing without error but not creating new file using ArcPy?

texts:
I'm trying to create a feature class from a table and I've gotten as far as you can see in my code below.  I noticed that after running the code through the arcpy.MakeXYEventLayer_management function, the code returned no errors but the January2018_points file never showed up in Output.gdb.  Because of this, the next section of code where I use Feature Class to Feature Class does not work, and it returns an error saying the input parameter does not exist or is not supported.

Is there something I'm missing here?

    # Import modules
    import arcpy
    from arcpy import env

    # Set workspace
    env.workspace = r"E:\GIS 4090\Projects\Project 1\Data"
    env.overwriteOutput = True

    # Import CSV file into a geodatabase using Copy Rows
    arcpy.CopyRows_management("January2018.csv",\
    r"E:\GIS 4090\Projects\Project 1\Data\Output.gdb\January2018")
    print("Table imported to gdb.")

    # Make XY Event later from the table
    arcpy.MakeXYEventLayer_management(r"Output.gdb\January2018", "XCoord",
    "YCoord",\
    r"Output\January2018_points")
    print("Temporary XY Event layer created.")

    # Convert XY Event layer to a feature class using Feature Class to
    Feature Class
    out_path = r"E:\GIS 4090\Projects\Project 1\Data\Output.gdb"

arcpy.FeatureClassToFeatureClass_conversion(r"Output.gdb\January2018_points", out_path, r"Output.gdb\January2018_points")
    print("XY Event Layer converted to feature class.")


texts: Atmospheric correction for Landsat 8 data in Google Earth Engine?

texts: What is the code for this? Already mosaic'ed my image


texts: Creating heat map animation based on point-shapefile with data attribute in ArcGIS

texts:
I want to create an animation of forest fires, using the NASA' MODIS data for the last 11 years. Their shapefile has acquisition data attribute — and it goes day by day all those years.

I've already tried to create an animation based on points only — looks good, but not so perfect and representative as heat maps could be.

So what I want to do is to create day by day animation (at least for one year, just for the start), but transforming these points to the heat maps and get the whole animation.

What advices can you give?

Maybe there's some drill with arcpy code that I can do?

Or there are other off-code algorithms to deal with?


texts: How do you emulate a right click on a tablet/mobile device using OpenLayers?

texts: I'm trying to make my OpenLayers mapping application work for mobile devices and one feature I haven't figured out how to emulate is the 'handleRightClicks'  Control Click option.  Anyone know how to make this work with tablets/mobile devices?


texts:
I am creating a directory which uses locations to categorize all listings. Instead of relying on user input, I'd like to provide a list of locations which the user can select. I have already tried GeoWorldMap which is provided by GeoBytes but it has certain issues. For example, it lists 275 countries. Other than that, their database is updated on a regular basis so that means if I maintain a local database, I'll need to update it every now and then to make sure it's up to date.

I want to know what developers use for this. Instead of a database, do they use some service to fetch a list of locations?

I am afraid manual user input is not an open option because categorization by location is going to be essential so I cannot trust or rely on what the user enters manually.

I have seen few sites doing this so I am sure there is a practical solution out there.


texts: While storing feature in ArcGIS 10 web application getting Class is not licensed for use

texts:
I have tried googling this issue

  Class is not licensed for use

but I could not found Any Solution .

I am using ArcGIS Server 10.0 of 32 bit and IBM DB2 32 bit,there is no problem while startEdting operation but  while storing feature(pFeature.store()) it is giving me this issue.

Has anyone already faced this issue?


texts:
I'm using QGIS 2.18.3 (Las Palmas).

I added the attribute Table in my print composer and would like to display it in 2 column (I have around 50 Cities I would like to display. Instead of having one long column that would need 2 pages, I would prefer splitting the attribute table in 2 columns).

It has worked before, but it was more by luck: for some reason my print composer automatically assumed I would need two pages, when I set the page-Number back to one, my attribute Table was split.

However, this is not working at the moment, and I would like to know if there is a better way of achieving what I need.

Is there a better way of working with attribute Tables in Print composers?


texts: Get values of terrain elevation of roads in QGIS?

texts:
I have one DEM raster layer with elevation values and one vector layer of the roads; same size, extension and CRS.

I want to do what this guy did here:

http://youarealegend.blogspot.de/2012/07/computing-roads-slopes-with-qgis.html

But instead of using slope I want to use simple elevation.

I converted my raster into polygons and tried to intersect (Vector > Geoprocessing > Intersection) it with the roads layer, but apparently is too much for Qgis to calculate, since it keeps crashing even before 2% of the processing. I would like to know another approach for doing this?

What I need in the end it's a table/data.frame that relates the "lines"(roads) and the elevation. Because I will use this in glm model.


texts: openlayers-3 mapserver WMS display issue

texts:
I am trying to include a mapserver WMS in Openlayers-3 using custom projection without proj4j.js, according to this. I have no image, no error message and Openlayers is getting mapserver .map according to debug tools, window network.

I suppose a problem of projection (?) but can't identify it.

<!doctype html>

<title>OpenLayers 3</title>


texts:
I am trying to figure out how to write metadata to ArcGIS Online Feature Service from ORACLE SDE and how to update the existing Feature Service metadata,  using FME or other tools.

I am able to read metadata from SDE using the "geodb_metadata_string" format parameter when setting the feature read mode to "Metadata" on the reader, but not being able to write to ArcGIS Online.

Is this even possible?


texts: Understanding QGIS scale dependent visibilty?

texts:
I cannot understand why the scale dependent visibility function does not want to work in QGIS.

I'm using the latest version (2.8.1 Wien) and want to insert diagrams (pie/histogram) to illustrate different statistical data.

It is a very simple and easy task to do, but for some reason the diagrams does not adhere to the minimum and maximum scale parameters I've entered.

Is there any assistance out there? Am I entering incorrect numbers? Is this scale based function different than rendering labels? What is it I am not understanding?


texts: Reducing the density of points based on zoom level

texts:
I have a Polyline like the one shown in the example, representing the path of a vehicle. GeoJSON markers show the vehicle's position over time.
How can I reduce the density of the markers when the user zooms out?

Example with Polyline Decorator


texts:
I tried to convert single LOD3 building using FME. I am able to convert successfully from DGN to cityGML format.But when I tried to view the same in FZKViewer/FME data inspector the performance is very slow. It takes more time to navigate.

Is there any way to control this?


texts: Control number of segments in ST_Buffer for geographies?

texts:
I've been attempting to construct circles centered around a list of points in PostGIS with good success. I've opted to circumvent having to look up UTM zones and using optimal SRID by simple letting ST_Buffer do the work for me by supplying a geography rather than a geometry.

However, in doing so I lose the ability to choose the number of segments for the resulting circle geometry, which defaults to 33 (not sure how this number is calculated, I make circle areas of radius 20 km). Is there any way to circumvent this and choose the number of segments myself? For my query of points I use:

update <table>
set poly = geometry(st_buffer(geography(point), 20000));

The output geometry will default to SRID 4326.


texts: Finding and deleting duplicate records in attribute table in QGIS

texts:
I have around 1 Lakh+ (100,000+) point features in my shapefile, which has some duplicate records in it. I had found the duplicate records present in the attribute table and print those in python console. But i want to select those duplicate features in attribute table so that i can delete those entries from attribute table. How should this can be done using python?
I had attached my snippet for printing the duplicate records present.

from collections import Counter
layer = iface.activeLayer()
idx = layer.fieldNameIndex('name')
attributeList = [feature.attributes()[idx] for feature in layer.getFeatures()]
print([attribute for attribute, count in Counter(attributeList).items() if count > 1])


texts: QGIS 2.14 Debian Sid - Couldn't load plugin 'processing'

texts:
The processing plugin can't load because python import error for QtWebKit.QWebView.

QGIS is installed from Debian official repositories. The plugin worked fine until 2.14.1.

The error message at startup :

Couldn't load plugin 'processing'

ImportError: No module named QtWebKit.QWebView
Traceback (most recent call last):
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 281, in loadPlugin
    __import__(packageName)
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 572, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/usr/share/qgis/python/plugins/processing/__init__.py", line 29, in
    from processing.tools.general import *
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 572, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/usr/share/qgis/python/plugins/processing/tools/general.py", line 28, in
    from processing.core.Processing import Processing
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 572, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/usr/share/qgis/python/plugins/processing/core/Processing.py", line 44, in
    from processing.gui.Postprocessing import handleAlgorithmResults
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 572, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/usr/share/qgis/python/plugins/processing/gui/Postprocessing.py", line 39, in
    from processing.gui.ResultsDialog import ResultsDialog
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 572, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/usr/share/qgis/python/plugins/processing/gui/ResultsDialog.py", line 38, in
    os.path.join(pluginPath, 'ui', 'DlgResults.ui'))
  File "/usr/lib/python2.7/dist-packages/PyQt4/uic/__init__.py", line 211, in loadUiType
    exec(code_string.getvalue(), ui_globals)
  File "", line 52, in
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 572, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
ImportError: No module named QtWebKit.QWebView

Python version: 2.7.12rc1 (default, Jun 13 2016, 09:20:59) [GCC 5.4.0 20160609]
QGIS version: 2.14.3-Essen Essen, exported

Python Path:
/usr/share/qgis/python
/home/pierma/.qgis2/python
/home/pierma/.qgis2/python/plugins
/usr/share/qgis/python/plugins
/usr/lib/python2.7
/usr/lib/python2.7/plat-x86_64-linux-gnu
/usr/lib/python2.7/lib-tk
/usr/lib/python2.7/lib-old
/usr/lib/python2.7/lib-dynload
/usr/local/lib/python2.7/dist-packages
/usr/lib/python2.7/dist-packages
/usr/lib/python2.7/dist-packages/PILcompat
/usr/lib/python2.7/dist-packages/gtk-2.0
/usr/lib/pymodules/python2.7
/home/pierma/.qgis2//python

Not sure if the issue comes from QGIS 2.14.2 update or from Qt/pyqt update. Any ideas to fix this?


texts:
Is it possible to only show one certain point of a layer in the print composer, but have them all shown in the QGis main window?

For an inventory of buildings, I want to have an overview of all buildings on the map in the main window, but in the layout, I want the atlas to only show one building at a time, so that there can't be any confusion about which building is being described.

I know how to insert a rule to only show one certain point ($id =  @atlas_featureid), which is great for the actual layout output. However, with this rule, I also only see one point in the main window and only have access to the other ones through the attribute tabel, instead of just being able to click on them (plus: I completely lose overview of the amount of points). At this moment, I just delete the rule as long as I'm working in the main window and not printing anything and then adding it, when I want to print my atlas. But this is of course a bit laborious, as we have to print the layout in different phases of the project and I have to add and delete that rule several times.


texts: GDAL C# Bindings for VRT don't write a VRT file

texts:
I have previously asked this question on stack overflow, and am now asking this more specialized audience.

I'm wondering why the GDAL C# library (v2.3.3) call to Gdal.wrapper_GDALBuildVRT_names(...) does not produce a physical destination VRT file.

When I issue the equivalent command in the OSGeo4W shell (to combine the 0000a013.gn1 and 0000b013.gn1 files), a VRT file is written to the supplied destination path (C:\Temp\z7_x43_y43.vrt):

gdalbuildvrt C:\Temp\z7_x43_y43.vrt F:\gnc\1\0000a013.gn1 F:\gnc\1\0000b013.gn1

However, when I issue the command in C#, no VRT file is written to the supplied destination path (I've made sure the folder exists and that the file isn't already there):

GdalConfiguration.ConfigureGdal();

var vrtFile = @"C:\Temp\z7_x43_y43.vrt";
var aeroFiles = new List<string>();
aeroFiles.Add(@"F:\gnc\1\0000a013.gn1");
aeroFiles.Add(@"F:\gnc\1\0000b013.gn1");

var vrtOptions = new GDALBuildVRTOptions(new [] { "-overwrite" });
var vrt = Gdal.wrapper_GDALBuildVRT_names(vrtFile, aeroFiles.ToArray(), vrtOptions, null, null);

The *.gn1 files are NITF, WGS84, degree unit files that work fine in other respects (e.g., the OSGeo4W command processes them perfectly).

The resulting vrt instance is correct - I can use it as a Dataset in other parts of my code, and can even cut PNG tiles from it doing something like this (pngFile and translateOptions are defined elsewhere):

Gdal.wrapper_GDALTranslate(pngFile, vrt, new GDALTranslateOptions(translateOptions.ToArray()), null, null);

But why is no VRT file ever written to disk? We have to supply a destination path (vrtFile), so doesn't that imply a destination file will be written?


texts: Why does "Points along geometry" in the processing toolbox take so long and is there something faster?

texts:
The points along geometry processing algorithm

  Creates points at regular intervals along line or polygon geometries. Created points will have new attributes added for the distance along the geometry and the angle of the line at the point.

I have a LineString with about 230,000 vertices and I picked an interval which would give me about 50,000 points in the end. The process takes 1+ hours (still going.) All I need are the equidistant points. I suspect that the extra "angle of the line at the point" is taking a lot of the time. The documentation does not give any visibility into how the algorithm works. Does anyone know if this is correct?

And is there a plugin which is faster for getting equidistant points? In the end I will convert them back into a linestring using the points to path tool. So if these is something that returns a reduced LineString instead of Points layer I could use that as well.

*QGIS verion 3.10.5, Windows


texts:
I have been trying to display a satellite image from Google Earth Engine in QGIS.

Then using both the QGIS and Google Earth Engine I want to process the satellite image.

Is it possible?


texts: GEE: resampling MODIS LST from 1km to 200km produces gradient

texts:
I am trying to resample raster from 1km resolution to 200m using this code:

var august = ee.ImageCollection(landsat)
  .filterDate('2019-08-01', '2019-08-30')
  .filterBounds(geometry);
  //Mean temperatre for August
var mean = august.reduce(ee.Reducer.mean());

// Get the projection information for a band. test 2
var Day_LST_mean = mean.select('LST_Day_1km_mean');
print('CRS:', Day_LST_mean.projection().crs());
//Resample to 200m
var Day_LST_200m = mean.resample('bilinear').reproject({
    crs: Day_LST_mean.projection().crs(),
    scale: 200
});
//testing resampling
//var mean_200m = mean.resample('bilinear').reproject(crs=crs, scale=200)

Map.addLayer(mean.clip(donbas), imageVisParam, "august mean 1km");
//print(mean);
Map.addLayer(Day_LST_200m.clip(donbas), imageVisParam, "august mean 200m");

What I get is a gradient that does not resemble the original.
Can someone clarify where is error in the logic?
Sorry if the question is dumb, I am just beginning my GEE journey.


texts:
I have a street map where one street(consisting of 2 lanes each) is shown as one line. If I click on that line, there are two objects (one object for each lane) listed like in the following image:

For each object, I do have the From/To Node which gives me the direction of the lane. What I want to do is to assign street signs (which are represented as blue points) to the nearest lane. That is why I need to move both Lane objects in the opposite direction so they approach the proper street sign. In the end, it should look like this.

So, with respect to the direction of the vector, the vector must always be shifted to the right.

Is there some sort of algorithm to achieve this?

As I have to do this which 6000 other streets, I will probably have to write a python script to do this.

Is there a possibility to move curved lines too without converting them to points?

What I have tried so far is a QGIS plugin called "Line Offset" which creates a parallel line of a selected polyline feature at a given distance. I do not know if this plugin solves my problem, but I would like to give it a try. I unchecked all layers in my layer window except the street map layer. Furthermore, I selected (in the attribute table) all objects I want to move. In the plugin's interface, I have the following options set:

But it always stops running and gives me the following error:

Maybe someone has already worked with this plugin and can figure out what I am doing wrong. It says "QgsRasterLayer" object has no attribute.. but isn't my streetmap layer a polyline-layer?

This doesn't make sense to me.

I have selected both the layer and the objects of the layer, but I still get the same error.


texts: What could cause a Spatial Join (Completely Within) to return false positives?

texts:
Im using ArcMap 10.3

I have a large feature class of Archaeological surveys that I need to Spatial Join to a target feature class of buffered work locations to find locations that are Completely Within a survey. This is part of a large model I have created for our team, but my coworker noticed it is returning some join features that do not contain the target features. What's more strange, is that when I switched the join option to just Within, which should be less limited in scope, these false positives disappear, although now I am left with other unwanted target features where work areas and surveys just partially overlap.

What could cause this strange error?

I have run the repair geometry tool on both feature classes involved. Are there other geometry errors that could be causing this?

The image above shows the two work areas that are falsely returning information as if they were entirely within the survey area in orange.


texts: Rescaling NDVI values 0-255 to [-1;1] using ArcGIS Raster Calculator?

texts:
How do I rescale NDVI values that are 0-255 to -1;1?

I'm trying to do it in ArcGIS Desktop.


texts: Help finding point in polygon via Spatialite with example karst zones data

texts:
I'm working on a school project involving Spatialite and point in polygon queries.
I've successfully followed this example:
http://blog.crazystyle.net/?p=255

I am now trying the same sequence of steps with the shapefiles for “Travis and Williamson Counties Karst Zones and Critical Salamander Habitat”  downloadable from here and viewable online here. I have manually determined that the GPS point 30.418986, -97.83754 should intersect the polygons Jollyville Plateau salamander subsurf crit. hab., karst zone 3, and Jollyville Plateau salamander surface crit. hab.

However, I cannot get the SQL query I'm running to return the data I expect. Please see the transcript below for the query syntax.

Am I running the correct query for this set of data?
Does the shapefile have the info I need for this sort of query? I notice the column names for this data are different from the column names in the example census data table.

Please note: I am entirely new to GIS. My project is a software project using GIS data/technology.

Here's a transcript of my spatialite work:

spatialite> spatialite karst.sqlite
SpatiaLite version ..: 4.3.0    Supported Extensions:
    - 'VirtualShape'    [direct Shapefile access]
    - 'VirtualDbf'      [direct DBF access]
    - 'VirtualXL'       [direct XLS access]
    - 'VirtualText'     [direct CSV/TXT access]
    - 'VirtualNetwork'  [Dijkstra shortest path]
    - 'RTree'       [Spatial Index - R*Tree]
    - 'MbrCache'        [Spatial Index - MBR cache]
    - 'VirtualSpatialIndex' [R*Tree metahandler]
    - 'VirtualElementary'   [ElemGeoms metahandler]
    - 'VirtualXPath'    [XML Path Language - XPath]
    - 'VirtualFDO'      [FDO-OGR interoperability]
    - 'VirtualGPKG' [OGC GeoPackage interoperability]
    - 'VirtualBBox'     [BoundingBox tables]
    - 'SpatiaLite'      [Spatial SQL - OGC]
PROJ.4 version ......: Rel. 4.9.1, 04 March 2015
GEOS version ........: 3.4.2-CAPI-1.8.2 r3921
LWGEOM version ......: 2.1.8
TARGET CPU ..........: x86_64-linux-gnu
the SPATIAL_REF_SYS table already contains some row(s)
SQLite version ......: 3.8.11.1
Enter ".help" for instructions
SQLite version 3.8.11.1 2015-07-29 20:00:57
Enter ".help" for instructions
Enter SQL statements terminated with a ";"
spatialite> .loadshp /mnt/Trav_Will_KarstZones_Veni_07 zones CP1252 26914
========
Loading shapefile at '/mnt/Trav_Will_KarstZones_Veni_07' into SQLite table 'zones'

BEGIN;
CREATE TABLE "zones" (
"PK_UID" INTEGER PRIMARY KEY AUTOINCREMENT,
"OBJECTID" INTEGER,
"ZONE_" INTEGER,
"ZONE_DESC" TEXT,
"LABEL" TEXT,
"MAP_ID" DOUBLE,
"Area_km2" DOUBLE,
"Shape_Leng" DOUBLE,
"Shape_Area" DOUBLE);
SELECT AddGeometryColumn('zones', 'Geometry', 26914, 'MULTIPOLYGON', 'XY');
COMMIT;

Inserted 310 rows into 'zones' from SHAPEFILE
spatialite>  select * from zones where within(GeomFromText('POINT(-97.837543 30.418986)'),zones.Geometry);
spatialite>


texts: Subtract two daily rasters based on name using ArcGIS ModelBuilder?

texts:
I'm not creative today it seems and I haven't found a way to do the following: I have two folders each containing 10,000+ raster files for daily values with names like 19810101_A and 19810101_B for the first day of January 1981 and so on. Now I want to substract the A-raster files in one folder from the B-ones in the other folder on a day-by-day basis and write them to a third folder.

How can I do this in ModelBuilder?


texts: arcgis javascript api map spatialreference changes when i add a layer

texts:
This one is strange: i have a map that takes features from a gp service and turns them into a graphics layer. If i add the layer, for some reason the map's spatialreference changes from 4686 to 4326, if i comment that line, the change does not happen.

The features exposed by the service come in 4686, below is the relevant part

Even worse, the table where data is stored, was created via arccatalog wid wkid 4686. So, i have absolutely no idea where this 4326 is coming from.

<!DOCTYPE html>

    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
    <meta name="viewport" content="initial-scale=1, maximum-scale=1,user-scalable=no"/>
    <title>test</title>

    <link rel="stylesheet" href="http://js.arcgis.com/3.12/dijit/themes/nihilo/nihilo.css">
    <link rel="stylesheet" href="http://js.arcgis.com/3.12/esri/css/esri.css">


texts: Looping through raster layers in folder with PyQGIS?

texts:
I'm looking to apply a mask to a number of rasters in PyQGIS.
(updated but still not working)

import processing
rstr = "/Users/rasterfilepath/"
shp = "/Users/shapefile.shp"
outputDir = "/Users/output/"

for lyr in rstr:
  processing.runandload("gdalogr:cliprasterbymasklayer", rstr, shp, "none", False, False, "", outputDir + lyr + ".tif")

What am I doing wrong?


texts:
I have two tables. One was created in excel, Table1, the other is an attribute table associated with a shapefile, Table2.

I want to join Table1 and Table2 into a single table, using the day of the year. When I try and run the "Join Field" tool, Table3 contains all the information from Table1, but Table2's time values have all been changed to midnight. When I perform a temporary join, it works just fine. I tried turning Table1's 'Date' from a string to a date, so it would match Table2, but it didn't fix anything.


texts:
Creating QGIS plugin via PluginBuilder and compiling resources via simple 'make' works straightforward under Linux (SuSe 13.2).
Under Windows 'make' (3.81) produces an error:

C:\Users\Jochen\.qgis2\python\plugins\ALBa>make
pyrcc4 -o resources.py  resources.qrc
process_begin: CreateProcess(NULL, pyrcc4 -o resources.py resources.qrc, ...) failed.
make (e=2): Das System kann die angegebene Datei nicht finden.
make: *** [resources.py] Fehler 2

[The system can not find the file specified]

What is my problem (besides using Win...) ?


texts:
I need to compare the values of two raster layers at a set of given points. The raster layers are at a national scale, split into tiles and are at 1m resolution (2+TB). The points number in the 100s of millions.
I have an implementation which works but is slow and consumes a lot of resources. Two virtual rasters of the layers are created and the points are split into ~1 million point chunks (arbitrarily selected). Rasterio is then used to sample the rasters at each point and the results (in vector format) are compared. The process takes about 20 mins / million points and tops out at about 10GB of RAM. I've parallelised the process but the limiting factor seems to be the RAM requirement. The CPU load is relatively low and the disk I/O does not appear to be a limiting factor.
The raster layers are already loaded in a local GeoServer. Using this to interface with the data should mitigate the RAM requirement and enable more parallel processes to run concurrently. However, I have been unable to find a way to sample a raster at a given point rather than downloading an area GeoTIFF.
I can go down the GeoTIFF route if that is the only option but it seem counterproductive to speeding up the process due to the fact that it must be written to disk first (even if its just a tempfile). So, I have two questions:

What is the fastest way to sample a raster (of this size) at a given set of points
Can it be done using a GeoServer WCS

NB This process is taking place on a scalable GCP VM. If required it can be scaled up to 96 vCPUs and/or 624GB of RAM but I'd prefer not to brute force it and find a more elegant way. The scaling could apply to the GeoServers specs too.


texts: Problem with calculating QGIS Zonal Statistics MEAN

texts:
I need to calculate some zonal statistics from raster files for some polygons in a shapefile. I am using the following code which works:

    eligibleLayer = QgsVectorLayer("/home/usr/Desktop/polygons.shp", "polygons", "ogr")
    for raster in findRasters(windMapsDirectory, '*.tif'):
        (infilepath, infilename) = os.path.split(raster)
        windMapName = infilename
        list1 = re.split('\W+',windMapName)
        list2 = re.split("_", list1[0])
        list3 = list2[2]
        string1 = ''.join(list3)
        if string1 == 'K':
            zoneStat = QgsZonalStatistics(eligibleLayer, windMapName, 'K_', 1)
            zoneStat.calculateStatistics(None)

The problem is that I only need the mean from the zonal statistics (no count and sum). I tried using:

        if string1 == 'K':
            zoneStat = QgsZonalStatistics(eligibleLayer, windMapName, 'K_', 1, QgsZonalStatistics.Mean)

However, I get the following error: type object 'QgsZonalStatistics' has no attribute 'Mean' . I am using QGIS 2.8.6.


texts: Converting xls into csv in QGIS with correct column definitions and exporting the data into postgis

texts:
I'm using QGIS 3.10.2 A Coruna.
After reading this tutorial (https://ieqgis.wordpress.com/2015/02/08/importing-csv-files-into-postgresql-using-the-db-manager-in-qgis/) I thought I had finally found a simple way of transforming xls-Data in csv and then loading this csv-Data into Postgis. All actions taken before in Loading data directly into pgAdmin were very straining because there were always problems with importing the data with the correct column defitions. When converting xls into csv via qgis all columns are converted into strings, which is incorrect.
Are there any thoughts how excel sheets can be transformed to csv with correct column definitions and imported with these correct formats into Postgis?


texts: Plotting polyline trail using slope raster with least change in slope and ArcMap

texts:
How do I plot a trail in ArcMap so that it runs from point A to point B?
This is part of a larger model I am working on that takes various variables into account and creates a cost layer that I then can use to run a cost path analysis to determine the "cheapest" path between 2 points. Up until now, slope hasn't been too important. I am just now working with slope in the model, and I am trying to think of a method that would create a preference map of sorts to influence the trail to be plotted perpendicular to the slope. It is fine if the trail has a change in elevation, as that is pretty typical, but I want the slope to influence the trail and cause the generated trail to more closely conform to be perpendicular to the slope.


texts:
I have 20 lines on one layer, and many of them overlap like a grid.

I have currently used the following:

[%'ROUTE #' || CASE WHEN $id = $atlasfeatureid THEN "Route #" END%]

When I use the atlas function to zoom to a specific line, it will still show the other lines within the zoom extent.

I only want the line that the atlas function is zooming to to be visible, but still allow me to flick through all 20 lines/atlas pages.

Is there any way to do this?


texts:
I've turned a shapefile into a postgis table using shp2pgsql-gui.

Then, i wanted to get the coordinates of geom column, for that i converted geom to WTK throught ST_AsText(geom) function.

The problem is, the coordinates format is unknown for me :

-69674.5303999996,-294276.2939

Anyone?


texts: MapServer PostGIS: WMS Layers Without Features

texts:
I've prepared an instance of PostGIS and MapServer on the same machine. I have loaded data into PostGIS and can view these features successfully in thick clients (ArcGIS Pro, QGIS).
My MapServer appears to be creating a successful WMS service with two layers. However, these two layers, called from PostGIS, are empty of any features. The same two layers were the features successfully confirmed to be in the PostGIS database, above. Nevertheless, my WMS layers from MapServer show up empty.
What am I missing?
My Mapfile:
https://pastebin.com/tRTJKTxA
My ms_error.txt:
https://pastebin.com/daJJMEkJ
My PostgreSQL log file:
https://pastebin.com/EN2LzheE
My WMS Returns:
https://pastebin.com/NR0W0Pg7


texts: What's this JSON format when generically connecting to Mifi/Jetpack GPS over Wi-Fi (NMEA)?

texts:
I just bought a Mifi 4G wifi modem, which Verizon, a US mobile company, has rebranded as a 'Jetpack'.

While playing with the settings, I found that it also has a GPS receiver built into it, and offers what it calls 'GPS over Wi-Fi (NMEA)'. It comes with a driver for connecting to Windows' GPS system, but no drivers for other operating systems. An extensive session of websearching has not found any Linux tools for connecting to this device, nor any information about connecting gpsd, the Linux GPS daemon, to an NMEA source via wifi.

But I've done a bit of hacking and have found that I can curl the current GPS data in JSON format:

$ curl http://192.168.1.1/gps/status/
{"gBaseAssetsPath":"09101217","gHtmlLang":"en-US","gIsLoggedIn":1,"gLangDirection":"ltr","gLocale":"en_US","gPageId":"gps_status","gSecureToken":"0123456789abcdef0123456789abcdef01234567","gpsIsOff":0,"gpsPrivacyCheckbox":1,"gpsStartupState":0,"gpsStatus":2,"gpsStatusAltitude":"1733 m","gpsStatusHeading":0,"gpsStatusLatitude":37.76520123,"gpsStatusLongitude":-112.24160123,"gpsStatusSatelliteCount":0,"gpsStatusTimestamp":435803696,"gpsStatusUncertainty":"10 m","haveInternet":1}

That's certainly good enough for me to write code to interpret it, but it'd be nice if I could get gpsd to connect to it as well, so I can use other Linux GPS tools.

Is 'GPS over Wi-Fi (NMEA)' a rebranded name for some protocol that I could websearch for?

I'm far from an expert, but that JSON format doesn't look remotely like NMEA, but maybe it's some other standard format, and I could websearch for that.

Any Linux advice about piping lat/long into gpsd so it's available to other apps would be welcome too, of course.


texts: Count rows of attribute table of a shapefile in ArcGIS with Python

texts:
I need to count the rows of a certain selection of the attribute table of a shapefile. I am very new to Python and have never executed it in ArcGIS.

When I run this code in the Python console of ArcMap

import arcpy

table = "C:/path/name.shp"
cursor = db.cursor()
count = cursor.execute("SELECT COUNT(*) FROM " + str(table) "where y = 0 and x = 0")
db.commit()
print cursor.fetchone()[0]
cursor.close()

I get the error Parsing error SyntaxError: invalid syntax (line 5)

If I run

import arcpy

table = "H:/Bearbeitung IS_SaND/Kartenbearbeitung/Kartendaten/mas_p.shp"
cursor = db.cursor()
count = cursor.execute("SELECT COUNT(*) FROM %s where y = 0 and x = 0" + str(table))
db.commit()
print cursor.fetchone()[0]
cursor.close()

I get error

Runtime error
Traceback (most recent call last):
  File "<string>", line 4, in <module>
NameError: name 'db' is not defined

How can I fix my code?


texts:
I'm trying to implement a WFS service in jsFiddle for an OpenLayers map.

It works fine here:
  http://openlayers.org/dev/examples/cql-format.html

My example of it is (the map loads but the filter doesn't):
  http://jsfiddle.net/L7sFY/

But I can't seem to get the filter to work properly.

WFS:
  http://en.wikipedia.org/wiki/Web_Feature_Service

XML WFS file:
http://demo.opengeo.org/geoserver/wfs?service=WFS&version=1.0.0&request=GetFeature&typename=topp:states&outputFormat=json&format_options=callback:loadFeatures


texts: Seeking data for urban and rural areas of United Kingdom?

texts: Is anyone aware of any (commercial) data sets that divide the UK into urban and rural areas based on admin data (e.g. outcode postcodes) and/or physical data (e.g. eastings and northing or lats and longs).


texts:
The WMS call GetLegendGraphic (in my case to Geoserver) returns a single image of icons and legends. How can I specify the width of the returned image? The WIDTH parameter is the width of the symbology icons, not the returned image.
The aim is to force word-wrapping of legend text for long layer names, so I hope that if the width can be specified, text will be wrapped.
I know I can get individual legend icons back via the RULE parameter, but would rather control the whole legend block if possible.


texts:
I downloaded part of a map using JOSM, and saved it in .osm. I am reading the same file using networkx in python. When I read the attribute of an edge, it gives me something like this:

{'ShpName': 'streets_all', 'Json': '{ "type": "LineString", "coordinates": [ [ 55.4509595, 25.2316337 ], [ 55.4508415, 25.2311484 ] ] }', 'Wkb': '\x00\x00\x00\x00\x02\x00\x00\x00\x02@K\xb9\xb9\nx)\x07@9;LX\x9e&\xc7@K\xb9\xb5,\x9d\x16\xfd@9;,\x8a\xa2\x85\xd1', 'id': 169108307, 'Wkt': 'LINESTRING (55.4509595 25.2316337,55.4508415 25.2311484)'}

My code to read the above is:

G_streets = nx.read_shp("D:\\osm_maps\\streets_all.shp", simplify=False)
for n1,n2 in G_streets.edges_iter():
    print G_streets.edge[n1][n2]

I would like to read information related to streets like: maxspeed, oneway, highway, name, ref. In JOSM, all of those are stored as tags. These tags do not appear on the above structure. Am I doing the right thing? and how can I read these tags?

UPDATE:
Please note that I use QGIS to import the .osm file and convert it to shp file.


texts: Multiples labels one on top another on bottom

texts: How can i do this labels? A labels in top another in bottom in the same layer?


texts:
I am trying to style an GeometryCollection geometry data type from PotGIS table in Geoserver. I have not found any examples on the Internet how to do it nor is it supported in Geoserver or not. When I try to style it using standard PolygonSymbolizer, LineSymbolizer and PointSymbolizer I get strange behaviours. For example, filled non-closed linestrings, non-exitsting points appear on linestring etc. This is one example of my SLD file and rendered WMS request.

Am I doing something wrong in my SLD?

Is GeometryCollection supported in Geoserver?

Example of SLD:

<?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>
       <Rule>
        <ogc:Filter>
         <ogc:PropertyIsEqualTo>
            <ogc:Function name="geometryType">
               <ogc:PropertyName>the_geom</ogc:PropertyName>
            </ogc:Function>
            <ogc:Literal>GeometryCollection</ogc:Literal>
         </ogc:PropertyIsEqualTo>
       </ogc:Filter>
          <PolygonSymbolizer>
            <Fill>
              <CssParameter name="fill">#FF0000</CssParameter>
            </Fill>
            <Stroke>
              <CssParameter name="stroke">#FFFFFF</CssParameter>
              <CssParameter name="stroke-width">6</CssParameter>
            </Stroke>
          </PolygonSymbolizer>
          <LineSymbolizer>
            <Stroke>
              <CssParameter name="fill">#00FF00</CssParameter>
              <CssParameter name="stroke-width">2</CssParameter>
              <CssParameter name="stroke-linejoin">round</CssParameter>
              <CssParameter name="stroke-linecap">round</CssParameter>
            </Stroke>
          </LineSymbolizer>
         <PointSymbolizer>
              <Graphic>
                <Mark>
                  <WellKnownName>circle</WellKnownName>
                  <Fill>
                    <CssParameter name="fill">#0000FF</CssParameter>
                  </Fill>
                </Mark>
              <Size>9</Size>
            </Graphic>
          </PointSymbolizer>
      </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>

Rendered WMS request with provided SLD file:


texts:
New to GDAL, have some tiff files and I'm trying to do the following:

"0 to NoData" to avoid black background in Mapbox:

gdalwarp -dstnodata 0 -dstalpha [input file] [output file]

The problem is, the output file size for this tasks is huge (≈150MB to ≈10GB).

The same thing occurs when I next add a color relief:

gdaldem color-relief -of GTiff [input file] [color text file] [output file]  -alpha

The file size roughly doubles...

My only guess is that the data is remaining as a value and is somehow being interpreted as such, so my only option is to compress the file using LZW afterwards.  However, this process takes a ridiculously long time to complete, so I'd like to find a more efficient way.  I typically don't see such great file inflation, but this dataset is new to me.

EDIT:

Here is the gdalinfo readout...

Driver: GTiff/GeoTIFF
Files: [input file]
Size is 129753, 78920
Coordinate System is:
PROJCS["NAD_1983_Albers",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.2572221010042,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Albers_Conic_Equal_Area"],
    PARAMETER["standard_parallel_1",29.5],
    PARAMETER["standard_parallel_2",45.5],
    PARAMETER["latitude_of_center",23],
    PARAMETER["longitude_of_center",-96],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-2353941.226999999955297,3160112.077999999746680)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  DataType=Generic
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-2353941.227, 3160112.078) (127d48'20.50"W, 47d51'28.35"N)
Lower Left  (-2353941.227,  792512.078) (119d57'48.72"W, 27d30'44.36"N)
Upper Right ( 1538648.773, 3160112.078) ( 74d45'30.86"W, 49d54'18.35"N)
Lower Right ( 1538648.773,  792512.078) ( 80d 8'39.87"W, 29d 2'53.17"N)
Center      ( -407646.227, 1976312.078) (100d52' 1.48"W, 40d42'15.31"N)
Band 1 Block=128x128 Type=Byte, ColorInterp=Gray
  NoData Value=0
  Metadata:
    RepresentationType=THEMATIC
  Image Structure Metadata:
    NBITS=2


texts: Appending data to SQL table in QGIS by matching field names

texts:
In QGIS 3+, is there a plugin that allows you append data from one feature into another without needing to move columns around, and map columns automatically based on field name?
I want to merge between two layers where both layers have the same column names, e.g. Feature 1: Length into Feature 2: Length. We currently change the order of columns of sourch feature to match the destination feature.


texts: Calculate Field not saving on joined layer

texts:
I'm joining a table to a layer file, calculating a field in the layer and then removing the join. The code below runs with no errors, but when I check the field in the attribute table, nothing has been updated. If I do this manually in ArcMap with Field Calculator, all is good. What am I missing?

checkds = 'AllCounties' # don't want to calculate values for this feature class
fields = {'FMCDL' : 'CityLeft', 'FMCDR' : 'CityRight'}
lyr = 'Working'
join_table = os.path.join(targetGDB, 'MunicipalCodes')
arcpy.MakeTableView_management(join_table, 'tblvw')
join_field = 'Code'
calcFrom = '!MunicipalCodes.Description!'
fds = 'Framework'
fclasses = getFeatureClasses(fds)

if fclasses: # make sure feature classes were returned
    if checkds in fclasses:
        fclasses.remove(checkds) # remove unwanted feature class

    for ds in fclasses: # each feature class
        working_lyr = ds #+ '_' + lyr # layer name must be unique for each feature class
        arcpy.MakeFeatureLayer_management(ds, working_lyr)

        for k, v in fields.iteritems(): # each  city left and right field
            try:
                print('\tAdding join...')
                arcpy.AddJoin_management(working_lyr, k, 'tblvw', join_field)
                arcpy.CalculateField_management(working_lyr, v, calcFrom, 'PYTHON_9.3')
                arcpy.RemoveJoin_management(working_lyr)


texts: Geoserver 2.4.0 via the web admin interface, not adding layers in Internet Explorer, but working in Firefox

texts: For some reason I cannot seem to add layers via Internet Explorer, but I can via Firefox 23. In addition the styles that I upload also appear broken in Internet Explorer 8, but do not in Firefox. Is there any sort of modification to enable it to work in any browser?


texts: Set the layers style programmatically using QGIS

texts:
I have a layer of GPS points (about 60 layers) loaded inside my .qgs file.
I want to set all their styles in the same manner using a .qml file without doing it one by one.

Can somebody illustrate initial Python code to do this? I'm relatively new to QGIS (2.16.1 Nodebo) since I came from ArcGIS doing simple batch processing through arcpy. In arcpy, I call the arcpy.mapping module to access my .mxd file and update my layer style through UpdateLayer function.


texts: How come when I try and export my time series in ArcMap it compresses the video to 3 seconds? My time change visualization is at least 10. What can I do to make my exported video longer?