texts:
I have generated a coarse heatmap of where data is available by 1 degree x 1 degree cells that I want to overlay on top of a world image.  Using a python script I merged two images (heat map and a world map -- generated with gdal_rasterize).

Reviewing the picture I have a vertical registration issue (I am guessing a flattening issue).   The picture can be sent if needed but I cannot post it at this time (I do not have a reputation >= 10).

How can I use gdal to merge these images correctly?

Thanks


texts: Scripting Address Coder actions in ArcPy?

texts:
I am currently using ArcGIS 10.3 for Desktop that comes with the "Address Coder" tool. This tool has some features that I need as part of my geocoding requirements: in specific, the option to append FIPS codes to the output. However, I see that this option is not available if I use the geocoding feature directly from ArcMap or ArcPy (arcpy.GeocodeAddresses_geocoding function).

Is there a separate scripting API (through ArcPy) available for Address Coder so that I can use its additional capabilities?


texts: calculate positive difference by column common value with sql

texts:
I have one table with a few columns (and around 9000 rows) see below. Now I would like to calculate the positive difference from R where Xn has the same value and when Xn has the same value as Yn.

So the search should take place within the column xn first where there is a value the same like in this example, value 1 then the R-diff should be calculated. After that the xn-value should be compared with the yn-value, when these values are the same the R-diff should be calculated.
Is this possible with sql, i tried a few things but so far could not figure it out...?


texts: How to detect coordinate reference system (CRS)? – I cannot display two vector layers in QGIS at the same time!

texts:
When I in QGIS try to display layer (a) together with another layer (b) that differ in their CRS it doesn’t work. Both are vector shapefiles. I have tried both with using QGIS ‘on the fly’ CRS transformation and by saving the two different layers in the CRS- format of the other layer. How do I display both layers at the same time?

About layer (a)

I have gotten a .shp file sent to me (together with .dbf,.sbn, .shx, and in addition .sbx, .TAB, .idm and .ind files – which I don’t really know what it is). There is no documentation, but I really want to use the map! QGIS displays it nicely as CRS: WGS 84, EPSG:4326. It maps Sweden. It’s metainfo is:

xMin,yMin 4062955.57,3164156.34 : xMax,yMax 4646582.49,4704921.56

About layer (b)

I trust the other layer (b) more- it is documented to be SWEREF99 TM, EPSG:3006 with metainfo:

xMin,yMin 269616.42,6137945.67 : xMax,yMax 749134.54,6908654.00

More on what happens

I cannot convert layer (a) to EPSG:3006 since I get an error message. When I convert layer (b) to EPSG:4326 it becomes weirdly squeezed, and QGIS don’t want to display the layers at the same time. Their coordinates differ heavily, so I guess it is not strange that QGIS don’t want to display them at the same time. I’m using QGIS 2.8.2 on Mac.

Speculation on "what’s wrong?"

My first guess is that I have not selected the right CRS for layer (a)- but how do I know what it is? I am completely new to all this.


texts:
I'm using QGIS 3.4 to consult the Italian cadastre maps, using the Government's WMS service available here: https://wms.cartografia.agenziaentrate.gov.it/inspire/wms/ows01.php

As the technical specifications of the WMS source explain (confirmed by GetCapabilities), the native coordinate system of their backend is EPSG:6706, so, at first, that's what I tried to use for my map. I set that as my project's CRS, enabled on-the-fly reprojection, imported some OSM tiles for visual reference (EPSG:3857).

I then created the WMS connection in QGIS and added the relevant layers. The outline of Italian provinces appeared; however, it overlayed Eritrea. Thinking I had made a mistake, I checked the CRS for the OSM tiles was correct, and it was.

So I made a test: I first checked that QGIS was making requests to the service using EPSG:6706. Right-clicked on the layer, clicked on Properties, and the WMS parameters confirmed the right CRS: crs=EPSG:6706&dpiMode=7&format=image/png&layers=province&styles&url=https://wms.cartografia.agenziaentrate.gov.it/inspire/wms/ows01.php.

I then tried to add the layer using a different CRS, just for the sake of trying. I opened the Data Source Manager, selected the WMS data source, clicked Connect, selected the same layer as before, changed the CRS to EPSG:25832 and added it. To my surprise, it worked! The outlines now overlayed the OSM tiles for Italy. Of course, the lines were a bit blurred, because my project CRS was still set to EPSG:6706 and so the WMS data was being reprojected on the fly. So I tried setting my project CRS to EPSG:25832 as well and everything looked fine.

At that point, I wrote the issue off as a mistake on the part of the WMS service, maybe they were simply returning the wrong data when making requests with EPSG:6706. That would make sense, because the issue was solved by simply switching to a different supported CRS.

Then I started fumbling around, and at some point I needed to find the GPS coordinates (i.e. EPSG:4326 coordinates) of something I was seeing on the map. So I noted the shown X,Y coordinates (in EPSG:25832, the project CRS), I opened a calculator, I set EPGS:25832 as the source CRS and EPSG:4326 as the destination CRS and I typed in the coordinates.

The coordinates didn't quite look right to me, so I pasted them in Google Maps and I ended up in... Eritrea.

I'm now convinced I'm doing something wrong; however, I simply do not know what to check. If I'm doing something wrong with the coordinates, how can the OSM tiles and the WMS outlines match perfectly, being in two completely different CRSs?


texts: QGIS data is not shown while opening a .qgis file

texts:
I am trying to load a module project with .qgis from my phone from a program called Locus GIS project called "power network", by me at house it loads on my win 10 pc.

But my ubuntu pc at work the data won't load. The layers load already and I can plot new data, points, elements but the saved "network power" won't load.
2021-04-06T16:02:22     WARNING    SQLite-Fehler:no such module: rtree
  SQL: SELECT "ROWID", AsBinary("geometry") FROM "power_lines_110kv" WHERE ROWID IN (SELECT pkid FROM "idx_power_lines_110kv_geometry" WHERE xmin = 10.37316050898623132 AND ymin = 50.98134542764638866)
2021-04-06T16:02:22     WARNING    SQLite-Fehler:no such module: rtree
  SQL: SELECT "ROWID", AsBinary("geometry") FROM "power_lines_110kv" WHERE ROWID IN (SELECT pkid FROM "idx_power_lines_110kv_geometry" WHERE xmin = 10.37316050898623132 AND ymin = 50.98134542764638866)
2021-04-06T16:02:22     WARNING    SQLite-Fehler:no such module: rtree
  SQL: SELECT "ROWID", AsBinary("geometry") FROM "power_lines_380_kv" WHERE ROWID IN (SELECT pkid FROM "idx_power_lines_380_kv_geometry" WHERE xmin = 10.37316050898623132 AND ymin = 50.98134542764638866)
2021-04-06T16:02:22     WARNING    SQLite-Fehler:no such module: rtree
  SQL: SELECT "ROWID", AsBinary("geometry") FROM "power_lines_380_kv" WHERE ROWID IN (SELECT pkid FROM "idx_power_lines_380_kv_geometry" WHERE xmin = 10.37316050898623132 AND ymin = 50.98134542764638866)
2021-04-06T16:02:22     WARNING    SQLite-Fehler:no such module: rtree
  SQL: SELECT "ROWID", AsBinary("geometry") FROM "power_poles" WHERE ROWID IN (SELECT pkid FROM "idx_power_poles_geometry" WHERE xmin = 10.37316050898623132 AND ymin = 50.98134542764638866)
2021-04-06T16:02:22     WARNING    SQLite-Fehler:no such module: rtree
  SQL: SELECT "ROWID", AsBinary("geometry") FROM "power_poles" WHERE ROWID IN (SELECT pkid FROM "idx_power_poles_geometry" WHERE xmin = 10.37316050898623132 AND ymin = 50.98134542764638866)


texts: Preserving Projection on SAGA-produced Rasters

texts:
I perform a fair amount of hydrological processing in QGIS and rely on SAGA's Topographical Wetness Index and slope-derived functions. In order to use them, SAGA demands that input rasters are warped from geographic to projected "meters-on-the-ground".

The problem arises from the fact that a successful operation like "Fill Sinks", "TWI" or "Relative Heights and Slope Positions" always produces a raster that QGIS loads and interprets as geographically projected (e.g. WGS 84).

I've tried setting my project to the projected CRS prior to running any processes but it seems that rasters produced in SAGA need to have their projections manually set in order to be loaded properly by QGIS. I find this particularly strange since SAGA operations typically call for projected rasters to begin with, so loading SAGA-produced rasters as geographically projected seems counterintuitive.

Am I missing an option that should be set within QGIS to enable the creation of a CRS-defining file or do SAGA-produced rasters need to always have their CRS's changed manually?


texts: Compression not applying when updating/exiting GeoTiff file in GDAL?

texts:
I have a GeoTiff file with dimensions 300 x 300 and 8760 bands (hourly for one full year). I am testing some code that opens the file, updates ONE band and exits - but I am finding that the file size almost doubles when I exit the file, despite the new band being the same data type etc. Code is as follows:

I am creating the GeoTiff file using the following creation options using many singular GeoTiff files:

gdal_merge.py -separate -co BIGTIFF=YES -co COMPRESS=DEFLATE -co PREDICTOR=1 -co TILED=YES -co BLOCKXSIZE=64 -co BLOCKYSIZE=64 -o block64.tif *.tif

The file size is 947MB.

I am then updating the file using the following code:

ds = gdal.Open("block64.tif", gdal.GA_Update)

# Some simple operation to change band values
arr = ds.GetRasterBand(11).ReadAsArray() / 2
ds.GetRasterBand(11).WriteArray(arr)
ds.FlushCache()
ds = None

The file size is now 1.9GB.

My preliminary conclusion is that the compression/creation options are not being applied once I've updated and exited the file. Is there any way to re-compress/reapply the options once I close the file?


texts: PyQGIS: How to color raster layer with existing colorBrewerMap

texts:
With pyqgis I have added multiple new raster layers to the legendinterface().
Because i like the ColorBrewer style i want to color some of them with the "Spectral" gradient.

How can I do this in python? The API-help only gives advice how to color with individual set colorramps.
I found out that I can query the existing Ramps with QgsColorBrewerPalette.listSchemes(), but how can I set this style for my layer?

qshade = QgsColorRampShader()
qshade.setColorRampType(QgsColorRampShader.INTERPOLATED)
qshade.setColorRampItemList("Spectral") # Similar to this?


texts:
I'm using a walkthrough for QGIS that requires I use the Map Management tool to open a zip file. The walkthrough makes it seem that this is a window that opens automatically when QGIS opens, which doesn't happen for me. I'm using the US Census Bureau's GUPS application.
This is the first time I've ever used QGIS... is there anyone that can give me some guidance to locate the Map Management tool? I can clarify/re-phrase this if necessary.


texts: How to create a spatial view when allowing multiple selection in qgis’s edit widget “value relation”?

texts:
I've a point layer in qgis and serveral non spatial tables in a spatialite database. In qgis the fields were filled using the edit widget "value relation" (multiple selection allowed for some fields). I also created a custom form and conected some fields with a "List widget" in QT-Designer as you can see here, where one can choose more than one entry for a field.
Now in my "master table (baumkataster)" I got the following entries for one point feature (see example below):

baumkataster:
baum_ID, boden_verd
1,{2}
2,{1,3}
3,{1,2,3}

In my query I joined the "master table" with a non spatial table (boden_ver):

boden_verd:
boden_verd_ID,boden_verd
1,value1
2,value2
3,value3

By now I've learned how to create a view to display the values not the key value. For testing I only do the SELECT statement. This works in case of single selection...

For example:

baum_ID,voden_verd
1,value2
2,value1
3,value3

...but it returns a empty table in case of multiple selection.
How can I fix this to get the following result?

baum_ID,boden_verd
1,{value2}
2,{value3,value3}
3,{value1,value2,value3}

I found the following code for a similar problem but it returns an error and I don't realy know whrer I'am wrong:

SELECT baumkataster.baum_ID AS baum_ID,
baumkataster.geom AS geom,

stuff((select  DISTINCT  ', ' +  boden_verd.boden_verd
             from boden_verd
             WHERE ','+baumkataster.boden_verd+',' LIKE     '%,'+boden_verd.boden_verd_ID+',%'
             for xml    path(''),type).value('.','nvarchar(max)'), 1, 2, '' )  AS    boden_verd

FROM baumkataster;

Thanks, Patrick


texts: Adding dates into existing feature class from Excel spreadsheet using ArcPy?

texts:
I've been trying to load data from an excel spreadsheet into existing feature classes, but I've been having some trouble with the dates.

I've written a script that compares the values in the spreadsheet & in the feature class.

arcpy.MakeFeatureLayer_management(engPath, 'engPath')
arcpy.SelectLayerByAttribute_management('engPath', "NEW_SELECTION", d)
EngCursor = arcpy.da.SearchCursor('engPath', '*')
for row in EngCursor:
    for item in Master_Dict.itervalues():
        print item[1], row[3]
        if item[1] == row[3]:
            i = 0
            j = 2
            while i < 34:
                if item[i] == row[j]:
                    print "Values are the same"
                else:
                    print "NOPE!", item[i], row[j]
                i+=1
                j+=1

When I run the script, I get values like this:

Excel Spreadsheet Value = "41526.0"
Feature Class Value = "2013-09-09 00:00:00"

Those two values are supposed to be identical, and the excel value for the cell does show
"09/09/2013".  I've checked the spreadsheet, and the column is set to a date data type.

I don't know how to tell arc that it is indeed looking at a date data type.


texts: How to find the closest stop point PostGIS

texts:
I have Point geometries created from lat,lon
I loaded all the 10,000 stops in a certain city using gtfs

i would like to query my db to find the closest stop reference to my point.

the staypoint table is

    CREATE TABLE stay_points
    (
        ad_id uuid NOT NULL,
        latitude numeric NOT NULL,
        longitude numeric NOT NULL,
        start_utc_timestamp timestamp without time zone NOT NULL,
        end_utc_timestamp timestamp without time zone NOT NULL,
        cluster_id numeric NOT NULL,
        radius numeric NOT NULL,
        geom geometry(Point,4326)
    )
    WITH (
        OIDS = FALSE
    )
    TABLESPACE pg_default;

and the stop table is :

    CREATE TABLE public.stops
(
    stop_name text COLLATE pg_catalog."default",
    stop_lon real NOT NULL,
    stop_lat real NOT NULL,
    stop_desc text COLLATE pg_catalog."default",
    stop_code text COLLATE pg_catalog."default",
    agency integer,
    zone_id text COLLATE pg_catalog."default",
    stop_url text COLLATE pg_catalog."default",
    location_type numeric,
    parent_station text COLLATE pg_catalog."default",
    stop_id text COLLATE pg_catalog."default"
)
WITH (
    OIDS = FALSE
)
TABLESPACE pg_default;

My current solution involves pandas and scipy kdtree and thought there must be a better way


texts:
After starting to use MOBAC to create offline background maps I found the aux.xml files very useful as QGIS automatically sets the CRS so this would be useful for my colleagues with smaller GIS knowledge.

Is it possible to generate such file for particular projection using WKT or .prj file? I am able to get both from spatialreference.org

Sample file from MOBAC for EPSG:3395 is here:

  PROJCS["World_Mercator",GEOGCS["GCS_WGS_1984",DATUM["WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_1SP"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["latitude_of_origin",0.0],UNIT["Meter",1.0]]

I use mostly only few CRS - epsg:4326, 32633, 900913, 3395 - so I could always create a copy of current *.aux.xml and rename it to fit the raster dataset name.


texts:
I did a SQL query from a Geopackage with QGIS database manager. The query is correct and it returns me the expected result with it's ID and the GEOM columns. The problem becomes when I load the query as a layer. The software creates the layer but is void as you can see:

Parameters and columns;

Any ideas?


texts:
I had a discussion about this topic with a co-worker.

We are working implementing a GPS validation system, and I told him the longitude -180 (or 180W) is not a valid longitude because a valid longitude is the same value as one of the angular values of the spherical coordinates, and the latitude is the another one.

  In spherical coordinates:

  Theta must be (-180,180]. Theta does not contain the -180 as a valid value

  Phi must be [-90,90]. Phi contain both values as valid ones

He said to me there must be 2 values of meridian 180 (positive and negative 180 degrees, +180 & -180, or 180E and 180W), and there is no document specifies this issue: [45,-180] is a valid location? Is not the exact position as [45,180] one?

Any online tool I found on the net or code available includes both meridians as different values, but in the real world, are the same!

My question is: Is there any document specifies the 180 meridian can be positive and negative or any document specifies the 180 meridian must be unique? This co-worker does not accept any math demonstrations, therefore any spherical coordinate to rectangular coordinate system will not be a valid reference for this person.


texts:
I am pulling State geometry data from a database (populated from Shapefiles) to find what state a coordinate pair is in.  I was originally using GeoAlchemy, but tens of thousands of queries to determine the state is quite expensive.

def get_state_from_coordinates(self, db_session, long, lat):
    return db_session.query(States)\
        .filter(func.ST_Contains(States.geom, func.ST_PointFromText('POINT({0} {1})'.format(long, lat), 4269)))\
        .first()

So my current solution is to pull all of the state data into my application and do the lookup locally.  I wanted to spike out a lookup and have gotten stuck.

states = db_session.query(States).all()

for state in states:
    point = Point(-81.4406293, 30.323164199999997)
    if state.geom is not None:
        poly = to_shape(state.geom)
        if point.within(poly):
            print('Point is in {0}'.format(state.name))

This should print Point is in Florida but the point is never found in the poly.  For those curious, after transforming the geom from the DB to a Shapely poly here are the poly.geoms: pastebin.


texts: Python Toolbox keep xml meta data formatted for version control

texts:
Is there a way to keep Python Toolboxes metadata formatted nicely for version control?

I use git to maintain my Python Toolboxes and when I make edits to the toolbox the metadata is updated and put all on one line. For example if I make a change to the tool alias I only want git to record the difference on that line.

Currently metadata looks like this

toolbox.tool.pyt.xml

<?xml version="1.0" encoding="UTF-8"?>
<metadata xml:lang="en"><Esri><CreaDate>20180824</CreaDate><CreaTime>10320700</CreaTime><ArcGISFormat>1.0</ArcGISFormat><SyncOnce>TRUE</SyncOnce><ModDate>20180828</ModDate><ModTime>111340</ModTime></Esri><toolbox name="partition" alias="partition"><arcToolboxHelpPath>c:\program files (x86)\arcgis\desktop10.3\Help\gp</arcToolboxHelpPath><toolsets/></toolbox><dataIdInfo><idCitation><resTitle>partition</resTitle></idCitation></dataIdInfo><distInfo><distributor><distorFormat><formatName>ArcToolbox Toolbox</formatName></distorFormat></distributor></distInfo></metadata>

I want this

toolbox.tool.pyt.xml

<?xml version="1.0" encoding="UTF-8"?>
<metadata xml:lang="en">
  <Esri>
    <CreaDate>20180824</CreaDate>
    <CreaTime>10320700</CreaTime>
    <ArcGISFormat>1.0</ArcGISFormat>
    <SyncOnce>TRUE</SyncOnce>
    <ModDate>20180828</ModDate>
    <ModTime>111340</ModTime>
  </Esri>
  <toolbox name="partition" alias="partition">
    <arcToolboxHelpPath>c:\program files (x86)\arcgis\desktop10.3\Help\gp</arcToolboxHelpPath>
    <toolsets/>
  </toolbox>
  <dataIdInfo>
    <idCitation>
      <resTitle>partition</resTitle>
    </idCitation>
  </dataIdInfo>
  <distInfo>
    <distributor>
      <distorFormat>
        <formatName>ArcToolbox Toolbox</formatName>
      </distorFormat>
    </distributor>
  </distInfo>
</metadata>


texts: QGIS function editor tab indents generating errors?

texts:
i'm trying the function editor in QGIS 2.10.1 (windows 32 bits).

I always get "indents" error as the function editor inserts TAB indents instead of spaces ...

How do i solve this issue ?


texts: Which impact has the "layer_coordinate_transform_info"-tag in the QGIS-project file

texts:
In QGIS there is a persistent bug which causes that the project file gets cluttered with old entrys in the "layer_coordinate_transform_info"-tag of layers which are not present in the project file anymore.

https://issues.qgis.org/issues/14733

At the moment I am cleaning old project files and also had to deal with entries which pointed to old destDatumTransform-ID's which did not exist in my srs.db anymore. These entries were also not updated when saving the project.
So I think this are two bugs in QGIS (not removing old entries and not updating entries).

As long as everyone who uses the project file has set the proper options in this section i would expect QGIS to use these settings:

Instead it seems that QGIS uses the entries of the project files "layer_coordinate_transform_info"-section even if the destDatumTransform-ID's used in the project file don't exist any more in the srs.db, so layers are displaced without any error message.

Is this a QGIS-bug? And is it save to just delete this section of the project file?


texts: When debugging ArcGIS add-in, can't use "Edit and Continue"

texts:
I've been working on this add-in for years, and just today set up a new computer with VS 2013 and ArcGIS 10.3.  When I run my add-in from the debugger, breakpoints can be hit, but I can't edit when code execution is paused - I get an error, "Changes are not allowed if the assembly has not been loaded".

I've done all the suggestions I could find in various forums:

Change build output path to bin\

verify that Configuration = Debug

verify that Platform = x86

Edit and Continue enabled (Tools > Options)

checked "Break all processes when one process breaks"

unchecked "Require source files to exactly match the original version"

Anything else I can try?

When I make a simple Windows Forms application, I can successfully "Edit and Continue". –  kev_gis 1 hour ago

I'm not sure how to get an ArcMap SDK installed for Visual Studio 2013. I installed "ArcGIS_Runtime_SDK_DotNet_1025.exe" as well as "ArcObjects_SDK_for_NET_Framework_1022_140467.exe". When I try to make a new Add-In in Visual Studio 2013, I do not see ArcGIS under the list of available templates.

After attempting to install the SDKs listed above, now I get a "No Debugging Information" error message when attempting to run.  "Debugging information for 'ArcMap.exe' cannot be found or does not match. Cannot find or open the PDB file."

I've gone through all of the steps I could find on this page, with no effect:
Why can't the breakpoint be "hit" when debugging an ArcGIS 10 Add-In?


texts: Changing the active tool in an ArcGIS Python Add-in Tool Palette?

texts:
I have an ArcGIS Python Add-in project with a tool palette that contains a dozen different tools. I would like to be able to change the 'active tool', i.e. the tool that is visible on the toolbar, when a button on the same toolbar is clicked.

The active tool changes when the user selects a tool from the dropdown palette, but I'm not sure if this behaviour can be replicated using Python.


texts: How to create a web-based GIS that can be queried by multiple users?

texts:
As part of my position, we're required to randomly select townships in our state, select forest, wetland, or grasslands and then throw random points into those habitats to create site for plant and bird samples.  We've created shapefiles and points in R, but what we'd like to do is create a web-based GIS where any user (with the appropriate security clearance) can select a township, display points and polygons in that township, and the random points contained therein, and generate a map that they can print and take out to the field with them, or share with other workers.

Our base layers would be comprised of roads, land use/ land cover, cities/towns, section lines, etc. and then add our custom layers.  We number each township (1-1765), then give them random ranks for each habitat.  Users would need to select "Habitat Type", "Rank", and "Township" and maybe some other criteria, which would then the GIS would zoom into the township that you selected.

In the past we've done this in a piecemeal fashion where one of us would generate individual maps from a local machine in ArcGIS, but it's pretty ponderous and just having one person do it disrupts workflow.  We'd like to get this set up so that each person could generate custom maps from their own work station (regardless of where that work station is located).

I've looked into GeoServer a little, but I was hoping for a little guidance on good places to start the process.  We're not complete noobs, but as scientists (and not web developers), this is a little out of our comfort zone.


texts: Create Polygons in Python that use the max extents to make a box

texts:
I have field collected data that is usually collected with 6 points to create a form of polygon. I can subset those 6 points and create the polygon no problem, with this code:
eastList = channelFilter['Easting'].tolist()

polygon_geom = Polygon(zip(eastList, northList))
crs = {'init': 'epsg:26919'}
polygon = gpd.GeoDataFrame(index=[0], crs=crs, geometry=[polygon_geom])
print(polygon.geometry)

polygon.to_file(filename = Polygon_Folder + "\\channelGeom.shp", driver = 'ESRI Shapefile'

ChannelFilter is simply the filter of coordinates / points I'm using. However, the order of points collected often goes Bottom middle, Bottom right, bottom left, top middle, top right, top left, meaning the lines cross over each other and give a convex (I think that's the word?) polygon.
Does anyone know a convenient way to take 6 points and create a polygon with those are maximum extents, regardless of the order?


texts: Access QGIS layer style manager in stand-alone pyqgis application

texts: I have a working stand-alone pyqgis application that can display vector and raster layers, the question is how to access the layer style manager as in QGIS's layer property, where users can change the symbology of a layer? I am hoping there is a class already defined in the QGIS core library that I can just call from my application.


texts: Other way to see average standard error in cross validation between model for geostatistical layer?

texts: Is there another way to compare average standard error in geostatistical wizard between model type? It seem not relevant to generate one by one the geostatistical layer to compare each model type just to see the lowest error.


texts: Force coordinate precision in QGIS when saving feature to PostGIS

texts:
Is it possible to configure QGIS in a way, that one can specifically set the desired coordinate precision which is saved for each vertex/point when writing to a PostGIS layer? I want to achieve that WGS84 coordinates are saved as

POINT(0.06288 0.24537)

instead of

POINT(0.0628843267570899 0.24537407891492)

Changing the coordinate display setting of the project does not do the trick.


texts: Indexing attribute field of shapefile in QGIS

texts:
Is it possible to index an attribute value of a shapefile in QGis, and if so save this as part of the project (similair to a spatial index)?

I have a large shapefile of land parcels (maybe about 1 000 000) and I would like to index them on their land parcel number so as to speed up queries on this number.


texts: Where to find topographic basemap with no labels

texts: I am new to QGIS and trying to find a simple topographical basemap for Britain and Ireland. It is a map of prehistoric sites, so I can't have any cities/streets/labels, just topography.


texts: I'm brand new to ArcGIS and I want to create a map for the county conservation board I'm interning at. Specifically right now, I need to add the street names on a map of black hawk county for my boss to take to the board meeting (my other bosses). I have a layer of the actually streets, but I need to label them. I know there is a simple solution to this, but like I said, I am brand new to this world.


texts: Invalid Input Data Error 000368 Running ArcPy tools in for loop

texts:
I am attempting to iterate over all the feature classes in a folder and I'm getting the error 000368. The other solutions in on SE were not helpful to me unfortunately, so I'm posting a separate question. search_grids is a list of point feature classes and searches is a another point feature class.

The points never actually intersect with each other, but I'm using the Select by location tool to select points within a specified distance of each other. This isn't the reason for the error though as I was able to get this to work in the python window in ArcMap.

My code is the following

arcpy.env.workspace = r"C:\workspace"
search_grids = arcpy.ListFeatureClasses()
searches = r"C:\workspace\search"
for grid in search_grids:
    arcpy.management.SelectLayerByLocation(grid,"INTERSECT",searches,"97.5 meters","NEW_SELECTION","INVERT")
    arcpy.management.SelectLayerByLocation(grid,"INTERSECT",searches,"2.5 meters","ADD_TO_SELECTION","NOT_INVERT")

Why am I getting this error?

  arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.
  ERROR 000368: Invalid input data.
  Failed to execute (SelectLayerByLocation)


texts: clip vector by polygons in python script

texts:
I have a big shapefile, and I want to divide in x by y tiles. For testing purposes let's say just 2 by 1. I take the bounds of the big shapefile, create a polygon with smaller dimensions, and try to use ogr2ogr.

            with fiona.open(os.path.join(subdir, file), 'r') as ds_in:
            print('primero', ds_in)
            print('bounds', ds_in.bounds)
            minx, miny, maxx, maxy = ds_in.bounds
            pol = Polygon([(minx, miny), (minx, maxy/2), (maxx/2, maxy/2), (maxx/2, miny)])
            subprocess.call(["ogr2ogr", "-f", "ESRI Shapefile", "-clipsrc", pol, 'results.shp', ds_in], shell=True)
            line = LineString(pol.coords).wkt
            subprocess.call(["ogr2ogr", "-f", "ESRI Shapefile", "-clipsrc", line, 'results.shp', ds_in], shell=True)

I tried with the Polygon object and with a string representation in WKT, noone of them work.

ogr2ogr is suposed to take shapefiles as inputs. Should I create a shapefile with the cutting polygon, and then another shp with the clipped polygons? I doubt that is efficient.

Which is the easiest way of doing this from a standalone python script?


texts: PostGIS automatically creating geometries after importing using QGIS DB Manager

texts:
I have a shapefile consisting of linestrings, in which some of the features have null geometry. Yet, after importing it to PostGIS, using QGIS' DB Manager, a geometry is created for those features (what it does is create a small segment on top of another existing linestring).

Has anyone ever seen this happen? This problem doesn't happen when I import the shapefile using OpenGeo's pgShapeLoader. Is it a bug of DB Manager, or am I not configuring it correctly for importing?

I'm using QGIS 2.2.0 for mac, and PostGIS 10.1.0


texts: Creating own OS Carto Basemap equivalent for use in QGIS

texts:
I am looking to create my own OS Carto Basemap which is currently a ArcGIS third party product. It gives you the ability to to work at different zoom levels within your project.

Is there a way I can mimic this for use in QGIS so that at different zoom levels the basemap refreshes? I understand that you can download the data from ordnance survey but do not know how to structure it to form the basemap.


texts:
I'd like to concatenate a few vrts together. So i have three vrts representing data from different years and to make my life easier would like just to have one vrt created from the three vrts.

Each of the vrts has 365 bands, and i want the outputted vrt to be the stacking of these bands so (3*365) bands.

<VRTDataset rasterXSize="2400" rasterYSize="2400">
  <VRTRasterBand dataType="UInt16" band="1">
   [from 2007.vrt]
  </VRTRasterBand>
    ...
  <VRTRasterBand dataType="UInt16" band="365">
    [from 2007.vrt]
    </VRTRasterBand>
  <VRTRasterBand dataType="UInt16" band="366">
    [from 2008.vrt]
  </VRTRasterBand>
    ...
  <VRTRasterBand dataType="UInt16" band="730">
    [from 2008.vrt]
  </VRTRasterBand>
  <VRTRasterBand dataType="UInt16" band="731">
   [from 2009.vrt]
  </VRTRasterBand>
   ...
  <VRTRasterBand dataType="UInt16" band="1095">
   [from 2009.vrt]
  </VRTRasterBand>
</VRTDataset >

I was hoping something like:

gdalbuildvrt 2007_2009.vrt 200{7,8,9}.vrt

would work. But this unfortunately not. Instead of stacking as above i get:

<VRTDataset rasterXSize="2400" rasterYSize="2400">
  <VRTRasterBand dataType="UInt16" band="1">
   [from 2007.vrt band 1]
   [from 2008.vrt band 1]
   [from 2009.vrt band 1]
  </VRTRasterBand>
    ...
  <VRTRasterBand dataType="UInt16" band="N">
   [from 2007.vrt band N]
   [from 2008.vrt band N]
   [from 2009.vrt band N]
  </VRTRasterBand>
</VRTDataset >

which gdalinfo only records as 365 bands.

Pictorially 2007_2009.vrt would be constructed like:

2007_2009.vrt-|
             -- 2007.vrt
             -- 2008.vrt
             -- 2009.vrt

I'm on a unix system


texts: Large amount of time publishing to ArcGIS Online?

texts:
Publishing <5Mb tile packages to ArcGIS Online is taking more than 20 minutes.  ESRI tech support has consistently given extremely vague answers that never solve and always reference server load.  <5Mb does not seems like much of a server load.

Anyone else having this same problem?


texts: Create letters as geometry polygons from SVG files?

texts: After searching a lot and not finding a definitive direction, is there a way to create geometry for a map that contains letters? I want text that is not a label or a symbol but lays flat as geometry on a map at a GPS location. So I am thinking that I can take SVG outline letters and convert them to polygons. Any advice on how to accomplish this? I want something like this:


texts: How to get elevation of a coordinate given a dtm file out of a PostGIS?

texts:
What would be the best way to get elevation of a coordinate out of a PostGIS database (SQL/plpgsql) given a GeoTIFF dtm file?

Importing as PostGIS Raster (which seems more costly) and calling ST_Value?

Use pgsql-http (https://github.com/pramsey/pgsql-http ) and a cgi which calls gdallocationinfo (GDAL)? Or using pgsql-http to call GDAL/WCS? Other?


texts: Replacing characters using an expression in QGIS

texts:
What regular expression can I use to change values in a string field? I need to change the letter B to a symbol - before the digits, and digits with letter F to make positive.
Like this:


texts: Problem reprojecting coordinates using DotSpatial

texts:
I'm trying to reproject some coordinates from EPSG:27700 to WGS84. I can do this fine in JS using the proj4 port, but I'm getting different results using DotSpatial.

Here is my .NET code:

var wgs84 = KnownCoordinateSystems.Geographic.World.WGS1984;

var osgb1936 = ProjectionInfo.FromProj4String("+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs");

var reproject = new Func<Point,Point>(input =>
{
    var arr = new double[] { input.Coordinate.X, input.Y};
    Reproject.ReprojectPoints(arr, null, osgb1936, wgs84, 0, 1);
    return new Point()
    {
        X = arr[0],
        Y = arr[1]
    };
});

var r = reproject(new Point(161887.0, 1030889.5));

My JS code is:

var osgb1936 = proj4.defs("EPSG:27700", "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs");

var result = proj4('EPSG:27700', 'WGS84', [ 161887.0, 1030889.5 ]);

In JS I'm getting the result [-6.1597047213905, 59.098069153809384] which I believe is correct, but in .NET/DotSpatial I'm getting the result [-162.1702802817330000 -41.6926543605886000].

What am I doing wrong?


texts:
Is it possible to make a single layer style similar to what the map below has?

I.e. first-order districts have different hues, but their constituent second-order districts have different shades.

I'd rather not use the workaround of splitting the 1st-order districts into their own layers, since I'm working on a map with somewhere around 150 districts.

example map: https://commons.wikimedia.org/wiki/File:Hungary_1941-44_Administrative_Map.png


texts: Exporting GRID to TIFF keeping (classified) symbology using ArcGIS for Desktop?

texts:
I want to export an esri grid-raster into a *.tiff without losing the symbology I created for the grid (the symbology is classified).

No matter how I try to solve it I get a stretched symbology in a greyscale. What I've tried so far:

Copy/export raster via: toolbox; TOC, rightclick --> data --> export data. When checking the 'Use Renderer' and 'Force RGB' boxes on the data export menu a lot of values turn into 'noData'. Furthermore, the colors are slightly different.
Created a layer-file from the original raster. Copied the grid and tried to set the source of the layer file to the newly copied raster. However, I could not set the new source.
Imported the symbology using a *.lyr file (Chris W's comment) with the folder button next to the save button upper right of the image: nothing happened.
Imported the symbology directly from my symbolized grid (Chris W's comment): nothing happened.

There must be a way to solve this but I cannot find a solution... Any hints?

Edit: Added clipped grid and *.lyr-file here.


texts: How thick in ground units is this line on a printed map?

texts:
Sometimes in discussing how accurate this that or the other thing needs to be someone will make the observation that "a line shift of 30m at 1:1million scale is less the width of the line itself when printed, so, no worries".  That's great! For that map that is, but what about this map which will be printed at 1:25,000?

This got me thinking that it would be useful to have a quick reference chart of "how thick is an Xpt line on the ground at Y scale at Z dpi?". A bit of searching reveals no such beast readily available online, so let's fix that!

I suggest a base line of 0.2 pt and 300 dpi since in offset printing that's the thinnest reasonable colour line and standard dpi.


texts: Creating feature class on ArcSDE schema with different user

texts:
I have created an enterprise geodatabase using ArcCatalog 10.5 on a PostgreSQL database. I have a database admin user named administrator and another user named sde. The geodatabase admin is sde user.
Geodatabase is created successfully. And I connected to database from ArcCatalog with administrator user. So I want to create a feature class but throws an error:

Failed to create Feature class.
Database user name and current schema
do not match. [mygdb.administrator.test]

Should I connect with sde user and create feature classes or other objects?
So If I need multiple users to create feature objects, all of them will use sde user. Is this right?


texts: Overpass / Overpy: Getting Way IDs from Nodes

texts:
I am struggling to grok Overpass as well as I'd like, so I hope someone might be able to help.

I have a route which I'm importing from a GPX file. This has been 'snapped' to OSM data, so I have the lat/lon of each point on the route. Each of these coincides with an OSM node, but I don't have the OSM node ID in my raw lat/lon data.

I also want to find which of these points are what I'm terming 'intersections', i.e. where the route passes from one OSM way to another. In order to do this, my plan was to go from nodes to (highway) ways and work out where the route changed way, perhaps with some refinements to e.g. take the name of the way into account.

So, using Overpy, I have found that the fastest way to go from the route points to OSM nodes is to use a union of multiple bounding-box queries:

TOL = 0.00001

query = """
<osm-script>
  <union>
"""
for index, row in df_snapped.iterrows():
    query += '    <bbox-query s="{}" w="{}" n="{}" e="{}"/>\n'.format(
         round(row.lat-TOL, 5)
        ,round(row.lon-TOL, 5)
        ,round(row.lat+TOL, 5)
        ,round(row.lon+TOL, 5)
    )

query += """
  </union>
  <union into="_">
    <item from="_" into="_"/>
    <recurse from="_" into="_" type="node-way"/>
  </union>
  <print e="" from="_" geometry="skeleton" ids="yes" limit="" mode="body" n="" order="id" s="" w=""/>
</osm-script>
"""

This results in a query similar to the following:

<osm-script>
  <union>
    <bbox-query s="51.48825" w="-2.62352" n="51.48827" e="-2.6235"/>
    <bbox-query s="51.48801" w="-2.62364" n="51.48803" e="-2.62362"/>
    <bbox-query s="51.4878" w="-2.62373" n="51.48782" e="-2.62371"/>
    <bbox-query s="51.48697" w="-2.62406" n="51.48699" e="-2.62404"/>
    <bbox-query s="51.48682" w="-2.62414" n="51.48684" e="-2.62412"/>
    <bbox-query s="51.4868" w="-2.62416" n="51.48682" e="-2.62414"/>
    <bbox-query s="51.48665" w="-2.62431" n="51.48667" e="-2.62429"/>
    <bbox-query s="51.48654" w="-2.62442" n="51.48656" e="-2.6244"/>
    <bbox-query s="51.48633" w="-2.62463" n="51.48635" e="-2.62461"/>

    ...

  </union>
  <union into="_">
    <item from="_" into="_"/>
    <recurse from="_" into="_" type="node-way"/>
  </union>
  <print e="" from="_" geometry="skeleton" ids="yes" limit="" mode="body" n="" order="id" s="" w=""/>
</osm-script>

The second <union> clause gives me the ways from the nodes. So far, so good.

My plan is then to create a python dictionary containing the node IDs as its keys, and sets of way IDs as values:

node_ways = {}

for idx, way in enumerate(result.ways):
    nodes = way.get_nodes(resolve_missing=True)
    for node in nodes:
        if node.id in node_ways:
            temp = node_ways[node.id]
            temp.add(way.id)
            node_ways[node.id] = temp
        else:
            node_ways[node.id] = set([way.id])

This works fine, but has two problems:

The use of resolve_missing=True takes around 3.5 minutes to resolve all the missing nodes for the ways, and I'd very much like to reduce this time. I have tried removing the resolve_missing=True and putting an exception handler around get_nodes, but that results in a set of nodes which is incomplete.
When I use resolve_missing=True, I get all the nodes on all the ways which intersect with the original set of nodes. This is more nodes than I need. I could reduce this down in the code, but ideally I'd rather get rid of the need to use resolve_missing=True somehow.

So, my question is: can I either adapt the original query to give me the output I want directly, or can I get what I need and avoid using resolve_missing=True by some other method?


texts:
I have some roads and some points containing elevation data.

Obviously the data does not cover the whole area, as it was created by a surveyor. Based on this data, I need to check, whether the inclination of the roads surpasses a specific threshold.

I used linear TIN-interpolation to create a DEM and based on this the inclination - which has a lot of erroneous triangles:

I tried IDW- and cubic interpolations, but these don't produce an output ("couldn't create layer").

Basically I either need an option/tool to create the DEM with fewer artifacts/errors, or to interpolate the data.


texts: MMQGIS export shapefile to KML not working, error message: 'Unknown geometry type 0'. But Point data being used?

texts:
I'm trying to convert a shapefile (created from a CSV) to a KML using the MMQGIS Plugin (as I want to show styling by attributes on Google Earth.) However, when clicking 'apply' it says 'Unknown Geometry type 0. Use convert geometry tool to change to point, line or polygon'.
The data I am using is Point data, so confused this is an error for the data?


texts: How to import data (with schema) from MSSQL to PostgreSQL?

texts:
I have an ArcSDE enterprise geodatabase in MS SQL Server 2008. I want to replace MS SQL Server with PostgreSQL.

I have setup and created the enterprise geodatabase in PostgreSQL and now want to import the data along with the schema to PostgreSQL.

I tried with export/import XML Workspace but I get the following error:

  import xml failed;
  The parameter is (or has an element that is) the wrong kind of geometry.


texts: GeoServer Connection Pool uses most available connections and not releasing them

texts:
We are using GeoServer with a geoserver cache for some layers. Most of our layer updates frequently and therefore we do not want to setup a cache for them (they do have in memory cache though).

First, we use different workspaces for different types/actions we want to perform and want to keep the different workspaces/layers separate, however they all connect to the same PostGIS database.

Now I am facing a problem where we have 15+ workspaces, each workspace has a connection pool (default 10) to the same PostGIS DB. So that means that our GeoServer instance is taking more connections that PostGIS has (100), also we do have other applications running on the same database.

The connection pool for each store has a max connections = 10, min connections=1, connection timeout = 10, evictor run periodcity = 30, max connection idle time = 30, evictor tests per run = 3.

I noticed that PostGIS has a lot of connections on idle and does not close them. What is the best way of connection GeoServer to a PostGIS database with the minimum connections? How can I force the idle connections to close without changing the max connections?


texts: Why shapefiles created with sf package (in R) do not display same CRS in ArcMap?

texts:
I have created a polygon shapefile in ArcMap using 7 points(lat,lon) then read this shapefile into R to obtain the CRS. I then created a shape file in R using the sf package (and the 7 points) and used the CRS obtained from the ArcMap shapefile. When I bring the shapefile created in R into ArcMap I get a warning regarding different coordinate systems. Why would this happen?

My steps:

7 points (lon and lat) where point 1 = point 7 (closed loop). Create a polygon shapefile in ArcMap and display. Properties in ArcMap show : Geographic Coordinate System: GCS_North_American_1983 Datum: D_North_American_1983

In R I use rgdal to read in shapefile and obtain CRS

cashes <- rgdal::readOGR(dsn="path",layer="shapefile")
cashes@proj4string
This results in : "+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0 "

Now I create the shapefile in R using the sf package (I could use %>%)

Rshapefile <- readxl::read_xlsx(path,sheet="cashes")
Rshapefile <- sf::st_as_sf(Rshapefile,coords=c("lon","lat"),crs="+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0 ")

Edit: At this point if i inspect Rshapefile i get this
proj4string:    +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
Note: +datum=NAD83 is dropped.

Rshapefile <- dplyr::summarize(Rshapefile,geometry = sf::st_combine(geometry))
Rshapefile <- sf::st_cast(Rshapefile,"POLYGON")

sf::st_write(Rshapefile,dsn=path,layer="Rshapefile",driver="ESRI Shapefile")

I then bring this Rshapefile into ArcMap and the coordinate systems are different. When I view the properties in ArcMap they show: Geographic Coordinate System: GRS 1980(IUGG, 1980) Datum: D_unknown

Why is this? I must be doing something wrong but I can't figure out what!


texts: ArcPy script using 'symbology.addAllValues()' failing

texts:
I'm having a problem trying to use symbology.addAllValues() on a layer with Unique Values Symbology. Unfortunately, I can't paste in actual code here since the machine I'm working on is classified, but essentially I am trying to do this:

source.symbology.valueField = "Name"
source.symbology.addAllValues()

it is failing with an error "in addAllValues return convertArcObjectToPythonObject ..."

The funny thing is that if I manually go into Arc and tell it to add all values, then run my script, it doesn't give the error.

Is addAllValues expecting a table of the values to already exist so it can update it, and it is choking because it doesn't find it?


texts: Where can I download Southern California Interstate Highway System in a Shapefile?

texts: Does anybody know where I can download the Interstate highway system (as a shapefile) for Southern California (the area around San Diego in particular)?


texts:
I have multiple shapefiles with different geometry and I am trying to add a label of layers.

For creating a map I am using GeoServer and OpenLayers.

Requirement is layer are visible on a max zoom level but labels are visible on min zoom level for same shapefile.

Ex : Layer gets visible at 1:6,00,000 but labels gets visible at 1:5000 for the same shapefile.

Please if anyone know the answer reply me as soon as possible.


texts: TypeError: ufunc 'expi' not supported for input type, and input can't be safely coerced to supported type according to casting rule ''safe''

texts:
Here is my code for raster calculation in ArcGIS. However, I am getting raster file for "uw" only. The python is unable to calculate exponential integral of raster uw. Thus, I am not getting the required result which I expect after calculation as raster4. Any help from the experts would be of great help to me.

import arcpy
from arcpy import env
from arcpy.sa import *
import numpy as np
from scipy import special as sp
from scipy.special import expi

##### parameters####
rw = 0.15
S = 0.003
t = 30

#### raster band selection ####
raster1 = Raster('allowable head')
raster2 = Raster('Transmissivity')

##### well function calculation ###
uw = rw**2 * S / (4 * raster2 * t)

#### calculate exponential integral of well function ######
W_uw = -sp.expi(-uw)

##### capacity calculation ####
raster4 = (raster1 * raster2 * 4 * np.pi * t) / W_uw

expi in the code refers to "exponential integral" mathematical function which is used to solve the Euler's series. The code is executing up to "uw" only. I receive a raster output of "uw" and nothing after that.

I edited the code as per the suggestion of @Hornbydd:

##### well function calculation ###
uw = rw**2 * S / (4 * raster2 * t)

#### Convert Raster to numpy array
arr = arcpy.RasterToNumPyArray(uw,nodata_to_value=0)
### expi of array
arrexpi = -sp.expi(-arr)

# Convert Array to raster (keep origin and cellsize the same as the input)
newRaster = arcpy.NumPyArrayToRaster(arrexpi,
                                     lowerleft,cellSize,value_to_nodata=0)

However, I still get the value of uw only and noting after that. How to set lowerleft value? Is it the X-cordinate and Y-cordinate? Please help.


texts: Moving cartographic representation of labels in ArcMap?

texts:
Cartographic representation in ArcMap is a nice way to temporary change geometry of features, so you can make adjustments just for map making.

In my mapping project I need to move some point features with representations and also, I need to automatically label them. It seems that labels are always positioned in relation to the original geometry, not where the representation was moved.

Is this how representations are supposed to work or is there a way to automatically label representations?


texts: Render 3D spikes based on population data in QGIS?

texts:
I believe it can be done with 3D Analyst in ArcGIS, but is it possible to achieve something like this in QGIS?

Or should I generate a raster and use it as a relief map in a 3D application?


texts:
I have a shapefile and I want to select some features and update specific field with pyqgis automate.

for example :

myshape:
myfield(values)
1
2
3
4
5
6
7
8

I want to select some features and update automate something like this (pseudocode):

newvalue = 1212

for feature in selected:
    i = feature[myfield]
    feature[myfield] = i + newvalue

i is order number to avoid duplicate.

Resulting in something like (provided 5...8 were selected):

myshape:
myfield(values)
1
2
3
4
1217
1218
1219
1220


texts: Alternative of iface for standalone PyQGIS application

texts:
To know the active layer and the attribute of the qgis project, I'm using the iface variable in the integrated Python console.
But I want to achieve this information through a standalone Python script. So how to retrieve the information (activelayer, attribute, features, etc) without the iface.
As per my understanding iface variable will work for the integrated Python console in QGIS.


texts: How do I correct the GeoTransform of a GRIB dataset if I scale/crop the data?

texts: I have modified the GRIB2 data to be cropped and a different size. I'd like to create a new GRIB2 file which contains this modified dataset. How do I calculate the new GeoTransform and create it using gdal? Do I need to modify the Projection too?


texts: Programmitcally Change Text Font for all Labels on MXD

texts:
I'd like to know if there is a way to change the text font for all layers and all label classes within a layer on a MXD.

I use a font not often found on all computers for created maps. I love the font, but it becomes an issue when publishing or sharing MXD's. Usually before I share or publish a MXD I need to go through all the layers and set the font to be something found everywhere.

Personally I'd like if the script auto found layers and label classes withing every layer in the TOC. I'm pretty good with python, but I haven't ever found a way to interact with the label's for a layer within ArcPy.

Does anyone have an idea of how to due this?


texts: Geostationary satellite projection to be introduced in GDAL

texts:
I'm working on raster files and their spatial reference system is:

PROJCS["Geostationary_Satellite",
    GEOGCS["GCS_Unknown datum based upon the custom ellipsoid",
        DATUM["D_Not specified (based on custom ellipsoid)",
            SPHEROID["Custom_ellipsoid",6378140,298.252981]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Geostationary_Satellite"],
    PARAMETER["central_meridian",0],
    PARAMETER["satellite_height",35785831],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]

How can I introduce it to GDAL? Is there proj4 or WKT equivalent to it ?


texts: How to Add a Legend in config.yaml File in GeoServer Printing

texts: I am unable to add a legend in config.yaml file in GeoServer.  Kindly, could any one share the code to show the legend in GeoServer print?


texts: Problems with QGIS plugin traceDigitize on Mac?

texts:
I'm trying to get the traceDigitize plugin to work with QGIS 2.0.1 on my Mac while following this tutorial but I get this error as soon as I ctrl-click:

An error has occured while executing Python code:

Traceback (most recent call last):
  File "/Users/oskarlin/.qgis2/python/plugins/traceDigitize/traceDigitize.py", line 135, in createFeature
    if (geom.isGeosValid()):
AttributeError: 'NoneType' object has no attribute 'isGeosValid'

Python version:
2.7.5 (default, Aug 25 2013, 00:04:04)
[GCC 4.2.1 Compatible Apple LLVM 5.0 (clang-500.0.68)]

QGIS version:
2.0.1-Dufour Dufour, f738351

Python path: ['/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/processing', '/Applications/QGIS.app/Contents/MacOS/../Resources/python', '/Users/oskarlin/.qgis2/python', '/Users/oskarlin/.qgis2/python/plugins', '/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins', '/Library/Frameworks/SQLite3.framework/Versions/B/Python/2.7', '/Library/Frameworks/GEOS.framework/Versions/3/Python/2.7', '/Library/Frameworks/GDAL.framework/Versions/1.10/Python/2.7/site-packages', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python27.zip', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-darwin', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-mac', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-mac/lib-scriptpackages', '/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-tk', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-old', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-dynload', '/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/PyObjC', '/Library/Python/2.7/site-packages', '/Applications/QGIS.app/Contents/Resources/python/plugins/fTools/tools']

Could this be a Mac specific problem or am I doing something wrong?

I'm using OSX version 10.9 (Mavericks) with QGIS 2.0.1 and GDAL Complete framework.


texts: How to link Python and FME to automate an operation of converting my IFC files to CityGML files

texts:
I am working on developing a web app for  3D data stockage, conversion and visualization .
In fact I need to develop a code that  automates the conversion from  IFC format to 3D CityGML format, in that way for every file IFC uploaded, I can easily have it in 3D CityGML format.
While searching, I was told that it is possible to achieve the aim by using both FME and Python and I wonder if you can help me knowing how to create the Conversion Python code that will help me convert all my IFC files. I mentioned that I already have done the FME conversion model and it is working (IFC to CityGML workspace).


texts: Using callback with python gdal.RasterizeLayer?

texts:
In a QGis python plugin i rasterize a vector layer with this code. All work fine except the callback: no error but nothing happen. The function callbackRasterize is never called

    gdal.RasterizeLayer(target_ds, [1], source_layer,
                        options=["ALL_TOUCHED=TRUE", attribute],
                        callback=self.callbackRasterize)

I found that using this code, the progress is visualized on standard output, but i can't use it with python:

    gdal.RasterizeLayer( .... callback=gdal.TermProgress)

How can I do to know when the rasterization process is over?


texts: I need to plot contour plot for a geodetic dataset.I have data for 94 locations of New Zealand. I have been able to generate a map a New Zealand using cartopy library of Python. I have used python's tricontour and tricontourf functions to interpolate the data to equally spaces grid points( for plotting contours). But I have been getting contours which are not smooth. Has anyone else used this function? Any other suggestions for mapping my data to grid points?


texts: Raster (.tif) files fail to load after ~15 minutes on ArcMap

texts:
When using ArcMap 10.6, the base mapping we use fails to draw after about 15 mins and also any other rasters (these are all .Tif files) fail to draw as well. It shows a black box or is just white, as shown in picture below.

I have to keep restarting the Mxd and this becomes quite a pain. Can anyone shed any light on what may be causing this?


texts: Symbolizing field with many values in Mapbox Studio?

texts:
I am using Mapbox Studio and trying to symbolize my Tileset which has about 200 unique values.

I am able to create individual 'stops' and manually assign a unique color, but looking for a more efficient way to symbolize each of these 200 values without having to touch each value manually.

What I've tried

I can easily create unique symbology in ArcGIS and QGIS, but Mapbox does not appear to accept such formats.
Searching Mapbox's documentation to see if I can apply a unique default color symbology for each unique value, but haven't found such a process.
Searched GIS Stack Exchange and StackOverflow, but haven't found such a question.


texts:
How do I add or append excel data to an existing attribute table?

The excel column headings will match the attribute table headings (in name at least). Simply want to expand the attribute table without having to manually type all the data in.


texts: Python Script to convert a mass of KMLs to a single fGDB

texts:
I have a mass of KML files, approx 350+ that I need in a single fGDB. I lifted a string of code from ESRI help to attempt to do this - I got it working to the point of creating all the individual gGDBs (one per KML) but its failing with an invalid character in a variable somewhere. Below is the code, and then the error that it throws out. I've added some print commands to show whats in the variables right before it fails. Help? Also, this is my first time trying to use python for actual work, so I apologize in advance for n00bish mistakes.

# Name: BatchKML_to_GDB.py
# Description: Converts a directory of KMLs and copies the output into a single fGDB.
#              A 2 step process: first convert the KML files, and then copy the featureclases

# Import system models
import arcpy, os

# Set workspace (where all the KMLs are)
 arcpy.env.workspace="C:\\Users\\mperez\\Desktop\\IT-Paths"

# Set local variables and location for the consolidated file geodatabase
outLocation = "V:\\2155\\active\\215510456\\GDB"
MasterGDB = r'AllKMLLayers.gdb'
MasterGDBLocation = os.path.join(outLocation, MasterGDB)

# Create the master FileGeodatabase
 arcpy.CreateFileGDB_management(outLocation, MasterGDB)

# Convert all KMZ and KML files found in the current workspace
 for kmz in arcpy.ListFiles('*.KM*'):
  print "CONVERTING: " + os.path.join(arcpy.env.workspace,kmz)
  arcpy.KMLToLayer_conversion(kmz, outLocation)

# Change the workspace to fGDB location
arcpy.env.workspace = outLocation

# Loop through all the FileGeodatabases within the workspace
wks = arcpy.ListWorkspaces('*', 'FileGDB')
# Skip the Master GDB
wks.remove(MasterGDBLocation)

for fgdb in wks:

  # Change the workspace to the current FileGeodatabase
  arcpy.env.workspace = fgdb

  # For every Featureclass inside, copy it to the Master and use the name from the original fGDB
  featureClasses = arcpy.ListFeatureClasses('*', '', 'Placemarks')
  for fc in featureClasses:
    print r"COPYING: " + fc + r" FROM: " + fgdb
    fcCopy = fgdb + os.sep + 'Placemarks' + os.sep + fc
    print r"fcCopy = " + fcCopy
    print r"MasterGDBLocation = " + MasterGDBLocation
    print r"fgdb = " + fgdb
    arcpy.FeatureClassToFeatureClass_conversion(fcCopy, MasterGDBLocation, fgdb[fgdb.rfind(os.sep)+1:-4])

# Clean up
del kmz, wks, fc, featureClasses, fgdb

Output from shell showing error:

COPYING: Points FROM: V:\2155\active\215510456\GDB\canal-move-2-golden-gate-blvd.gdb
fcCopy = V:\2155\active\215510456\GDB\canal-move-2-golden-gate-blvd.gdb\Placemarks\Points
MasterGDBLocation = V:\2155\active\215510456\GDB\AllKMLLayers.gdb
fgdb = V:\2155\active\215510456\GDB\canal-move-2-golden-gate-blvd.gdb

Traceback (most recent call last):
  File "C:\Users\mperez\Desktop\batch_kml_to_gdb.py", line 46, in <module>
    arcpy.FeatureClassToFeatureClass_conversion(fcCopy, MasterGDBLocation, fgdb[fgdb.rfind(os.sep)+1:-4])
  File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\conversion.py", line 1547, in FeatureClassToFeatureClass
    raise e
ExecuteError: ERROR 000354: The name contains invalid characters
Failed to execute (FeatureClassToFeatureClass).


texts: Where can I get detailed Bathymetry data for Ontario

texts: I need detailed Bathymetric data covering water bodies (including smaller ones) in the Province of Ontario.  Our Provincial Government has it, but the use of the data are restricted due to liability issues.  Does anyone know place to download this data, or can anyone suggest a way to interpolate lake depth contours (ie - would it be feasible to get a Google Earth Pro licence for this)?  I have access to 3D Analyst and Spatial Analyst.


texts:
we developing a project which contains a lot of data ,
using WMS layers for displaying the points,
Also we used WMS to return points information in a popup window ,
however to highlight selected point we used WFS layer.

The problem now is:
sometimes the popup displayed without highlighting the selected point and vice versa ..!

how could i combine between these tow features (highlight selected point and show the popup for it)?

notes:
i think one reason is the clicking event defined different,
WFS using:
 selectControl.events.register("featureselected", this, function(e) {
  select.addFeatures([e.feature]);}

and WMS using:
popupwindow1=new OpenLayers.Control.WMSGetFeatureInfo({
     layers: [layer1],
     autoActivate: true,
     maxFeatures: 3,
    infoFormat: "application/vnd.ogc.gml",
 eventListeners: {....}

The second problem is :
how could i return features from WFS layer without adding the vector layer to the map?


texts:
There is a point layer with points in a irregular grid. I need to separate them in 11 more or less identical cluster so that each cluster has more ore less the same amount of points in one area. Is there a smart way in QGIS or SpatiaLite?
I need more or less the same amount of points in the clusters of points. That is not possible with k-means.


texts: Changing color of columns names in query builder

texts: Does anybody know if it is possible to change the colors of the names of the columns (attributes) when writing a query in QGIS? (e.g. when selecting by expression in the attributes table).


texts: Finding "negative contours" when making contour lines from DEM

texts:
I have a large 2m DEM which consists mostly of mountains. But here and there there are depressions in the terrain which I need to find, or even better, automatically find in contour lines I'm planning to do.

Here's a sample image of contours from a small a part of the region, I have (manually) marked the depression contours as yellow.

In my final map I'm going to give the depression another style than the other contour lines so that the person using this map doesn't think it might be a hill!

Now I hope this is possible to achieve automatically in one way or another. The contour lines are generated using gdal_contour and has an ELEV attribute field.


texts: Choosing colour-ramp to use for elevation?

texts:
This is maybe rather subjective as it concerns style, but I would like to have a few decent colour ramps for my DEM when it is layed over my hill-shading.

I have created the Swiss hillshade model as shown here, but I'm not too convinced that the colour scheme is good.

Any examples or screenshots would be great.


texts: Why moving a label in QGIS sets the label rotation to 0° and how to avoid that?

texts:
I have some line labels that are parallel to the line. When I manually move a label with the "moove label" tool, it sets the rotation of the label back to 0 Degrees. It is not parallel to the line anymore. Why is QGIS doing this and how do I stop this?

I'm using QGIS 3.4.13


texts: Editing with Map Topology performance in ArcGIS Pro?

texts:
I have been working on some basic cartography for a fantasy world I am creating, and ArcMap 10.4 was being slow on editing some of the large island polygons I am creating. I figured that as all I am doing is basic editing I would try ArcGIS Pro for speed and to be able to use my 16 gb of RAM. I love how responsive the program is generally, but editing is frustrating. When I click on "Move" to move a feature, for example, the program doesn't allow me to actually move or rotate any feature until several minutes of the mouse cursor with the spinning wheel. All I have is a basic polygon layer for islands and a basemap. Why is it being so slow?

It's in a File Geodatabase. They do have large amounts of vertices because several of the polygons are the size of large islands several hundreds of miles across. I turned off Map Topology and things seem to have sped up a lot -- I bet that's the problem!


texts: Where has the 'Set to map canvas Extent' on QGIS Version gone?

texts:
I am using QGIS Version 3.10.0 . When using Map Composer I can not find the 'Set to map canvas Extent button'. When clicking on the extent tab The image below shows all the options I have.

Does anyone know why it is not there anymore?


texts: Specifying a file extension for the file output using a script from the processing toolbox in QGIS

texts:
I'm looking to specify the file extension of the output file created by running my script from the processing toolbox. Currently, the user must type in the extension in themselves, which I would like to avoid if possible. Ideally, there would be a drop-down, and among the options I would like to have the .pdf extension type.

My goal is to increase the usability of the script. It's not a major issue to  have the user type in .some extension manually for their file output, but I believe a drop-down would be easier to use than having to type in a file extension each time.

CODE

To create the output option where the user has to choose a path for their output file, I use the following code in the script:

##PDF_file=output file

This code works as long as the user specifies and manually types in the .pdf extension illustrated below.

Choosing output file location


texts:
I am using QGIS 2.18.16 on Windows 10.

What I want is a label to appear within the boundaries of a polygon.

I select a layer, then select command Layer>Labelling, and I get:

From here, I can change "No Labels" to "Show labels for this layer".

But I can't actually create label text.

I find a clue in a prior question:

“Labelling polygon using QGIS?”

The answer was:  The layer you are trying to label must have at least one attribute field populated with data. You can then turn on labelling in Layer Properties and select which attribute field to use for your labels.

I can enter command Layer>Open Attribute Table, and what I get is:

If I need something in the attribute table in order to apply a label, how do I accomplish that?


texts: Detecting click outside map control using OpenLayers 2?

texts:
My problem is the following:

User loads form containing map control and then clicks the edit icon and  starts editing a shape on the map.
The shape is in edit mode as you can see the handles.
Unfortunately when the user clicks the Save button the map updates are not captured.
The current work-around is to click in a blank area outside the shape. This causes the handles to disappear and the Save operation succeeds.

My question: what event should I register in the map object that would detect change of focus. That is, moving out of the map control.

If could handle this event I thought I could deactivate the map (remove the handles by turning off the edit mode) and therefore allow the map to be saved correctly.

Here's a picture illustrating my problem:


texts: ogr2ogr append 1 row to a File GDB feature class's data table

texts:
I export a polygon layer from Oracle to a FileGDB using ogr2ogr. I need to append one row to the data table, with '0' for numeric fields and 'No Elements' for the one string field.  In ArcGIS, I run a model that appends a feature class with the one observation to the target feature class (the one exported from Oracle).

How can I perform the same function using ogr2ogr?  Would I append the 'no_element' feature class to the target feature class, and if so, what is the ogr2ogr syntax?  Both feature classes are in the same File GDB.

Or, is there something I could add to my ogr2ogr syntax when I export my data from Oracle; I've included my ogr2ogr for my Oracle export.

     ogr2ogr -update -overwrite -nlt polygon -nln EO_DSMRE  -a_srs EPSG:3089 -progress -f FileGDB G:/Office/Heritage/DataMngt/dpaml_app/KSNPC_Update/DSMRE.gdb OCI:biotics_dlink/XXXXXXXXXXXX@biotics5:EO_SHAPE -sql "select  EO_SHAPE.SHAPE SHAPE,  EO_SHAPE.EO_ID EO_ID, to_char( EO_SHAPE.EO_ID) EO_ID_STR, CASE WHEN (DSMRE_EO_VW.ELCODE LIKE 'P%' OR DSMRE_EO_VW.ELCODE LIKE 'N%') THEN 'Plant' WHEN (DSMRE_EO_VW.HIGHER_TAX_CODE LIKE 'IMGASF%' OR DSMRE_EO_VW.ELCODE LIKE 'IMBIV%' OR DSMRE_EO_VW.ELCODE LIKE 'AF%') THEN 'Aquatic Animal' WHEN DSMRE_EO_VW.ELCODE LIKE 'AMACC%' THEN 'Bat' WHEN DSMRE_EO_VW.ELCODE LIKE 'C%' THEN 'Community' ELSE 'Terrestrial Animal' END TAX_TYPE  from  EO_SHAPE,  DSMRE_EO_VW where  DSMRE_EO_VW.EO_ID =  EO_SHAPE.EO_ID"

***Update

I have created a table in the Oracle data base with the fields 'SHAPE, EO_ID, TAX_TYPE'.  The values for these fields are 'NULL, 0, "No Elements"'.  I also created a view of the data that I have been exporting.  Then I tried this:

    select no_element_table.* from no_element_table union select shape_table.* from shape_table

That did not work, SQL Developer spits out ORA-22901: cannot compare VARRAY or LOB attributes of an object type.

So at this point, now that I have created the 'No_Elements' table, I just -update -append that to the target table.

    ogr2ogr -update -append ......


texts: Classify Data by Method Regardless of Data Range

texts:
I want to apply the same classification method to multiple Rasters in ModelBuilder Reclassify tool. The task is to apply the same classification method (percentile splits) across multiple rasters without setting the data value splits manually as I will always need the same 20/20/20/20/20 split of data.

How do you input a pre-defined classification method automatically e.g. Method: Quantile with 5 classes?


texts: Can a Relate within an MXD be detected programmatically using ArcPy?

texts:
I would like to know if it is possible detect relates in an MXD using arcpy.

Similar problem as the question below, but for relates instead of joins.

Detecting Join programmatically using ArcPy?

More specifically, relates listed in the "Joins & Relates" tab of the "Layer Properties" and not geodatabase "Relationship Classes".

I have checked the ArcGIS documentation and have done several Google searches without any luck.


texts: Change CRS of large GeoTIFF and extract windows from large raster in Python

texts:
This is the first time I am working with GeoTIFF formats and I want to do this in Python. I am trying to use rasterio and pyproj for the following tasks.

I have already made some first steps based on numerous Stackexchange questions. Now I am stuck at a point, where I am lost - too many suggestions how to approach this - which appear to be close to my problem, but not close enough. I could possibly copy-paste more together, but I am afraid that I am not approaching a clean solution.

Therefore I share a jupyter notebook via GitHub where you can see what I have already done and where I formulate the questions right where they come up.
https://github.com/bebissig/geotiff_basic_manipulations/blob/master/copernicus_preprocess.ipynb

Overview:
I have a GeoTIFF (~5GB) from the Copernicus database (https://land.copernicus.eu/imagery-in-situ/eu-dem/eu-dem-v1.1?tab=download) that contains elevation data of central Europe and need to achieve the following goals.

Goal/Issues:

(These are the main goals, specific questions are in the notebook.)

First; The data comes in EPSG:3035 projection. I need to extract data from Switzerland in projection EPSG:21781. So in my thinking I would be "nice" to translate the whole raster to EPSG:21781. How to do that, do I only have to change the .crs accordingly? Or also manually change the .transform? Or would you suggest to approach this totally differently?
Second; Say I have coordinates of a point in the destination projection (EPSG:21781). I want to extract a window around such a given point say coords+-2000m (in vertical and horizontal direction). Ideally, this would be still a GeoTiff, such that it still contains x,y coordinates in that window (in destination projeciton EPSG:21781). (This is also in view that possibly, later I need to add other data-bands from another database so, I guess I want to stick to a GeoFormat as long as possible).

That's actually it. If I had these windows I could read the data to np-arrays and go on with processing.

--> Se my answer below for a snippet how I made it work.


texts:
I have some autoCAD data in a locally made projection so no epsg code, no wkt..

I need to move the coordinates/.dxf data from the local projection to the national grid projection ( EPSG:3067 ETRS-TM35FIN )

I have been told that the creator of the local projection uses a helmert transform to move the data.

" We are using a simple Helmert transform for the conversion from the old local to the national system:

A-parameter : 0.996176526906 (Rotation and scaling)

B-parameter : 0.091026561878 (Rotation and scaling)

C-Parameter: 0.19725000 (z)

Origin (local x,y,z): 87581.87975000,70988.64475000,-8.16575000

New location (ETRS x,y,z): 6693404.95075000,239228.43500000,-7.96850000"

Is there any way i can do this in Qgis to make the old data in the old projection match up with the data in the new projection - EPSG:3067 ?


texts:
How do i change the map interface in arcgis online and portal? I have been trying to find a clear answer online for about a week and I have yet to find anything helpful from ESRI or any forum.

Basically I wanna go from this kind of interfance

To this

I dont have a lot of experience with arcgis Portal so please bear with me!


texts:
I'm trying to join attributes by location - my Join layer is one I created from a 10m fixed width buffer, my geometric predicate is only 'within'.  When I run the algorithm, I get an error of "-2147483646 See log for more details".  If I reverse the parameters and use my Join layer as the target and use 'contains' instead of 'within', the algorithm runs as expected (but then I'm not getting the values on the layer I want...)

I've read other articles on here and tried the suggestions (check geometries) but with the same results.


texts: Something changed in QgsVectorDataProvider().changeAttributeValues() in QGIS 2.8?

texts:
I have a plugin that works with QGIS 2.2 but does not work with 2.8.
Only deleteFeatures() works fine, but changeAttributeValues(), changeGeometryValues() and addFeatures() not.

I receive r == False on the code below.
The layer is an SpatiaLite table.

for feat_vias in lay_vias.getFeatures(req):

    if  feat_vias.attributes()[idx_pkvia] == div_rowid:
        fid = feat_vias.id()
        attrs = { idx_setvia : s_geocodi }
        r = lay_vias.dataProvider().changeAttributeValues({ fid : attrs })

I'm using QGIS 2.8.1 - Windows

In another question (Speed of editing attributes in QGIS from a Python plugin) I found an Matthias Kuhn's answer that solve this question too.

I put the commit sentence outside the loop as below. Work fine!

for feat_vias in lay_vias.getFeatures(req):

    if  feat_vias.attributes()[idx_pkvia] == div_rowid:
        fid = feat_vias.id()
        attrs[fid] = { idx_setvia : s_geocodi }

r = lay_vias.dataProvider().changeAttributeValues(attrs)


texts: Using OpenLayers Selectfeature get mouse position?

texts: How do I get mouse coordinates after a SelectFeature, which returns the feature by itself?


texts: How can I access individual layer in Geoserver similar to the arcgis rest-end urls?

texts:
I have got Geoserver installed with map data in postgis.

When I use arcgis server to publish mapService and featureService, I can get a list of layers published under the url like following.

https://"ARCGIS-SERVER".com/arcgis/rest/services/"servcie folder"/"servcie folder"_gis_fs/FeatureServer

Further, I can access individual layers within the service and can see the list of attribute fields and data type information within the data.

Is it possible to get a similar information from geoserver? if yes, How?


texts: Postgre/PostGIS - How to find points inside an area

texts:
I have a database of clients/users with their coordinates.
I would like to select all clients within the marked area.

Query:

SELECT * FROM clients
INNER JOIN polygons
ON polygons.polygon @> clients.geo_point

works fine for areas where there is not empty area inside.
If there is en empty space inside like on the picture, blank area is selected as well.

If the area looks like this:

Everything works fine.

Which query and which spatial function should I use to select clients correctly?


texts:
I have created mapbooks using DDP and am running into an exporting issue. I have 40 mapbooks of 40 different cities based on a my layer file called "Custom Grid".  When I go to export it gives me the option to either export all of the mapbooks together as one single PDF, or every page of each mapbook as an individual pdf.

There is a field in my layer file called "City_ID" in which I wish to append each mapbook page by.

I was wondering if there was a way to export an individual PDF for each mapbook (each city) without having to manually select each city (select by location or attributed) and then export one by one. Is there a simple arcpy script for this?


texts:
I have a problem with lowess function to smooth NDVI as a time in R. Here is the scripts I run:

#Create a time vector to act as x
time<-sort(sample(97:297, nlayers(ndvi2000crop1)))

#locally-weighted polynomial regression of ndvi ~ time
t.lowess <- function(x,...) { if (is.na(x[1])) {NA} else {lowess(x,y,...)$y } }
f.pred <- calc(ndvi2000crop1, t.lowess)

When I calculate the function on image, an error appeared something like this "Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply): cannot use this function". Can you tell me where this error comes from in this part?


texts: Displaying attribute value from related layer in QGIS Atlas?

texts:
I work on QGIS 3.4.4 Madeira. I have two related layers in my project :"towns" (polygon) and "interest points" (points), the common field is "town name".

In the Atlas, my cover layer is "towns".

I would like to put labels wich display values from different fields of the child layer ("interest points") corresponding to the current entity of the Atlas. It seems easy but I can't find the correct expression.

My cover layer in the Atlas is "towns" and in each town, I have several "interest points". The layer "interest points" have many fields and I would like to insert data from a given field of the child layer using an expression in the labelling (text) tool and not the table tool. So I need an expression to insert the value inside a sentence for example.

I tried relation-aggregate function, it returns values from child layer but it necessarily "calculates" (min, max, sum…). And I only want to display the child related value.


texts: MapInfo - how to *publish* a WFS service?

texts:
I seem to have my Pitney & Bowes technologies in a scramble.  We installed the demo version of MapInfo hoping to build a map and publish it as a WFS service.  We built the map but found no publish options.  We next tried MapXtreme and setup a 'Workspace' but there are again, no publish options.  I was expecting it to register assemblies in the GAC (e.g. MapInfo.CoreEngine) but it did not.

QUESTION: What P&B technology is used to publish a WFS service?

ANSWER: See page 400 on this PDF


texts:
Is there an alternative program to ogr2ogr when trying to convert file formats?
I've played around with ogr2ogr quite a bit lately (on Windows) but I can't get around some issues ..
For example when converting shapefile to geojson it destroys the attribute table due to the wrong encoding format. Specifying it won't help either because this feature doesn't seem to work on Windows machines ..

So I'm getting tired of this program and want to use something different. And I don't want to use QGIS, because it takes a lot of time to start and I don't want to wait 5 minutes to convert just a single file ..

Do you know any alternatives? It needs to run offline, so no web services.


texts:
When I try to include the "fwp layer" from geoserver 2.4 (Ubuntu 12.4 VServer) into my code it doesn't show up. But the geoserver layer preview (OL) shows the right image.
The GeoExplorer preview shows only pink tiles with a questionnmark.

The fwp layer is a Geotiff file which I had transformed with
gdaladdo -r average FWP_neun1.tif 2 4 8 16

Here is the actually website: WEB GIS where it should work.

Thanks a lot for your help!


texts: Hatching over simple fill in ArcGIS 10.3

texts: I have a polygon layer that is symbolized into four categories. Two of these categories are above a threshold limit and two are below. I would like to further distinguish the two categories above the threshold limit but putting "10% simple hatch" OVER the polgyon fill colour from the colour ramp. How do I combine symbol layers so that a solid background colour appears under hatching?


texts: how to find that two polygon boundaries are connected?

texts:
I am working on a Geo database where there is a layer with many polygons. i want to find that which polygons boundaries are connected so that i can create new layer from that polygons and which polygon boundary is not connected with first polygon, like there is space between them, so i want to not include that polygon in layer. This process will be achieved through python programming

I am using python in Arcgis 10.2.


texts: Leaflet: A way to show information of overlapping polygon of the same layer?

texts: I have a polygon layer with overlapping polygons. When clicking on a spot with more than one polygon eg. overlapping polygons, the popup shows the property of only one feature (the topmost?). Is there a way to show the content of the same attribute of several objects in one popup?


texts:
I am creating an animated time series map using Carto VL and have included a progress bar that updates as the animation progresses. The bar updates, but is not interactive — for instance, a viewer would not be able to use the bar to advance or move backward in the animation. What could I change to make the progress bar interactive?

Here is the code I used to build the progress bar:

Progress: <input type="range" id="js-progress-range" min="-21000" max="0" step="500"></p>

...and here is how I update the progress bar using my animation:

function updateProgress () {
            $progressRange.value = alnusViz.variables.animation.getProgressValue();
}
  setInterval(updateProgress, 100);


texts: Are nested parenthesis allowed in MultiPoint WKT?

texts:
I'm writing a basic parser for common WKT.

The Wikipedia entry for Well-known text says that both of these are equivalent.

MULTIPOINT (10 40, 40 30, 20 20, 30 10)
MULTIPOINT ((10 40), (40 30), (20 20), (30 10))

The first one is standard, but I can't find official documentation stating that the nested parenthesis is actually valid.

Does anybody know if the second version is valid?


texts: How to use Python modules as data providers in QGIS?

texts:
I don't know if it has been discussed in the past. Qgis data providers can only load native compiled libraries (.so, .dll, etc) that offer the required interface.
Is there anything blocking to the possibility of having python data providers?
Would it be possible to have them loaded along the native ones, or dynamically through some API method?


texts: OpenLayers3 with google map local tiles server

texts:
I've got a problem with openlayers3 and my local tiles server.
I have downloaded some offline google map tiles from zoomLevel 0 to zoomLevel 16. All of the tiles are 256 * 256. And they are distributed in some local folders, which name from "0" to "16".
Then i copied these tiles into my local apache server path, and i can access these tiles with my web browser such as google Chrome,and the url is "http://localhost:8089/Maps/15/26933/12364.png".
I created a html with openlayers3 , and tried to access the tiles on my apache server. The content in html is below:

    var TileLayer = function (options) {
        var layer = new ol.layer.Tile({
            extent: ol.proj.transformExtent(options.mapExtent, options.fromProject,

options.toProject),
            source: new ol.source.XYZ({
                attributions: [options.attribution],
                url: options.url,
                tilePixelRatio: options.tilePixelRatio, // THIS IS IMPORTANT
                minZoom: options.mapMinZoom,
                maxZoom: options.mapMaxZoom
            })
        });
        return layer;
    }

    var defaults = {
        url: 'http://localhost:8089/Maps/{z}/{x}/{y}.png',
        mapExtent: [115.89752197265625,39.63319206567459,116.9000244140625,40.34445080136368],
        mapMinZoom: 1,
        mapMaxZoom: 16,
        attribution: new ol.Attribution({
            html: 'Tiles © GoogleMap'
        }),
        tilePixelRatio: 1,
    // all of the project can not work
        fromProject: "EPSG:102100",
        //fromProject: "EPSG:900913",
        toProject: "EPSG:3857"
    };
    var layer = new TileLayer(defaults);

    var map = new ol.Map({
     layers: [layer],
     target: 'map',
     view: new ol.View({
     center: [115.886536, 39.842286],
     zoom: 8
     })
     });

When i browse to this html the following error happens:
Failed to load resource: the server responded with a status of 404 (Not Found)
http://localhost:8089/Maps/10/509/512.png Failed to load resource: the server responded with a status of 404 (Not Found).

I have searched for a while, but i can't find some useful information about my problems.
Any help?


texts: How to display OSM relations (in QGIS) from osm2pgsql populated PostGIS DB

texts:
I have a PostGIS DB populated with OSM data (osm2pgsql, standart.style). Now I would like to display a hiking trail (http://www.openstreetmap.org/relation/51855).
It would also be fine, if I just could extract/export the relation as a multiline.

Is that possible?


texts: Extracting atmosphere opacity using Google Earth Engine?

texts:
I met some doubts when I try to extract Landsat surface reflection product's atmosphere opacity.

the code  is :

var inCollection7 = ee.ImageCollection("LANDSAT/LE07/C01/T1_SR").filter(ee.Filter.calendarRange(2016,2018,'year')).filter(ee.Filter.calendarRange(6,8,'month')).filterBounds(geometry);

var list=inCollection7.toList(10)

var test_img7= ee.Image(list.get(0))

print('test_img7',test_img7)

Map.addLayer(test_img7.select('sr_atmos_opacity'),{},'Landsat7_Opacity')

the question is :

why the Landsat7's atmosphere opacity' value is not from 0 to 0.3 or more（like the data description）？

[![error][2]][2]


texts: QGIS How to link Symbology to Named Style in Library

texts:
Is there any way, how to 'link' symbology to the named style in Library instead of only 'copying' the symbol definition?
Example:
 - I'm styling roads, in Symbol Selector I pick a style from Library with name i.e. 'HIGHWAY', the symbol specification is copied accordingly, ok
 - Then I open Settings->Style Manager and change the line style 'HIGHWAY' i.e from black to red
 - all the symbology created from style 'HIGHWAY' just stay black instead of changing to red
 - I want them to change if I change the 'master style HIGHWAY' in the library

Any plugin creating the 'symbology link'please?


texts: Date output does not show "time" in ModelBuilder?

texts:
I have developed a Python script using ArcPy with one of the output variable being of Date type. In the actual Python code, the variable is a Python datetime object which has both the date and time sections.

A snippet of Python script "ExtractNetCDFToRaster.py" is as below. The purpose of the code is to extract a raster layer out of a netCDF file with the raster and the date and time (in the "time" dimension) as output.

import sys, os
from os.path import basename
import arcpy
import datetime
from dateutil import tz

# Convert UTC datetime to local datetime
def convertToLocalDT(string):
    from_zone = tz.tzutc()
    to_zone = tz.tzlocal()

    # 20/11/2017 1:00:00 PM
    dtUTC = datetime.datetime.strptime(string, "%d/%m/%Y %I:%M:%S %p")
    dtLocal = dtUTC.replace(tzinfo=from_zone).astimezone(to_zone)

    return dtLocal

...
# get dimension values from the time dimension
dimension_values = nc_FP.getDimensionValue("time", 0)

# the dimension_values is of a unicode type of data, e.g. u'20/11/2017 1:00:00 PM'
strDimensionValues = dimension_values.encode('utf-8')
arcpy.AddMessage("dimension_values: " + strDimensionValues)
dtDimensionValues = convertToLocalDT(strDimensionValues)

...

# Set dtDimensionValues as output which has an alias of "outDateTime" in the tool Properties window.
arcpy.SetParameter(5, dtDimensionValues)

When I imported this script into a Model as a custom tool and ran the tool, the variable only output the date but not the time section. For example, if the actual datetime is "2017-11-21 00:00:00" the outDateTime only shows "21/11/2017".


texts:
I have 235 raster tile 20km x 20km. I need to create a grid to overlay across them...ie the grid itself will be 20km x 20km, with the lines of the grid on the boundaries of the tiles.

I am assuming the vector grid tool is the key as per this question...How to emulate a specific vector grid in QGIS?

However they have created the base map as vector..mine are rasters.

i have search for the tutorial but cannot find it.


texts:
I have small problem with upgrading OpenGeo Suite from 4.0.2 to 4.1
Everything went smooth after updating. Exept dashboard.
I had OpenGeo's Geoserver installed on port 8180. And previous dashboard was starting also on 8180 (http://localhost:8180).
Latest version installed on port 8080. I changed port as before (C:\Program Files (x86)\Boundless\OpenGeo\jetty\start.ini). But my new dashboard is starting on port 8080. And previous one one port 8180.

My question is, how to "delete" previous dashboard and change new one's port?
Thanks in advance for answer.


texts: Classification of buildings using USGS Landsat 8 images on Google Earth Engine

texts:
I am trying to classify built-up and non-built-up areas in the district of Ambala in India within the time period 2014-17 using a simple classifier model I had created.

The model uses 20,000 points from this fusion table. The classifier takes these bands: ['B1','B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B10', 'B11','NDBI','NDVI'], trains a CART model on 3000 points (1200 built up+1800 non built up). I had also tried to account for cloud cover by sorting and filtering according to my value of CLOUD_COVER. In the final result I get a mask s.t white is for built up areas and black for non-built up. The problem is there is a clear discrepancy in my results for 4 consecutive years, 2014-2017. Here are my classification results for the 4 years.

Now , the white portion unnaturally increases from 2014-15, and then subsequently decreases from 2015-2016 and 2017. Such kind of pattern in building areas is not possible, and the weird part is that I tried the same code for different districts and the pattern of discrepancy is the same among all. I wish to know what could be the cause of this issue? Is it some kind of problem with my code or the LANDSAT 8 images for 2015? And how can I solve this?

Here is my code (for reference):

//Loading India image, the extracting data for Haryana (a state in India) and then subsequently Ambala (a district in Haryana)
var bands = ['B1','B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B10', 'B11'];
var india = ee.FeatureCollection('ft:1UDdgOCf8DoRJ9bVm-UVbR6CqxtkJToLQjTFd0r0Z','geometry').filter(ee.Filter.eq('Name','India')).geometry();
var india_image = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA').filterBounds(india).filterDate('2014-02-01','2014-09-01').sort('CLOUD_COVER').limit(500).mosaic();
var fc = ee.FeatureCollection('ft:1PA2zwArj8EsplrX9eMxJ2H_TICyyx855KPnbJhC1','geometry')
var district = fc.filter(ee.Filter.eq('name','Ambala'));
var haryana_image = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA').filterBounds(district).filterDate('2014-02-01','2014-04-01').sort('CLOUD_COVER').limit(4).mosaic();
var input = haryana_image;
input = addBands(input.select(bands));
india_image = addBands(india_image);

//Loading the points from the fusion table and training the classifier
var ft = ee.FeatureCollection('ft:1fWY4IyYiV-BA5HsAKi2V9LdoQgsbFtKK2BoQiHb0');
var ft_builtup = ft.filter(ee.Filter.eq('class',1)).limit(1200);
var ft_nonbuiltup = ft.filter(ee.Filter.eq('class',2)).limit(1800);
ft = ft_builtup.merge(ft_nonbuiltup);
var new_bands = ['B1','B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B10', 'B11','NDBI','NDVI'];
function addBands(image){
  var ndvi = image.normalizedDifference(['B4', 'B3']).rename('NDVI');
  var ndbi = image.normalizedDifference(['B5', 'B4']).rename('NDBI');
  var ndwi = image.normalizedDifference(['B6', 'B6']).rename('NDWI');
  return image.addBands(ndvi).addBands(ndbi).addBands(ndwi);
}

// Load a Landsat 8 image to be used for prediction.
var training = india_image.sampleRegions(ft,['class'],30);
var trained = ee.Classifier.cart().train(training, 'class', new_bands);

input = input.clip(district);
input = input.classify(trained);
input = input.expression('LC==1?1:0',{'LC':input.select('classification')});
Map.addLayer(input.clip(district));


texts: Need to generate ids to connect points and polygons

texts:
I have a simple question as I am just getting my head into GIS.

I am trying to compare the demographics of charter schools to traditional public schools. I have universe listings (lat,long) of US public schools in a vector layer with additional layers that include separate listings of

independent charter schools (not part of a school district) and
district linked (tabular district/school links) charter schools.

I want to compare the charters to comparable local schools, but this requires exploring multiple matches- district public schools to charter schools in the same districts and  then compare schools in a 10 mile buffer of each charter school. Since the independent charters are not linked to a district I would like to use QGIS to connect points to the polygons they are in.

What I cannot figure out is, after I have the schools in a shapefile, how can I generate new ID's to nest 1) charters within the district geographic areas (polygons) and a second id for each charter school that includes all the traditional public schools within the 10 m buffer?

Any ideas?
Thanks in advance.


texts: How to generate automatically the name of output parameters in a python script tool?

texts:
I have some Python script tools in ArcGIS for Desktop.

I need to reproduce what all geoprocessing tools do when the first input parameter is introduced: autogenerate a workspace and basename for the outputs.

How can I do that?